! MODULE RKSUITE ! CONTAINS SUBROUTINE SETUP(NEQ,TSTART,YSTART,TEND,TOL,THRES,METHOD,TASK, & ERRASS,HSTART,WORK,LENWRK,MESAGE) C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C If you are not familiar with the code SETUP and how it is used in C conjunction with UT or CT to solve initial value problems, you should study C the document file rksuite.doc carefully before attempting to use the code. C The following "Brief Reminder" is intended only to remind you of the C meaning, type, and size requirements of the arguments. C C The environmental parameters OUTCH, MCHEPS, and DWARF are used in the C following description. To find out their values C C CALL ENVIRN(OUTCH,MCHEPS,DWARF) C C INPUT VARIABLES C C NEQ - INTEGER C The number of differential equations in the system. C Constraint: NEQ >= 1 C TSTART - DOUBLE PRECISION C The initial value of the independent variable. C YSTART(*) - DOUBLE PRECISION array of length NEQ C The vector of initial values of the solution components. C TEND - DOUBLE PRECISION C The integration proceeds from TSTART in the direction of C TEND. You cannot go past TEND. C Constraint: TEND must be clearly distinguishable from TSTART C in the precision available. C TOL - DOUBLE PRECISION C The relative error tolerance. C Constraint: 0.01D0 >= TOL >= 10*MCHEPS C THRES(*) - DOUBLE PRECISION array of length NEQ C THRES(L) is the threshold for the Ith solution component. C Constraint: THRES(L) >= SQRT(DWARF) C METHOD - INTEGER C Specifies which Runge-Kutta pair is to be used. C = 1 - use the (2,3) pair C = 2 - use the (4,5) pair C = 3 - use the (7,8) pair C TASK - CHARACTER*(*) C Only the first character of TASK is significant. C TASK(1:1) = `U' or `u' - UT is to be used C = `C' or `c' - CT is to be used C Constraint: TASK(1:1) = `U'or `u' or`C' or `c' C ERRASS - LOGICAL C = .FALSE. - do not attempt to assess the true error. C = .TRUE. - assess the true error. Costs roughly twice C as much as the integration with METHODs 2 and C 3, and three times with METHOD = 1. C HSTART - DOUBLE PRECISION C 0.0D0 - select automatically the first step size. C non-zero - try HSTART for the first step. C C WORKSPACE C C WORK(*) - DOUBLE PRECISION array of length LENWRK C Do not alter the contents of this array after calling SETUP. C C INPUT VARIABLES C C LENWRK - INTEGER C Length of WORK(*): How big LENWRK must be depends C on the task and how it is to be solved. C C LENWRK = 32*NEQ is sufficient for all cases. C C If storage is a problem, the least storage possible C in the various cases is: C C If TASK = `U' or `u', then C if ERRASS = .FALSE. and C METHOD = 1, LENWRK must be at least 10*NEQ C = 2 20*NEQ C = 3 16*NEQ C if ERRASS = .TRUE. and C METHOD = 1, LENWRK must be at least 15*NEQ C = 2 32*NEQ C = 3 21*NEQ C C If TASK = `C' or `c', then C if ERRASS = .FALSE. and C METHOD = 1, LENWRK must be at least 10*NEQ C = 2 14*NEQ C = 3 16*NEQ C if ERRASS = .TRUE. and C METHOD = 1, LENWRK must be at least 15*NEQ C = 2 26*NEQ C = 3 21*NEQ C C Warning: To exploit the interpolation capability C of METHODs 1 and 2, you have to call INTRP. This C subroutine requires working storage in addition to C that specified here. C C MESAGE - LOGICAL C Specifies whether you want informative messages written to C the standard output channel OUTCH. C = .TRUE. - provide messages C = .FALSE. - do not provide messages C C In the event of a "catastrophic" failure to call SETUP correctly, the C nature of the catastrophe is reported on the standard output channel, C regardless of the value of MESAGE. Unless special provision was made C in advance (see rksuite.doc), the computation then comes to a STOP. C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C .. Scalar Arguments .. DOUBLE PRECISION HSTART, TEND, TOL, TSTART INTEGER LENWRK, METHOD, NEQ LOGICAL ERRASS, MESAGE CHARACTER*(*) TASK C .. Array Arguments .. DOUBLE PRECISION THRES(*), WORK(*), YSTART(*) C .. Common Block for Problem Definition .. DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR INTEGER NEQN COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN SAVE /RKCOM1/ C .. Common Block to hold Problem Status .. DOUBLE PRECISION T, H, TOLD, HOLD INTEGER NFCN, SVNFCN, OKSTP, FLSTP LOGICAL FIRST, LAST COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP, & FIRST, LAST SAVE /RKCOM2/ C .. Common Block for General Workspace Pointers .. INTEGER PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP, & PRSTGS, PRINTP, LNINTP COMMON /RKCOM3/ PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP, & PRSTGS, PRINTP, LNINTP SAVE /RKCOM3/ C .. Common Block to hold Formula Definitions .. DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6), & E(7) INTEGER PTR(13), NSTAGE, METHD, MINTP LOGICAL INTP COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD, & MINTP, INTP SAVE /RKCOM4/ C .. Common Block to hold Formula Characterisitcs .. DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4 INTEGER ORDER, LSTSTG, MAXTRY, NSEC LOGICAL FSAL COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY, & NSEC, FSAL SAVE /RKCOM5/ C .. Common Block for Global Error Assessment .. DOUBLE PRECISION MAXERR, LOCMAX INTEGER GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR, & PRZYNU LOGICAL ERASON, ERASFL COMMON /RKCOM6/ MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP, & PRZERS, PRZERR, PRZYNU, ERASON, ERASFL SAVE /RKCOM6/ C .. Common Block for Environment Parameters .. DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY INTEGER OUTCH COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY, & OUTCH SAVE /RKCOM7/ C .. Common Block for Integrator Options .. LOGICAL MSG, UTASK COMMON /RKCOM8/ MSG, UTASK SAVE /RKCOM8/ C .. Common Block for Error Message .. CHARACTER*80 REC(10) COMMON /RKCOM9/ REC SAVE /RKCOM9/ C .. Parameters .. CHARACTER*6 SRNAME PARAMETER (SRNAME='SETUP') INTEGER MINUS1 LOGICAL TELL PARAMETER (MINUS1=-1,TELL=.FALSE.) DOUBLE PRECISION ONE, ZERO, PT01 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0,PT01=0.01D+0) C .. Local Scalars .. DOUBLE PRECISION HMIN INTEGER FLAG, FREEPR, IER, L, LINTPL, LREQ, NREC, VECSTG LOGICAL LEGALT, REQSTG CHARACTER TASK1 C .. External Subroutines .. EXTERNAL CONST, MCONST, RKMSG, RKSIT C .. Intrinsic Functions .. INTRINSIC ABS, MAX, SIGN C .. Executable Statements .. C C Clear previous flag values of subprograms in the suite. C IER = MINUS1 CALL RKSIT(TELL,SRNAME,IER) C IER = 1 NREC = 0 C C Fetch output channel and machine constants; initialise common C block /RKCOM7/ C CALL MCONST(METHOD) C C Check for valid input of trivial arguments TASK1 = TASK(1:1) LEGALT = TASK1 .EQ. 'U' .OR. TASK1 .EQ. 'u' .OR. & TASK1 .EQ. 'C' .OR. TASK1 .EQ. 'c' IF (.NOT.LEGALT) THEN IER = 911 NREC = 3 WRITE (REC,'(A/A,A,A/A)') &' ** You have set the first character of ', &' ** TASK to be ''',TASK1,'''. It must be one of ', &' ** ''U'',''u'',''C'' or ''c''.' ELSE IF (NEQ.LT.1) THEN IER = 911 NREC = 1 WRITE (REC,'(A,I6,A)') &' ** You have set NEQ = ',NEQ,' which is less than 1.' ELSE IF (METHOD.LT.1 .OR. METHOD.GT.3) THEN IER = 911 NREC = 1 WRITE (REC,'(A,I6,A)') &' ** You have set METHOD = ',METHOD,' which is not 1, 2, or 3.' ELSE IF (TSTART.EQ.TEND) THEN IER = 911 NREC = 1 WRITE (REC,'(A,D13.5,A)') &' ** You have set TSTART = TEND = ',TSTART,'.' ELSE IF ((TOL.GT.PT01) .OR. (TOL.LT.RNDOFF)) THEN IER = 911 NREC = 2 WRITE (REC,'(A,D13.5,A/A,D13.5,A)') &' ** You have set TOL = ',TOL,' which is not permitted. The', &' ** range of permitted values is (',RNDOFF,',0.01D0).' ELSE L = 1 20 CONTINUE IF (THRES(L).LT.TINY) THEN IER = 911 NREC = 2 WRITE (REC,'(A,I6,A,D13.5,A/A,D13.5,A)') &' ** You have set THRES(',L,') to be ',THRES(L),' which is ', &' ** less than the permitted minimum,',TINY,'.' END IF L = L + 1 IF (IER.NE.911 .AND. L.LE.NEQ) GO TO 20 END IF C C Return if error detected C IF (IER.NE.1) GO TO 80 C C Set formula definitions and characteristics by means of arguments C in the call list and COMMON blocks /RKCOM4/ and /RKCOM5/ C CALL CONST(METHOD,VECSTG,REQSTG,LINTPL) C C Set options in /RKCOM8/ UTASK = TASK1 .EQ. 'U' .OR. TASK1 .EQ. 'u' MSG = MESAGE C C Initialise problem status in /RKCOM1/ and /RKCOM2/ NEQN = NEQ TSTRT = TSTART TND = TEND T = TSTART TOLD = TSTART DIR = SIGN(ONE,TEND-TSTART) C C In CT the first step taken will have magnitude H. If HSTRT = ABS(HSTART) C is not equal to zero, H = HSTRT. If HSTRT is equal to zero, the code is C to find an on-scale initial step size H. To start this process, H is set C here to an upper bound on the first step size that reflects the scale of C the independent variable. UT has some additional information, namely the C first output point, that is used to refine this bound in UT when UTASK C is .TRUE.. If HSTRT is not zero, but it is either too big or too small, C the input HSTART is ignored and HSTRT is set to zero to activate the C automatic determination of an on-scale initial step size. C HSTRT = ABS(HSTART) HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTART),ABS(TEND))) IF (HSTRT.GT.ABS(TEND-TSTART) .OR. HSTRT.LT.HMIN) HSTRT = ZERO IF (HSTRT.EQ.ZERO) THEN H = MAX(ABS(TEND-TSTART)/RS3,HMIN) ELSE H = HSTRT END IF HOLD = ZERO TOLR = TOL NFCN = 0 SVNFCN = 0 OKSTP = 0 FLSTP = 0 FIRST = .TRUE. LAST = .FALSE. C C WORK(*) is partioned into a number of arrays using pointers. These C pointers are set in /RKCOM3/. PRTHRS = 1 C the threshold values PRERST = PRTHRS + NEQ C the error estimates PRWT = PRERST + NEQ C the weights used in the local error test PRYOLD = PRWT + NEQ C the previous value of the solution PRSCR = PRYOLD + NEQ C scratch array used for the higher order C approximate solution and for the previous C value of the derivative of the solution PRY = PRSCR + NEQ C the dependent variables PRYP = PRY + NEQ C the derivatives PRSTGS = PRYP + NEQ C intermediate stages held in an internal C array STAGES(NEQ,VECSTG) C FREEPR = PRSTGS + VECSTG*NEQ C C Allocate storage for interpolation if the TASK = `U' or `u' was C specified. INTP and LINTPL returned by CONST indicate whether there C is an interpolation scheme associated with the pair and how much C storage is required. C PRINTP = 1 LNINTP = 1 IF (UTASK) THEN IF (INTP) THEN LNINTP = LINTPL*NEQ IF (REQSTG) THEN PRINTP = FREEPR FREEPR = PRINTP + LNINTP ELSE PRINTP = PRSTGS FREEPR = MAX(PRINTP+VECSTG*NEQ,PRINTP+LNINTP) END IF END IF END IF C C Initialise state and allocate storage for global error assessment C using /RKCOM6/ GNFCN = 0 MAXERR = ZERO LOCMAX = TSTART ERASON = ERRASS ERASFL = .FALSE. IF (ERRASS) THEN C C Storage is required for the stages of a secondary integration. The C stages of the primary intergration can only be overwritten in the C cases where there is no interpolant or the interpolant does not C require information about the stages (e.g. METHOD 3 and METHOD 1, C respectively). IF (.NOT.REQSTG) THEN PRZSTG = PRSTGS ELSE PRZSTG = FREEPR FREEPR = PRZSTG + VECSTG*NEQ END IF PRZY = FREEPR PRZYP = PRZY + NEQ PRZERS = PRZYP + NEQ PRZERR = PRZERS + NEQ PRZYNU = PRZERR + NEQ FREEPR = PRZYNU + NEQ ELSE PRZSTG = 1 PRZY = 1 PRZYP = 1 PRZERS = 1 PRZERR = 1 PRZYNU = 1 END IF C LREQ = FREEPR - 1 C C Check for enough workspace and suitable range of integration C IF (LENWRK.LT.LREQ) THEN IER = 911 NREC = 2 WRITE (REC,'(A/A,I6,A,I6,A)') &' ** You have not supplied enough workspace. You gave LENWRK ', &' ** as', LENWRK, ', but it must be at least ',LREQ,'.' ELSE HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTART),ABS(TEND))) IF (ABS(TEND-TSTART).LT.HMIN) THEN IER = 911 NREC = 4 WRITE (REC,'(A/A/A,D13.5/A,D13.5,A)') &' ** You have set values for TEND and TSTART that are not ', &' ** clearly distinguishable for the method and the precision ', &' ** of the computer being used. ABS(TEND-TSTART) is ', &ABS(TEND-TSTART), &' ** but should be at least ',HMIN,'.' END IF END IF C C Return if error detected C IF (IER.NE.1) GO TO 80 C C Initialize elements of the workspace DO 40 L = 1, NEQ WORK(PRTHRS-1+L) = THRES(L) WORK(PRY-1+L) = YSTART(L) 40 CONTINUE C C Initialize the global error to zero when ERRASS = .TRUE. IF (ERRASS) THEN DO 60 L = 1, NEQ WORK(PRZERR-1+L) = ZERO 60 CONTINUE END IF C 80 CONTINUE C CALL RKMSG(IER,SRNAME,NREC,FLAG) C RETURN END SUBROUTINE SETUP C SUBROUTINE UT(F,TWANT,TGOT,YGOT,YPGOT,YMAX,WORK,UFLAG) C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C If you are not familiar with the code UT and how it is used in C conjunction with SETUP to solve initial value problems, you should study C the document file rksuite.doc carefully before proceeding further. The C following "Brief Reminder" is intended only to remind you of the meaning, C type, and size requirements of the arguments. C C NAME DECLARED IN AN EXTERNAL STATEMENT IN THE CALLING PROGRAM: C C F - name of the subroutine for evaluating the differential C equations. C C The subroutine F must have the form C C SUBROUTINE F(T,Y,YP) C DOUBLE PRECISION T,Y(*),YP(*) C Given input values of the independent variable T and the solution C components Y(*), for each L = 1,2,...,NEQ evaluate the differential C equation for the derivative of the Ith solution component and place the C value in YP(L). Do not alter the input values of T and Y(*). C RETURN C END C C INPUT VARIABLE C C TWANT - DOUBLE PRECISION C The next value of the independent variable where a C solution is desired. C C Constraints: TWANT must lie between the previous value C of TGOT (TSTART on the first call) and TEND. TWANT can be C equal to TEND, but it must be clearly distinguishable from C the previous value of TGOT (TSTART on the first call) in C the precision available. C C OUTPUT VARIABLES C C TGOT - DOUBLE PRECISION C A solution has been computed at this value of the C independent variable. C YGOT(*) - DOUBLE PRECISION array of length NEQ C Approximation to the true solution at TGOT. Do not alter C the contents of this array C YPGOT(*) - DOUBLE PRECISION array of length NEQ C Approximation to the first derivative of the true C solution at TGOT. C YMAX(*) - DOUBLE PRECISION array of length NEQ C YMAX(L) is the largest magnitude of YGOT(L) computed at any C time in the integration from TSTART to TGOT. Do not alter C the contents of this array. C C WORKSPACE C C WORK(*) - DOUBLE PRECISION array as used in SETUP C Do not alter the contents of this array. C C OUTPUT VARIABLE C C UFLAG - INTEGER C C SUCCESS. TGOT = TWANT. C = 1 - Complete success. C C "SOFT" FAILURES C = 2 - Warning: You are using METHOD = 3 inefficiently C by computing answers at many values of TWANT. If C you really need answers at so many specific points, C it would be more efficient to compute them with C METHOD = 2. To do this you would need to restart C from TGOT, YGOT(*) by a call to SETUP. If you wish C to continue as you are, you may. C = 3 - Warning: A considerable amount of work has been C expended. If you wish to continue on to TWANT, just C call UT again. C = 4 - Warning: It appears that this problem is "stiff". C You really should change to another code that is C intended for such problems, but if you insist, you can C continue with UT by calling it again. C C "HARD" FAILURES C = 5 - You are asking for too much accuracy. You cannot C continue integrating this problem. C = 6 - The global error assessment may not be reliable beyond C the current point in the integration. You cannot C continue integrating this problem. C C "CATASTROPHIC" FAILURES C = 911 - The nature of the catastrophe is reported on C the standard output channel. Unless special C provision was made in advance (see rksuite.doc), C the computation then comes to a STOP. C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C .. Scalar Arguments .. DOUBLE PRECISION TGOT, TWANT INTEGER UFLAG C .. Array Arguments .. DOUBLE PRECISION WORK(*), YGOT(*), YMAX(*), YPGOT(*) C .. Subroutine Arguments .. EXTERNAL F C .. Common Block for Problem Definition .. DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR INTEGER NEQN COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN SAVE /RKCOM1/ C .. Common Block to hold Problem Status .. DOUBLE PRECISION T, H, TOLD, HOLD INTEGER NFCN, SVNFCN, OKSTP, FLSTP LOGICAL FIRST, LAST COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP, & FIRST, LAST SAVE /RKCOM2/ C .. Common Block for General Workspace Pointers .. INTEGER PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP, & PRSTGS, PRINTP, LNINTP COMMON /RKCOM3/ PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP, & PRSTGS, PRINTP, LNINTP SAVE /RKCOM3/ C .. Common Block to hold Formula Definitions .. DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6), & E(7) INTEGER PTR(13), NSTAGE, METHD, MINTP LOGICAL INTP COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD, & MINTP, INTP SAVE /RKCOM4/ C .. Common Block to hold Formula Characterisitcs .. DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4 INTEGER ORDER, LSTSTG, MAXTRY, NSEC LOGICAL FSAL COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY, & NSEC, FSAL SAVE /RKCOM5/ C .. Common Block for Environment Parameters .. DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY INTEGER OUTCH COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY, & OUTCH SAVE /RKCOM7/ C .. Common Block for Integrator Options .. LOGICAL MSG, UTASK COMMON /RKCOM8/ MSG, UTASK SAVE /RKCOM8/ C .. Common Block for Error Message .. CHARACTER*80 REC(10) COMMON /RKCOM9/ REC SAVE /RKCOM9/ C .. Parameters .. CHARACTER*6 SRNAME PARAMETER (SRNAME='UT') LOGICAL ASK, TELL PARAMETER (ASK=.TRUE.,TELL=.FALSE.) INTEGER MINUS1, MINUS2 PARAMETER (MINUS1=-1,MINUS2=-2) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D+0) C .. Local Scalars .. DOUBLE PRECISION HMIN, TLAST, TNOW, UTEND INTEGER CFLAG, IER, L, NREC, STATE LOGICAL BADERR, GOBACK C .. External Subroutines .. EXTERNAL CHKFL, CT, INTRP, RESET, RKMSG, RKSIT C .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN C .. Save statement .. SAVE UTEND, TLAST C .. Executable Statements .. IER = 1 NREC = 0 GOBACK = .FALSE. BADERR = .FALSE. C C Is it permissible to call UT? C CALL RKSIT(ASK,'SETUP',STATE) IF (STATE.EQ.911) THEN IER = 912 NREC = 1 WRITE (REC,'(A)') &' ** A catastrophic error has already been detected elsewhere.' GO TO 100 END IF IF (STATE.EQ.MINUS1) THEN IER = 911 NREC = 1 WRITE (REC,'(A)') &' ** You have not called SETUP, so you cannot use UT.' GO TO 100 END IF IF (.NOT.UTASK) THEN IER = 911 NREC = 2 WRITE (REC,'(A/A)') &' ** You have called UT after you specified in SETUP that ', &' ** you were going to use CT. This is not permitted.' GO TO 100 END IF CALL RKSIT(ASK,SRNAME,STATE) IF (STATE.EQ.5 .OR. STATE.EQ.6) THEN IER = 911 NREC = 1 WRITE (REC,'(A/A)') &' ** This routine has already returned with a hard failure.', &' ** You must call SETUP to start another problem.' GO TO 100 END IF STATE = MINUS2 CALL RKSIT(TELL,SRNAME,STATE) C IF (FIRST) THEN C C First call. C C A value of TND is specified in SETUP. When INTP = .FALSE., as with C METHD = 3, output is obtained at the specified TWANT by resetting TND C to TWANT. At this point, before the integration gets started, this can C be done with a simple assignment. Later it is done with a call to RESET. C The original TND is SAVEd as a local variable UTEND. C UTEND = TND IF (.NOT.INTP) TND = TWANT C C The last TGOT returned is SAVEd in the variable TLAST. T (a variable C passed through the common block RKCOM2) records how far the integration C has advanced towards the specified TND. When output is obtained by C interpolation, the integration goes past the TGOT returned (T is closer C to the specified TND than TGOT). Initialize these variables and YMAX(*). TLAST = TSTRT TGOT = TSTRT DO 20 L = 1, NEQN YMAX(L) = ABS(WORK(PRY-1+L)) 20 CONTINUE C C If the code is to find an on-scale initial step size H, a bound was placed C on H in SETUP. Here the first output point is used to refine this bound. IF (HSTRT.EQ.ZERO) THEN H = MIN(ABS(H),ABS(TWANT-TSTRT)) HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTRT),ABS(TND))) H = MAX(H,HMIN) END IF C ELSE C C Subsequent call. C IF (TLAST.EQ.UTEND) THEN IER = 911 NREC = 3 WRITE (REC,'(A/A/A)') &' ** You have called UT after reaching TEND. (Your last ', &' ** call to UT resulted in TGOT = TEND.) To start a new ', &' ** problem, you will need to call SETUP.' GO TO 100 END IF C END IF C C Check for valid TWANT. C IF (DIR*(TWANT-TLAST).LE.ZERO) THEN IER = 911 NREC = 4 WRITE (REC,'(A/A/A/A)') &' ** You have made a call to UT with a TWANT that does ', &' ** not lie between the previous value of TGOT (TSTART ', &' ** on the first call) and TEND. This is not permitted. ', &' ** Check your program carefully.' GO TO 100 END IF IF (DIR*(TWANT-UTEND).GT.ZERO) THEN HMIN = MAX(TINY,TOOSML*MAX(ABS(TWANT),ABS(UTEND))) IF (ABS(TWANT-UTEND).LT.HMIN) THEN IER = 911 NREC = 5 WRITE (REC,'(A/A/A/A)') &' ** You have made a call to UT with a TWANT that does ', &' ** not lie between the previous value of TGOT (TSTART on ', &' ** the first call) and TEND. This is not permitted. TWANT ', &' ** is very close to TEND, so you may have meant to set ', &' ** it to be TEND exactly. Check your program carefully. ' ELSE IER = 911 NREC = 4 WRITE (REC,'(A/A/A/A)') &' ** You have made a call to UT with a TWANT that does ', &' ** not lie between the previous value of TGOT (TSTART ', &' ** on the first call) and TEND. This is not permitted. ', &' ** Check your program carefully.' END IF GO TO 100 END IF IF (.NOT.INTP) THEN HMIN = MAX(TINY,TOOSML*MAX(ABS(TLAST),ABS(TWANT))) IF (ABS(TWANT-TLAST).LT.HMIN) THEN IER = 911 NREC = 4 WRITE (REC,'(A/A/A/A,D13.5,A)') &' ** You have made a call to UT with a TWANT that is not ', &' ** sufficiently different from the last value of TGOT ', &' ** (TSTART on the first call). When using METHOD = 3, ', &' ** it must differ by at least ',HMIN,'.' GO TO 100 END IF C C We have a valid TWANT. There is no interpolation with this METHD and C therefore we step to TWANT exactly by resetting TND with a call to RESET. C On the first step this matter is handled differently as explained above. C IF (.NOT.FIRST) THEN CALL RESET(TWANT) CALL CHKFL(ASK,BADERR) IF (BADERR) GO TO 100 END IF END IF C C Process output, decide whether to take another step. C 40 CONTINUE C IF (INTP) THEN C C Interpolation is possible with this METHD. The integration has C already reached T. If this is past TWANT, GOBACK is set .TRUE. and C the answers are obtained by interpolation. C GOBACK = DIR*(T-TWANT) .GE. ZERO IF (GOBACK) THEN CALL INTRP(TWANT,'Both solution and derivative',NEQN,YGOT, & YPGOT,F,WORK,WORK(PRINTP),LNINTP) CALL CHKFL(ASK,BADERR) IF (BADERR) GO TO 100 TGOT = TWANT END IF ELSE C C Interpolation is not possible with this METHD, so output is obtained C by integrating to TWANT = TND. Both YGOT(*) and YPGOT(*) are then C already loaded with the solution at TWANT by CT. C GOBACK = T .EQ. TWANT IF (GOBACK) TGOT = TWANT END IF C C Updating of YMAX(*) is done here to account for the fact that when C interpolation is done, the integration goes past TGOT. Note that YGOT(*) C is not defined until CT is called. YMAX(*) was initialized at TSTRT C from values stored in WORK(*), so only needs to be updated for T C different from TSTRT. IF (T.NE.TSTRT) THEN DO 60 L = 1, NEQN YMAX(L) = MAX(YMAX(L),ABS(YGOT(L))) 60 CONTINUE END IF C C If done, go to the exit point. IF (GOBACK) GO TO 100 C C Take a step with CT in the direction of TND. On exit, the solution is C advanced to TNOW. The way CT is written, the approximate solution at C TNOW is available in both YGOT(*) and in WORK(*). If output is obtained by C stepping to the end (TNOW = TWANT = TND), YGOT(*) can be returned directly. C If output is obtained by interpolation, the subroutine INTRP that does this C uses the values in WORK(*) for its computations and places the approximate C solution at TWANT in the array YGOT(*) for return to the calling program. C The approximate derivative is handled in the same way. TNOW is output from C CT and is actually a copy of T declared above in a common block. C CALL CT(F,TNOW,YGOT,YPGOT,WORK,CFLAG) IER = CFLAG C C A successful step by CT is indicated by CFLAG = 1 or = 2. IF (CFLAG.EQ.1) THEN GO TO 40 ELSE IF (CFLAG.EQ.2) THEN C C Supplement the warning message written in CT. NREC = 3 WRITE (REC,'(A/A/A)') &' ** The last message was produced on a call to CT from UT. ', &' ** In UT the appropriate action is to change to METHOD = 2,', &' ** or, if insufficient memory is available, to METHOD = 1. ' ELSE IF (CFLAG.LE.6) THEN NREC = 1 WRITE (REC,'(A)') &' ** The last message was produced on a call to CT from UT.' ELSE BADERR = .TRUE. END IF TGOT = T C C Update YMAX(*) before the return. DO 80 L = 1, NEQN YMAX(L) = MAX(YMAX(L),ABS(YGOT(L))) 80 CONTINUE C C Exit point for UT. C 100 CONTINUE C IF (BADERR) THEN IER = 911 NREC = 4 WRITE (REC,'(A/A/A/A)') &' ** An internal call by UT to a subroutine resulted in an ', &' ** error that should not happen. Check your program ', &' ** carefully for array sizes, correct number of arguments,', &' ** type mismatches ... .' END IF C TLAST = TGOT C C All exits are done here after a call to RKMSG to report C what happened and set UFLAG. C CALL RKMSG(IER,SRNAME,NREC,UFLAG) C RETURN END SUBROUTINE UT C SUBROUTINE STAT(TOTFCN,STPCST,WASTE,STPSOK,HNEXT) C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C If you are not familiar with the code STAT and how it is used in C conjunction with the integrators CT and UT, you should study the C document file rksuite.doc carefully before attempting to use the code. C The following "Brief Reminder" is intended only to remind you of the C meaning, type, and size requirements of the arguments. C C STAT is called to obtain some details about the integration. C C OUTPUT VARIABLES C C TOTFCN - INTEGER C Total number of calls to F in the integration so far -- C a measure of the cost of the integration. C STPCST - INTEGER C Cost of a typical step with this METHOD measured in C calls to F. C WASTE - DOUBLE PRECISION C The number of attempted steps that failed to meet the C local error requirement divided by the total number of C steps attempted so far in the integration. C STPSOK - INTEGER C The number of accepted steps. C HNEXT - DOUBLE PRECISION C The step size the integrator plans to use for the next step. C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C .. Scalar Arguments .. DOUBLE PRECISION HNEXT, WASTE INTEGER STPCST, STPSOK, TOTFCN C .. Common Block to hold Problem Status .. DOUBLE PRECISION T, H, TOLD, HOLD INTEGER NFCN, SVNFCN, OKSTP, FLSTP LOGICAL FIRST, LAST COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP, & FIRST, LAST SAVE /RKCOM2/ C .. Common Block to hold Formula Characterisitcs .. DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4 INTEGER ORDER, LSTSTG, MAXTRY, NSEC LOGICAL FSAL COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY, & NSEC, FSAL SAVE /RKCOM5/ C .. Common Block for Global Error Assessment .. DOUBLE PRECISION MAXERR, LOCMAX INTEGER GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR, & PRZYNU LOGICAL ERASON, ERASFL COMMON /RKCOM6/ MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP, & PRZERS, PRZERR, PRZYNU, ERASON, ERASFL SAVE /RKCOM6/ C .. Common Block for Integrator Options .. LOGICAL MSG, UTASK COMMON /RKCOM8/ MSG, UTASK SAVE /RKCOM8/ C .. Common Block for Error Message .. CHARACTER*80 REC(10) COMMON /RKCOM9/ REC SAVE /RKCOM9/ C .. Parameters .. CHARACTER*6 SRNAME PARAMETER (SRNAME='STAT') LOGICAL ASK INTEGER MINUS1, MINUS2 PARAMETER (ASK=.TRUE.,MINUS1=-1,MINUS2=-2) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D+0) C .. Local Scalars .. INTEGER FLAG, IER, NREC, STATE C .. External Subroutines .. EXTERNAL RKMSG, RKSIT C .. Intrinsic Functions .. INTRINSIC DBLE C .. Executable Statements .. C IER = 1 NREC = 0 C C Is it permissible to call STAT? C CALL RKSIT(ASK,SRNAME,STATE) IF (STATE.EQ.911) THEN IER = 912 NREC = 1 WRITE (REC,'(A)') &' ** A catastrophic error has already been detected elsewhere.' GO TO 20 END IF IF (STATE.EQ.MINUS2) THEN IER = 911 NREC = 3 WRITE (REC,'(A/A/A)') &' ** You have already made a call to STAT after a hard ', &' ** failure was reported from the integrator. You cannot', &' ** call STAT again.' GO TO 20 END IF CALL RKSIT(ASK,'CT',STATE) IF (STATE.EQ.MINUS1) THEN IER = 911 NREC = 1 IF (UTASK) THEN WRITE (REC,'(A)') &' ** You have not called UT, so you cannot use STAT.' ELSE WRITE (REC,'(A)') &' ** You have not called CT, so you cannot use STAT.' END IF GO TO 20 END IF C C Set flag so that the routine can only be called once after a hard C failure from the integrator. IF (STATE.EQ.5 .OR. STATE.EQ.6) IER = MINUS2 C TOTFCN = SVNFCN + NFCN IF (ERASON) TOTFCN = TOTFCN + GNFCN STPCST = COST STPSOK = OKSTP IF (OKSTP.LE.1) THEN WASTE = ZERO ELSE WASTE = DBLE(FLSTP)/DBLE(FLSTP+OKSTP) END IF HNEXT = H C 20 CONTINUE C CALL RKMSG(IER,SRNAME,NREC,FLAG) C RETURN END SUBROUTINE STAT C SUBROUTINE GLBERR(RMSERR,ERRMAX,TERRMX,WORK) C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C If you are not familiar with the code GLBERR and how it is used in C conjunction with UT and CT to solve initial value problems, you should C study the document file rksuite.doc carefully before attempting to use C the code. The following "Brief Reminder" is intended only to remind you C of the meaning, type, and size requirements of the arguments. C C If ERRASS was set .TRUE. in the call to SETUP, then after any call to UT C or CT to advance the integration to TNOW or TWANT, the subroutine GLBERR C may be called to obtain an assessment of the true error of the integration. C At each step and for each solution component Y(L), a more accurate "true" C solution YT(L), an average magnitude "size(L)" of its size, and its error C abs(Y(L) - YT(L))/max("size(L)",THRES(L)) C are computed. The assessment returned in RMSERR(L) is the RMS (root-mean- C square) average of the error in the Lth solution component over all steps C of the integration from TSTART through TNOW. C C OUTPUT VARIABLES C C RMSERR(*) - DOUBLE PRECISION array of length NEQ C RMSERR(L) approximates the RMS average of the true error C of the numerical solution for the Ith solution component, C L = 1,2,...,NEQ. The average is taken over all steps from C TSTART to TNOW. C ERRMAX - DOUBLE PRECISION C The maximum (approximate) true error taken over all C solution components and all steps from TSTART to TNOW. C TERRMX - DOUBLE PRECISION C First value of the independent variable where the C (approximate) true error attains the maximum value ERRMAX. C C WORKSPACE C C WORK(*) - DOUBLE PRECISION array as used in SETUP and UT or CT C Do not alter the contents of this array. C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C .. Scalar Arguments .. DOUBLE PRECISION ERRMAX, TERRMX C .. Array Arguments .. DOUBLE PRECISION RMSERR(*), WORK(*) C .. Common Block for Problem Definition .. DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR INTEGER NEQN COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN SAVE /RKCOM1/ C .. Common Block to hold Problem Status .. DOUBLE PRECISION T, H, TOLD, HOLD INTEGER NFCN, SVNFCN, OKSTP, FLSTP LOGICAL FIRST, LAST COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP, & FIRST, LAST SAVE /RKCOM2/ C .. Common Block for Global Error Assessment .. DOUBLE PRECISION MAXERR, LOCMAX INTEGER GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR, & PRZYNU LOGICAL ERASON, ERASFL COMMON /RKCOM6/ MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP, & PRZERS, PRZERR, PRZYNU, ERASON, ERASFL SAVE /RKCOM6/ C .. Common Block for Integrator Options .. LOGICAL MSG, UTASK COMMON /RKCOM8/ MSG, UTASK SAVE /RKCOM8/ C .. Common Block for Error Message .. CHARACTER*80 REC(10) COMMON /RKCOM9/ REC SAVE /RKCOM9/ C .. Parameters .. CHARACTER*6 SRNAME PARAMETER (SRNAME='GLBERR') LOGICAL ASK PARAMETER (ASK=.TRUE.) INTEGER MINUS1, MINUS2 PARAMETER (MINUS1=-1,MINUS2=-2) C .. Local Scalars .. INTEGER FLAG, IER, L, NREC, STATE C .. External Subroutines .. EXTERNAL RKMSG, RKSIT C .. Intrinsic Functions .. INTRINSIC SQRT C .. Executable Statements .. C IER = 1 NREC = 0 C C Is it permissible to call GLBERR? C CALL RKSIT(ASK,SRNAME,STATE) IF (STATE.EQ.911) THEN IER = 912 NREC = 1 WRITE (REC,'(A)') &' ** A catastrophic error has already been detected elsewhere.' GO TO 40 END IF IF (STATE.EQ.MINUS2) THEN IER = 911 NREC = 3 WRITE (REC,'(A/A/A)') &' ** You have already made a call to GLBERR after a hard ', &' ** failure was reported from the integrator. You cannot', &' ** call GLBERR again.' GO TO 40 END IF CALL RKSIT(ASK,'CT',STATE) IF (STATE.EQ.MINUS1) THEN IER = 911 NREC = 1 IF (UTASK) THEN WRITE (REC,'(A)') &' ** You have not yet called UT, so you cannot call GLBERR.' ELSE WRITE (REC,'(A)') &' ** You have not yet called CT, so you cannot call GLBERR.' END IF GO TO 40 END IF C C Set flag so that the routine can only be called once after a hard C failure from the integrator. IF (STATE.EQ.5 .OR. STATE.EQ.6) IER = MINUS2 C C Check that ERRASS was set properly for error assessment in SETUP. C IF (.NOT.ERASON) THEN IER = 911 NREC = 3 WRITE (REC,'(A/A/A)') &' ** No error assessment is available since you did not ', &' ** ask for it in your call to the routine SETUP.', &' ** Check your program carefully.' GO TO 40 END IF C C Check to see if the integrator has not actually taken a step. C IF (OKSTP.EQ.0) THEN IER = 911 NREC = 3 WRITE (REC,'(A/A/A)') &' ** The integrator has not actually taken any successful ', &' ** steps. You cannot call GLBERR in this circumstance. ', &' ** Check your program carefully.' GO TO 40 END IF C C Compute RMS error and set output variables. C ERRMAX = MAXERR TERRMX = LOCMAX DO 20 L = 1, NEQN RMSERR(L) = SQRT(WORK(PRZERR-1+L)/OKSTP) 20 CONTINUE C 40 CONTINUE C CALL RKMSG(IER,SRNAME,NREC,FLAG) C RETURN END SUBROUTINE GLBERR C SUBROUTINE CT(F,TNOW,YNOW,YPNOW,WORK,CFLAG) C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C If you are not familiar with the code CT and how it is used in C conjunction with SETUP to solve initial value problems, you should study C the document file rksuite.doc carefully before attempting to use the code. C The following "Brief Reminder" is intended only to remind you of the C meaning, type, and size requirements of the arguments. C C NAME DECLARED IN AN EXTERNAL STATEMENT IN THE CALLING PROGRAM: C C F - name of the subroutine for evaluating the differential C equations. C C The subroutine F must have the form C C SUBROUTINE F(T,Y,YP) C DOUBLE PRECISION T,Y(*),YP(*) C Using the input values of the independent variable T and the solution C components Y(*), for each L = 1,2,...,NEQ evaluate the differential C equation for the derivative of the Lth solution component and place the C value in YP(L). Do not alter the input values of T and Y(*). C RETURN C END C C OUTPUT VARIABLES C C TNOW - DOUBLE PRECISION C Current value of the independent variable. C YNOW(*) - DOUBLE PRECISION array of length NEQ C Approximation to the true solution at TNOW. C YPNOW(*) - DOUBLE PRECISION array of length NEQ C Approximation to the first derivative of the C true solution at TNOW. C C WORKSPACE C C WORK(*) - DOUBLE PRECISION array as used in SETUP C Do not alter the contents of this array. C C OUTPUT VARIABLE C C CFLAG - INTEGER C C SUCCESS. A STEP WAS TAKEN TO TNOW. C = 1 - Complete success. C C "SOFT" FAILURES C = 2 - Warning: You have obtained an answer by integrating C to TEND (TNOW = TEND). You have done this at least C 100 times, and monitoring of the computation reveals C that this way of getting output has degraded the C efficiency of the code. If you really need answers at C so many specific points, it would be more efficient to C get them with INTRP. (If METHOD = 3, you would need C to change METHOD and restart from TNOW, YNOW(*) by a C call to SETUP.) If you wish to continue as you are, C you may. C = 3 - Warning: A considerable amount of work has been C expended. To continue the integration, just call C CT again. C = 4 - Warning: It appears that this problem is "stiff". C You really should change to another code that is C intended for such problems, but if you insist, you C can continue with CT by calling it again. C C "HARD" FAILURES C = 5 - You are asking for too much accuracy. You cannot C continue integrating this problem. C = 6 - The global error assessment may not be reliable beyond C the current point in the integration. You cannot C continue integrating this problem. C C "CATASTROPHIC" FAILURES C = 911 - The nature of the catastrophe is reported on C the standard output channel. Unless special C provision was made in advance (see rksuite.doc), C the computation then comes to a STOP. C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C .. Scalar Arguments .. DOUBLE PRECISION TNOW INTEGER CFLAG C .. Array Arguments .. DOUBLE PRECISION WORK(*), YNOW(*), YPNOW(*) C .. Subroutine Arguments .. EXTERNAL F C .. Common Block for Problem Definition .. DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR INTEGER NEQN COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN SAVE /RKCOM1/ C .. Common Block to hold Problem Status .. DOUBLE PRECISION T, H, TOLD, HOLD INTEGER NFCN, SVNFCN, OKSTP, FLSTP LOGICAL FIRST, LAST COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP, & FIRST, LAST SAVE /RKCOM2/ C .. Common Block for General Workspace Pointers .. INTEGER PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP, & PRSTGS, PRINTP, LNINTP COMMON /RKCOM3/ PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP, & PRSTGS, PRINTP, LNINTP SAVE /RKCOM3/ C .. Common Block to hold Formula Characterisitcs .. DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4 INTEGER ORDER, LSTSTG, MAXTRY, NSEC LOGICAL FSAL COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY, & NSEC, FSAL SAVE /RKCOM5/ C .. Common Block for Global Error Assessment .. DOUBLE PRECISION MAXERR, LOCMAX INTEGER GNFCN, PRZSTG, PRZY, PRZYP, PRZERS, PRZERR, & PRZYNU LOGICAL ERASON, ERASFL COMMON /RKCOM6/ MAXERR, LOCMAX, GNFCN, PRZSTG, PRZY, PRZYP, & PRZERS, PRZERR, PRZYNU, ERASON, ERASFL SAVE /RKCOM6/ C .. Common Block for Environment Parameters .. DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY INTEGER OUTCH COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY, & OUTCH SAVE /RKCOM7/ C .. Common Block for Integrator Options .. LOGICAL MSG, UTASK COMMON /RKCOM8/ MSG, UTASK SAVE /RKCOM8/ C .. Common Block for Error Message .. CHARACTER*80 REC(10) COMMON /RKCOM9/ REC SAVE /RKCOM9/ C .. Parameters .. CHARACTER*6 SRNAME PARAMETER (SRNAME='CT') LOGICAL ASK, TELL PARAMETER (ASK=.TRUE.,TELL=.FALSE.) INTEGER MINUS1, MINUS2 PARAMETER (MINUS1=-1,MINUS2=-2) INTEGER MAXFCN PARAMETER (MAXFCN=5000) DOUBLE PRECISION ZERO, PT1, PT9, ONE, TWO, HUNDRD PARAMETER (ZERO=0.0D+0,PT1=0.1D+0,PT9=0.9D+0,ONE=1.0D+0, & TWO=2.0D+0,HUNDRD=100.0D+0) C .. Local Scalars .. DOUBLE PRECISION ALPHA, BETA, ERR, ERROLD, HAVG, HMIN, HTRY, TAU, & TEMP1, TEMP2, YPNORM INTEGER IER, JFLSTP, L, NREC, NTEND, POINT, STATE, YNEW, & YPOLD LOGICAL CHKEFF, FAILED, MAIN, PHASE1, PHASE2, PHASE3, & TOOMCH C .. External Subroutines .. EXTERNAL RKMSG, RKSIT, STEP, STIFF, TRUERR C .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SIGN C .. Save statement .. SAVE JFLSTP, NTEND, ERROLD, HAVG, PHASE2, YNEW, & YPOLD, CHKEFF C .. Executable Statements .. C IER = 1 NREC = 0 C C Is it permissible to call CT? C CALL RKSIT(ASK,'SETUP',STATE) IF (STATE.EQ.911) THEN IER = 912 NREC = 1 WRITE (REC,'(A)') &' ** A catastrophic error has already been detected elsewhere.' GO TO 180 END IF IF (STATE.EQ.MINUS1) THEN IER = 911 NREC = 1 WRITE (REC,'(A)') &' ** You have not called SETUP, so you cannot use CT.' GO TO 180 END IF IF (UTASK) THEN CALL RKSIT(ASK,'UT',STATE) IF (STATE.NE.MINUS2) THEN IER = 911 NREC = 2 WRITE (REC,'(A/A)') &' ** You have called CT after you specified in SETUP that ', &' ** you were going to use UT. This is not permitted.' UTASK = .FALSE. GO TO 180 END IF END IF CALL RKSIT(ASK,SRNAME,STATE) IF (STATE.EQ.5 .OR. STATE.EQ.6) THEN IER = 911 NREC = 3 WRITE (REC,'(A/A)') &' ** CT has already returned with a flag value of 5 or 6.', &' ** You cannot continue integrating this problem. You must ', &' ** call SETUP to start another problem.' GO TO 180 END IF C IF (FIRST) THEN C C First call in an integration -- initialize everything. C CHKEFF = .FALSE. NTEND = 0 JFLSTP = 0 C C A scratch area of WORK(*) starting at PRSCR is used to hold two C arrays in this subroutine: the higher order approximate solution at C the end of a step and the approximate derivative of the solution at C the end of the last step. To make this clear, local pointers YNEW and C YPOLD are used. YNEW = PRSCR YPOLD = PRSCR C C For this first step T was initialized to TSTRT in SETUP and the C starting values YSTART(*) were loaded into the area of WORK(*) reserved C for the current solution approximation starting at location PRY. The C derivative is now computed and stored in WORK(*) starting at PRYP. C Subsequently these arrays are copied to the output vectors YNOW(*) C and YPNOW(*). CALL F(T,WORK(PRY),WORK(PRYP)) NFCN = NFCN + 1 DO 20 L = 1, NEQN YNOW(L) = WORK(PRY-1+L) YPNOW(L) = WORK(PRYP-1+L) 20 CONTINUE C C Set dependent variables for error assessment. IF (ERASON) THEN DO 40 L = 1, NEQN WORK(PRZY-1+L) = YNOW(L) WORK(PRZYP-1+L) = YPNOW(L) 40 CONTINUE END IF C C The weights for the control of the error depend on the size of the C solution at the beginning and at the end of the step. On the first C step we do not have all this information. Whilst determining the C initial step size we initialize the weight vector to the larger of C abs(Y(L)) and the threshold for this component. DO 60 L = 1, NEQN WORK(PRWT-1+L) = MAX(ABS(YNOW(L)),WORK(PRTHRS-1+L)) 60 CONTINUE C C If HSTRT is equal to zero, the code is to find an on-scale initial step C size H. CT has an elaborate scheme of three phases for finding such an H, C and some preparations are made earlier. In SETUP an upper bound is placed C on H that reflects the scale of the independent variable. When UTASK is C .TRUE., UT refines this bound using the first output point. Here in CT C PHASE1 applies a rule of thumb based on the error control, the order of the C the formula, and the size of the initial slope to get a crude approximation C to an on-scale H. PHASE2 may reduce H in the course of taking the first C step. PHASE3 repeatedly adjusts H and retakes the first step until H is C on scale. C C A guess for the magnitude of the first step size H can be provided to SETUP C as HSTART. If it is too big or too small, it is ignored and the automatic C determination of an on-scale initial step size is activated. If it is C acceptable, H is set to HSTART in SETUP. Even when H is supplied to CT, C PHASE3 of the scheme for finding an on-scale initial step size is made C active so that the code can deal with a bad guess. C PHASE1 = HSTRT .EQ. ZERO PHASE2 = PHASE1 PHASE3 = .TRUE. IF (PHASE1) THEN H = ABS(H) YPNORM = ZERO DO 80 L = 1, NEQN IF (ABS(YNOW(L)).NE.ZERO) THEN YPNORM = MAX(YPNORM,ABS(YPNOW(L))/WORK(PRWT-1+L)) END IF 80 CONTINUE TAU = TOLR**EXPON IF (H*YPNORM.GT.TAU) H = TAU/YPNORM HMIN = MAX(TINY,TOOSML*MAX(ABS(TSTRT),ABS(TND))) H = DIR*MAX(H,HMIN) PHASE1 = .FALSE. END IF C ELSE C C Continuation call C IF (LAST) THEN IER = 911 NREC = 3 WRITE (REC,'(A,D13.5,A/A/A)') &' ** You have already reached TEND ( = ', TND, ').', &' ** To integrate further with the same problem you must ', &' ** call the routine RESET with a new value of TEND.' GO TO 180 END IF END IF C C Begin computation of a step here. C FAILED = .FALSE. C 100 CONTINUE H = SIGN(ABS(H),DIR) C C Reduce the step size if necessary so that the code will not step C past TND. "Look ahead" to prevent unnecessarily small step sizes. LAST = DIR*((T+H)-TND) .GE. ZERO IF (LAST) THEN H = TND - T ELSE IF (DIR*((T+TWO*H)-TND).GE.ZERO) THEN H = (TND-T)/TWO END IF C C When the integrator is at T and attempts a step of H, the function C defining the differential equations will be evaluated at a number of C arguments between T and T+H. If H is too small, these arguments cannot C be clearly distinguished in the precision available. C HMIN = MAX(TINY,TOOSML*MAX(ABS(T),ABS(T+H))) IF (ABS(H).LT.HMIN) THEN IER = 5 NREC = 3 WRITE (REC,'(A/A,D13.5,A,D13.5/A)') &' ** In order to satisfy your error requirements CT would ', &' ** have to use a step size of ',H,' at TNOW = ',T, &' ** This is too small for the machine precision.' GO TO 180 END IF C C Monitor the impact of output on the efficiency of the integration. C IF (CHKEFF) THEN NTEND = NTEND + 1 IF (NTEND.GE.100 .AND. NTEND.GE.OKSTP/3) THEN IER = 2 NREC = 6 WRITE (REC,'(A/A/A/A/A/A)') &' ** More than 100 output points have been obtained by ', &' ** integrating to TEND. They have been sufficiently close ', &' ** to one another that the efficiency of the integration has ', &' ** been degraded. It would probably be (much) more efficient ', &' ** to obtain output by interpolating with INTRP (after ', &' ** changing to METHOD=2 if you are using METHOD = 3).' NTEND = 0 GO TO 180 END IF END IF C C Check for stiffness and for too much work. Stiffness can be C checked only after a successful step. C IF (.NOT.FAILED) THEN C C Check for too much work. TOOMCH = NFCN .GT. MAXFCN IF (TOOMCH) THEN IER = 3 NREC = 3 WRITE (REC,'(A,I6,A/A/A)') &' ** Approximately ',MAXFCN,' function evaluations have been ', &' ** used to compute the solution since the integration ', &' ** started or since this message was last printed.' C C After this warning message, NFCN is reset to permit the integration C to continue. The total number of function evaluations in the primary C integration is SVNFCN + NFCN. SVNFCN = SVNFCN + NFCN NFCN = 0 END IF C C Check for stiffness. NREC is passed on to STIFF because when C TOOMCH = .TRUE. and stiffness is diagnosed, the message about too C much work is augmented inside STIFF to explain that it is due to C stiffness. CALL STIFF(F,HAVG,JFLSTP,TOOMCH,MAXFCN,WORK,IER,NREC) C IF (IER.NE.1) GO TO 180 END IF C C Take a step. Whilst finding an on-scale H (PHASE2 = .TRUE.), the input C value of H might be reduced (repeatedly), but it will not be reduced C below HMIN. The local error is estimated, a weight vector is formed, C and a weighted maximum norm, ERR, of the local error is returned. C The variable MAIN is input as .TRUE. to tell STEP that this is the C primary, or "main", integration. C C H resides in the common block /RKCOM2/ which is used by both CT and STEP; C since it may be changed inside STEP, a local copy is made to ensure C portability of the code. C MAIN = .TRUE. HTRY = H CALL STEP(F,NEQN,T,WORK(PRY),WORK(PRYP),WORK(PRSTGS),TOLR,HTRY, & WORK(PRWT),WORK(YNEW),WORK(PRERST),ERR,MAIN,HMIN, & WORK(PRTHRS),PHASE2) H = HTRY C C Compare the norm of the local error to the tolerance. C IF (ERR.GT.TOLR) THEN C C Failed step. Reduce the step size and try again. C C First step: Terminate PHASE3 of the search for an on-scale step size. C The step size is not on scale, so ERR may not be accurate; C reduce H by a fixed factor. Failed attempts to take the C first step are not counted. C Later step: Use ERR to compute an "optimal" reduction of H. More than C one failure indicates a difficulty with the problem and an C ERR that may not be accurate, so reduce H by a fixed factor. C IF (FIRST) THEN PHASE3 = .FALSE. ALPHA = RS1 ELSE FLSTP = FLSTP + 1 JFLSTP = JFLSTP + 1 IF (FAILED) THEN ALPHA = RS1 ELSE ALPHA = SAFETY*(TOLR/ERR)**EXPON ALPHA = MAX(ALPHA,RS1) END IF END IF H = ALPHA*H FAILED = .TRUE. GO TO 100 END IF C C Successful step. C C Predict a step size appropriate for the next step. After the first C step the prediction can be refined using an idea of H.A. Watts that C takes account of how well the prediction worked on the previous step. BETA = (ERR/TOLR)**EXPON IF (.NOT.FIRST) THEN TEMP1 = (ERR**EXPON)/H TEMP2 = (ERROLD**EXPON)/HOLD IF (TEMP1.LT.TEMP2*HUNDRD .AND. TEMP2.LT.TEMP1*HUNDRD) THEN BETA = BETA*(TEMP1/TEMP2) END IF END IF ALPHA = RS3 IF (SAFETY.LT.BETA*ALPHA) ALPHA = SAFETY/BETA C C On the first step a search is made for an on-scale step size. PHASE2 C of the scheme comes to an end here because a step size has been found C that is both successful and has a credible local error estimate. Except C in the special case that the first step is also the last, the step is C repeated in PHASE3 as long as an increase greater than RS2 appears C possible. An increase as big as RS3 is permitted. A step failure C terminates PHASE3. C IF (FIRST) THEN PHASE2 = .FALSE. PHASE3 = PHASE3 .AND. .NOT. LAST .AND. (ALPHA.GT.RS2) IF (PHASE3) THEN H = ALPHA*H GO TO 100 END IF END IF C C After getting on scale, step size changes are more restricted. ALPHA = MIN(ALPHA,RS) IF (FAILED) ALPHA = MIN(ALPHA,ONE) ALPHA = MAX(ALPHA,RS1) HOLD = H H = ALPHA*H C C For the diagnosis of stiffness, an average accepted step size, HAVG, C must be computed and SAVEd. IF (FIRST) THEN HAVG = HOLD ELSE HAVG = PT9*HAVG + PT1*HOLD END IF C FIRST = .FALSE. ERROLD = ERR TOLD = T C C Take care that T is set to precisely TND when the end of the C integration is reached. IF (LAST) THEN T = TND ELSE T = T + HOLD END IF C C Increment counter on accepted steps. Note that successful steps C that are repeated whilst getting on scale are not counted. OKSTP = OKSTP + 1 C C Advance the current solution and its derivative. (Stored in WORK(*) C with the first location being PRY and PRYP, respectively.) Update the C previous solution and its derivative. (Stored in WORK(*) with the first C location being PRYOLD and YPOLD, respectively.) Note that the previous C derivative will overwrite YNEW(*). C DO 120 L = 1, NEQN WORK(PRYOLD-1+L) = WORK(PRY-1+L) WORK(PRY-1+L) = WORK(YNEW-1+L) WORK(YPOLD-1+L) = WORK(PRYP-1+L) 120 CONTINUE C IF (FSAL) THEN C C When FSAL = .TRUE., YP(*) is the last stage of the step. POINT = PRSTGS + (LSTSTG-1)*NEQN DO 140 L = 1, NEQN WORK(PRYP-1+L) = WORK(POINT-1+L) 140 CONTINUE ELSE C C Call F to evaluate YP(*). CALL F(T,WORK(PRY),WORK(PRYP)) NFCN = NFCN + 1 END IF C C If global error assessment is desired, advance the secondary C integration from TOLD to T. C IF (ERASON) THEN CALL TRUERR(F,NEQN,WORK(PRY),TOLR,WORK(PRWT),WORK(PRZY), & WORK(PRZYP),WORK(PRZERR),WORK(PRZYNU),WORK(PRZERS), & WORK(PRZSTG),IER) IF (IER.EQ.6) THEN C C The global error estimating procedure has broken down. Treat it as a C failed step. The solution and derivative are reset to their values at C the beginning of the step since the last valid error assessment refers C to them. OKSTP = OKSTP - 1 ERASFL = .TRUE. LAST = .FALSE. T = TOLD H = HOLD DO 160 L = 1, NEQN WORK(PRY-1+L) = WORK(PRYOLD-1+L) WORK(PRYP-1+L) = WORK(YPOLD-1+L) 160 CONTINUE IF (OKSTP.GT.1) THEN NREC = 2 WRITE (REC,'(A/A,D13.5,A)') &' ** The global error assessment may not be reliable for T past ', &' ** TNOW = ',T,'. The integration is being terminated.' ELSE NREC = 2 WRITE (REC,'(A/A)') &' ** The global error assessment algorithm failed at the start', &' ** the integration. The integration is being terminated.' END IF GO TO 180 END IF END IF C C C Exit point for CT C 180 CONTINUE C C Set the output variables and flag that interpolation is permitted C IF (IER.LT.911) THEN TNOW = T LAST = TNOW .EQ. TND CHKEFF = LAST DO 200 L = 1, NEQN YNOW(L) = WORK(PRY-1+L) YPNOW(L) = WORK(PRYP-1+L) 200 CONTINUE IF (IER.EQ.1) THEN STATE = MINUS2 CALL RKSIT(TELL,'INTRP',STATE) END IF END IF C C Call RKMSG to report what happened and set CFLAG. C CALL RKMSG(IER,SRNAME,NREC,CFLAG) C RETURN END SUBROUTINE CT C SUBROUTINE INTRP(TWANT,REQEST,NWANT,YWANT,YPWANT,F,WORK,WRKINT, & LENINT) C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C If you are not familiar with the code INTRP and how it is used in C conjunction with CT to solve initial value problems, you should study the C document file rksuite.doc carefully before attempting to use the code. The C following "Brief Reminder" is intended only to remind you of the meaning, C type, and size requirements of the arguments. C C When integrating with METHOD = 1 or 2, answers may be obtained inexpensively C between steps by interpolation. INTRP is called after a step by CT from a C previous value of T, called TOLD below, to the current value of T to get C an answer at TWANT. You can specify any value of TWANT you wish, but C specifying a value outside the interval [TOLD,T] is likely to yield C answers with unsatisfactory accuracy. C C INPUT VARIABLE C C TWANT - DOUBLE PRECISION C The value of the independent variable where a solution C is desired. C C The interpolant is to be evaluated at TWANT to approximate the solution C and/or its first derivative there. There are three cases: C C INPUT VARIABLE C C REQEST - CHARACTER*(*) C Only the first character of REQEST is significant. C REQEST(1:1) = `S' or `s'- compute approximate `S'olution C only. C = `D' or `d'- compute approximate first C `D'erivative of the solution only. C = `B' or `b'- compute `B'oth approximate solution C and first derivative. C Constraint: REQEST(1:1) must be `S',`s',`D',`d',`B' or `b'. C C If you intend to interpolate at many points, you should arrange for the C the interesting components to be the first ones because the code C approximates only the first NWANT components. C C INPUT VARIABLE C C NWANT - INTEGER C Only the first NWANT components of the answer are to be C computed. C Constraint: NEQ >= NWANT >= 1 C C OUTPUT VARIABLES C C YWANT(*) - DOUBLE PRECISION array of length NWANT C Approximation to the first NWANT components of the true C solution at TWANT when REQESTed. C YPWANT(*) - DOUBLE PRECISION array of length NWANT C Approximation to the first NWANT components of the first C derivative of the true solution at TWANT when REQESTed. C C NAME DECLARED IN AN EXTERNAL STATEMENT IN THE PROGRAM CALLING INTRP: C C F - name of the subroutine for evaluating the differential C equations as provided to CT. C C WORKSPACE C C WORK(*) - DOUBLE PRECISION array as used in SETUP and CT C Do not alter the contents of this array. C C WRKINT(*) - DOUBLE PRECISION array of length LENINT C Do not alter the contents of this array. C C LENINT - INTEGER C Length of WRKINT. If C METHOD = 1, LENINT must be at least 1 C = 2, LENINT must be at least NEQ+MAX(NEQ,5*NWANT) C = 3--CANNOT BE USED WITH THIS SUBROUTINE C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C .. Scalar Arguments .. DOUBLE PRECISION TWANT INTEGER LENINT, NWANT CHARACTER*(*) REQEST C .. Array Arguments .. DOUBLE PRECISION WORK(*), WRKINT(*), YPWANT(*), YWANT(*) C .. Subroutine Arguments .. EXTERNAL F C .. Common Block for Problem Definition .. DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR INTEGER NEQN COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN SAVE /RKCOM1/ C .. Common Block for General Workspace Pointers .. INTEGER PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP, & PRSTGS, PRINTP, LNINTP COMMON /RKCOM3/ PRTHRS, PRERST, PRWT, PRYOLD, PRSCR, PRY, PRYP, & PRSTGS, PRINTP, LNINTP SAVE /RKCOM3/ C .. Common Block to hold Formula Definitions .. DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6), & E(7) INTEGER PTR(13), NSTAGE, METHD, MINTP LOGICAL INTP COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD, & MINTP, INTP SAVE /RKCOM4/ C .. Common Block for Integrator Options .. LOGICAL MSG, UTASK COMMON /RKCOM8/ MSG, UTASK SAVE /RKCOM8/ C .. Common Block for Error Message .. CHARACTER*80 REC(10) COMMON /RKCOM9/ REC SAVE /RKCOM9/ C .. Parameters .. CHARACTER*6 SRNAME PARAMETER (SRNAME='INTRP') LOGICAL ASK INTEGER PLUS1, MINUS1, MINUS2 PARAMETER (ASK=.TRUE.,PLUS1=1,MINUS1=-1,MINUS2=-2) C .. Local Scalars .. INTEGER FLAG, ICHK, IER, NREC, NWNTSV, STARTP, STATE, & STATE1 LOGICAL ININTP, LEGALR CHARACTER REQST1 C .. External Subroutines .. EXTERNAL EVALI, FORMI, RKMSG, RKSIT C .. Intrinsic Functions .. INTRINSIC MAX C .. Save statement .. SAVE NWNTSV, ININTP, STARTP C .. Data statements .. DATA NWNTSV/MINUS1/ C .. Executable Statements .. C IER = 1 NREC = 0 C C Is it permissible to call INTRP? C CALL RKSIT(ASK,'CT',STATE) IF (STATE.EQ.911) THEN IER = 912 NREC = 1 WRITE (REC,'(A)') &' ** A catastrophic error has already been detected elsewhere.' GO TO 20 END IF IF (UTASK) THEN CALL RKSIT(ASK,'UT',STATE1) IF (STATE1.NE.MINUS2) THEN IER = 911 NREC = 2 WRITE (REC,'(A/A)') &' ** You have called INTRP after you specified to SETUP ', &' ** that you were going to use UT. This is not permitted.' GO TO 20 END IF END IF IF (STATE.EQ.MINUS1) THEN IER = 911 NREC = 1 WRITE (REC,'(A)') &' ** You have not called CT, so you cannot use INTRP.' GO TO 20 END IF IF (STATE.GT.PLUS1) THEN IER = 911 NREC = 2 WRITE (REC,'(A/A)') &' ** CT has returned with a flag value greater than 1.', &' ** You cannot call INTRP in this circumstance.' GO TO 20 END IF C C Check input C REQST1 = REQEST(1:1) LEGALR = REQST1 .EQ. 'S' .OR. REQST1 .EQ. 's' .OR. & REQST1 .EQ. 'D' .OR. REQST1 .EQ. 'd' .OR. & REQST1 .EQ. 'B' .OR. REQST1 .EQ. 'b' IF (.NOT.LEGALR) THEN IER = 911 NREC = 3 WRITE (REC,'(A/A,A,A/A)') &' ** You have set the first character of ', &' ** REQEST to be ''',REQST1,'''. It must be one of ', &' ** ''S'',''s'',''D'',''d'',''B'' or ''b''.' GO TO 20 END IF C IF (NWANT.GT.NEQN) THEN IER = 911 NREC = 3 WRITE (REC,'(A,I6,A/A,I6,A/A)') &' ** You have specified the value of NWANT to be ',NWANT,'. This', &' ** is greater than ',NEQN,', which is the number of equations ', &' ** in the system being integrated.' GO TO 20 ELSE IF (NWANT.LT.1) THEN IER = 911 NREC = 3 WRITE (REC,'(A,I6,A/A/A)') &' ** You have specified the value of NWANT to be ',NWANT,', but ', &' ** this is less than 1. You cannot interpolate a zero or ', &' ** negative number of components.' GO TO 20 END IF C IF (METHD.EQ.1) THEN IF (LENINT.LT.1) THEN IER = 911 NREC = 2 WRITE (REC,'(A,I6,A/A)') &' ** You have specified LENINT to be ',LENINT,'.', &' ** This is too small. LENINT must be at least 1.' GO TO 20 END IF STARTP = 1 ELSE IF (METHD.EQ.2) THEN ICHK = NEQN + MAX(NEQN,5*NWANT) IF (LENINT.LT.ICHK) THEN IER = 911 NREC = 3 WRITE (REC,'(A,I6,A/A/A,I6,A)') &' ** You have specified LENINT to be ',LENINT,'. This is too', &' ** small. NINT must be at least NEQ + MAX(NEQ, 5*NWANT) ', &' ** which is ', ICHK,'.' GO TO 20 END IF STARTP = NEQN + 1 ELSE IF (METHD.EQ.3) THEN IER = 911 NREC = 5 WRITE (REC,'(A/A/A/A/A)') &' ** You have been using CT with METHOD = 3 to integrate your ', &' ** equations. You have just called INTRP, but interpolation ', &' ** is not available for this METHOD. Either use METHOD = 2, ', &' ** for which interpolation is available, or use RESET to make', &' ** CT step exactly to the points where you want output.' GO TO 20 END IF C C Has the interpolant been initialised for this step? C CALL RKSIT(ASK,SRNAME,STATE) ININTP = STATE .NE. MINUS2 C C Some initialization must be done before interpolation is possible. C To reduce the overhead, the interpolating polynomial is formed for C the first NWANT components. In the unusual circumstance that NWANT C is changed while still interpolating within the span of the current C step, the scheme must be reinitialized to accomodate the additional C components. C IF (.NOT.ININTP .OR. NWANT.NE.NWNTSV) THEN C C At present the derivative of the solution at the previous step, YPOLD(*), C is stored in the scratch array area starting at PRSCR. In the case of C METHD = 1 we can overwrite the stages. C IF (METHD.EQ.1) THEN CALL FORMI(F,NEQN,NWANT,WORK(PRY),WORK(PRYP),WORK(PRYOLD), & WORK(PRSCR),WORK(PRSTGS),.NOT.ININTP, & WORK(PRSTGS),WORK(PRSTGS)) ELSE CALL FORMI(F,NEQN,NWANT,WORK(PRY),WORK(PRYP),WORK(PRYOLD), & WORK(PRSCR),WORK(PRSTGS),.NOT.ININTP,WRKINT, & WRKINT(STARTP)) END IF C C Set markers to show that interpolation has been initialized for C NWANT components. NWNTSV = NWANT ININTP = .TRUE. END IF C C The actual evaluation of the interpolating polynomial and/or its first C derivative is done in EVALI. C IF (METHD.EQ.1) THEN CALL EVALI(WORK(PRY),WORK(PRYP),WORK(PRSTGS),TWANT,REQEST, & NWANT,YWANT,YPWANT) ELSE CALL EVALI(WORK(PRY),WORK(PRYP),WRKINT(STARTP),TWANT,REQEST, & NWANT,YWANT,YPWANT) END IF C 20 CONTINUE C CALL RKMSG(IER,SRNAME,NREC,FLAG) C RETURN END SUBROUTINE INTRP C SUBROUTINE RESET(TENDNU) C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C If you are not familiar with the code RESET and how it is used in C conjunction with CT to solve initial value problems, you should study the C document file rksuite.doc carefully before attempting to use the code. The C following "Brief Reminder" is intended only to remind you of the meaning, C type, and size requirements of the arguments. C C The integration using CT proceeds from TSTART in the direction of TEND, and C is now at TNOW. To reset TEND to a new value TENDNU, just call RESET with C TENDNU as the argument. You must continue integrating in the same C direction, so the sign of (TENDNU - TNOW) must be the same as the sign of C (TEND - TSTART). To change direction you must restart by a call to SETUP. C C INPUT VARIABLE C C TENDNU - DOUBLE PRECISION C The new value of TEND. C Constraint: TENDNU and TNOW must be clearly distinguishable C in the precision used. The sign of (TENDNU - TNOW) must be C the same as the sign of (TEND - TSTART). C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C .. Scalar Arguments .. DOUBLE PRECISION TENDNU C .. Common Block for Problem Definition .. DOUBLE PRECISION TSTRT, TND, DIR, HSTRT, TOLR INTEGER NEQN COMMON /RKCOM1/ TSTRT, TND, DIR, HSTRT, TOLR, NEQN SAVE /RKCOM1/ C .. Common Block to hold Problem Status .. DOUBLE PRECISION T, H, TOLD, HOLD INTEGER NFCN, SVNFCN, OKSTP, FLSTP LOGICAL FIRST, LAST COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP, & FIRST, LAST SAVE /RKCOM2/ C .. Common Block to hold Formula Characterisitcs .. DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4 INTEGER ORDER, LSTSTG, MAXTRY, NSEC LOGICAL FSAL COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY, & NSEC, FSAL SAVE /RKCOM5/ C .. Common Block for Environment Parameters .. DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY INTEGER OUTCH COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY, & OUTCH SAVE /RKCOM7/ C .. Common Block for Integrator Options .. LOGICAL MSG, UTASK COMMON /RKCOM8/ MSG, UTASK SAVE /RKCOM8/ C .. Common Block for Error Message .. CHARACTER*80 REC(10) COMMON /RKCOM9/ REC SAVE /RKCOM9/ C .. Parameters .. CHARACTER*6 SRNAME PARAMETER (SRNAME='RESET') LOGICAL ASK INTEGER MINUS1, MINUS2 PARAMETER (ASK=.TRUE.,MINUS1=-1,MINUS2=-2) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D+0) C .. Local Scalars .. DOUBLE PRECISION HMIN, TDIFF INTEGER FLAG, IER, NREC, STATE, STATE1 C .. External Subroutines .. EXTERNAL RKMSG, RKSIT C .. Intrinsic Functions .. INTRINSIC ABS, MAX C .. Executable Statements .. IER = 1 NREC = 0 C C Is it permissible to call RESET? C CALL RKSIT(ASK,'CT',STATE) IF (STATE.EQ.911) THEN IER = 912 NREC = 1 WRITE (REC,'(A)') &' ** A catastrophic error has already been detected elsewhere.' GO TO 20 END IF IF (UTASK) THEN CALL RKSIT(ASK,'UT',STATE1) IF (STATE1.NE.MINUS2) THEN IER = 911 NREC = 2 WRITE (REC,'(A/A)') &' ** You have called RESET after you specified to SETUP that ', &' ** you were going to use UT. This is not permitted.' GO TO 20 END IF END IF IF (STATE.EQ.MINUS1) THEN IER = 911 NREC = 1 WRITE (REC,'(A)') &' ** You have not called CT, so you cannot use RESET.' GO TO 20 END IF IF (STATE.EQ.5 .OR. STATE.EQ.6) THEN IER = 911 NREC = 2 WRITE (REC,'(A,I1,A/A)') &' ** CT has returned with CFLAG = ',STATE,'.', &' ** You cannot call RESET in this circumstance.' GO TO 20 END IF C C Check value of TENDNU C IF (DIR.GT.ZERO .AND. TENDNU.LE.T) THEN IER = 911 NREC = 4 WRITE (REC,'(A/A,D13.5/A,D13.5,A/A)') &' ** Integration is proceeding in the positive direction. The ', &' ** current value for the independent variable is ',T, &' ** and you have set TENDNU = ',TENDNU,'. TENDNU must be ', &' ** greater than T.' ELSE IF (DIR.LT.ZERO .AND. TENDNU.GE.T) THEN IER = 911 NREC = 4 WRITE (REC,'(A/A,D13.5/A,D13.5,A/A)') &' ** Integration is proceeding in the negative direction. The ', &' ** current value for the independent variable is ',T, &' ** and you have set TENDNU = ',TENDNU,'. TENDNU must be ', &' ** less than T.' ELSE HMIN = MAX(TINY,TOOSML*MAX(ABS(T),ABS(TENDNU))) TDIFF = ABS(TENDNU-T) IF (TDIFF.LT.HMIN) THEN IER = 911 NREC = 4 WRITE (REC,'(A,D13.5,A/A,D13.5,A/A/A,D13.5,A)') &' ** The current value of the independent variable T is ',T,'.', &' ** The TENDNU you supplied has ABS(TENDNU-T) = ',TDIFF,'.', &' ** For the METHOD and the precision of the computer being ', &' ** used, this difference must be at least ',HMIN,'.' END IF END IF IF (IER.EQ.911) GO TO 20 C TND = TENDNU LAST = .FALSE. C 20 CONTINUE C CALL RKMSG(IER,SRNAME,NREC,FLAG) C RETURN END SUBROUTINE RESET C SUBROUTINE MCONST(METHOD) C************************************************ C**** NOT A DESIGNATED USER-CALLABLE ROUTINE **** C************************************************ C C Purpose: Sets machine-dependent global quantities C C Common: Initializes: /RKCOM7/ OUTCH, MCHEPS, DWARF, RNDOFF, C SQRRMC, CUBRMC, TINY C Reads: none C Alters: none C C Comments: C ========= C OUTCH, MCHEPS, DWARF are pure environmental parameters with values C obtained from a call to ENVIRN. The other quantities set depend on C the environmental parameters, the implementation, and, possibly, C METHOD. At present the METHODs implemented in the RK suite do not C influence the values of these quantities. C OUTCH - Standard output channel C MCHEPS - Largest positive number such that 1.0D0 + MCHEPS = 1.0D0 C DWARF - Smallest positive number C RNDOFF - 10 times MCHEPS C SQRRMC - Square root of MCHEPS C CUBRMC - Cube root of MCHEPS C TINY - Square root of DWARF C C .. Scalar Arguments .. INTEGER METHOD C .. Common Block for Environment Parameters .. DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY INTEGER OUTCH COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY, & OUTCH SAVE /RKCOM7/ C .. Parameters .. DOUBLE PRECISION TEN, THIRD PARAMETER (TEN=10.0D+0,THIRD=1.0D+0/3.0D+0) C .. External Subroutines .. EXTERNAL ENVIRN C .. Intrinsic Functions .. INTRINSIC SQRT C .. Executable Statements .. C CALL ENVIRN(OUTCH,MCHEPS,DWARF) C RNDOFF = TEN*MCHEPS SQRRMC = SQRT(MCHEPS) CUBRMC = MCHEPS**THIRD TINY = SQRT(DWARF) C RETURN END SUBROUTINE MCONST C SUBROUTINE ENVIRN(OUTCH,MCHEPS,DWARF) C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C The RK suite requires some environmental parameters that are provided by C this subroutine. The values provided with the distribution codes are those C appropriate to the IEEE standard. They must be altered, if necessary, to C those appropriate to the computing system you are using before calling the C codes of the suite. C C ================================================================ C ================================================================ C TO MAKE SURE THAT THESE MACHINE AND INSTALLATION DEPENDENT C QUANTITIES ARE SPECIFIED PROPERLY, THE DISTRIBUTION VERSION C WRITES A MESSAGE ABOUT THE MATTER TO THE STANDARD OUTPUT CHANNEL C AND TERMINATES THE RUN. THE VALUES PROVIDED IN THE DISTRIBUTION C VERSION SHOULD BE ALTERED, IF NECESSARY, AND THE "WRITE" AND C "STOP" STATEMENTS COMMENTED OUT. C ================================================================ C ================================================================ C C OUTPUT VARIABLES C C OUTCH - INTEGER C Standard output channel C MCHEPS - DOUBLE PRECISION C MCHEPS is the largest positive number such that C 1.0D0 + MCHEPS = 1.0D0. C DWARF - DOUBLE PRECISION C DWARF is the smallest positive number. C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C .. Scalar Arguments .. INTEGER OUTCH DOUBLE PRECISION DWARF, MCHEPS C .. Executable Statements .. C C The following six statements are to be Commented out after verification that C the machine and installation dependent quantities are specified correctly. C If you pass copies of RKSUITE on to others, please give them the whole C distribution version of RKSUITE, and in particular, give them a version C of ENVIRN that does not have the following six statements Commented out. C WRITE(*,*) ' Before using RKSUITE, you must verify that the ' C WRITE(*,*) ' machine- and installation-dependent quantities ' C WRITE(*,*) ' specified in the subroutine ENVIRN are correct, ' C WRITE(*,*) ' and then Comment these WRITE statements and the ' C WRITE(*,*) ' STOP statement out of ENVIRN. ' C STOP C C The following values are appropriate to IEEE arithmetic with the typical C standard output channel. C OUTCH = 6 MCHEPS = 1.11D-16 DWARF = 2.23D-308 C C------------------------------------------------------------------------------ C If you have the routines D1MACH and I1MACH on your system, you could C replace the preceding statements by the following ones to obtain the C appropriate machine dependent numbers. The routines D1MACH and I1MACH C are public domain software. They are available from NETLIB. C .. Scalar Arguments .. C INTEGER OUTCH C DOUBLE PRECISION DWARF, MCHEPS C .. External Functions .. C INTEGER I1MACH C DOUBLE PRECISION D1MACH C .. Executable Statements .. C C OUTCH = I1MACH(2) C MCHEPS = D1MACH(3) C DWARF = D1MACH(1) C C If you have the NAG Fortran Library available on your system, you could C replace the preceding statements by the following ones to obtain the C appropriate machine dependent numbers. C C .. Scalar Arguments .. C INTEGER OUTCH C DOUBLE PRECISION DWARF, MCHEPS C .. External Functions .. C DOUBLE PRECISION X02AJF, X02AMF C .. Executable Statements .. C C CALL X04AAF(0,OUTCH) C MCHEPS = X02AJF() C DWARF = X02AMF() C C If you have the IMSL MATH/LIBRARY available on your system, you could C replace the preceding statements by the following ones to obtain the C appropriate machine dependent numbers. C C .. Scalar Arguments .. C INTEGER OUTCH C DOUBLE PRECISION DWARF, MCHEPS C .. External Functions .. C DOUBLE PRECISION DMACH C .. Executable Statements .. C C CALL UMACH(2,OUTCH) C MCHEPS = DMACH(4) C DWARF = DMACH(1) C------------------------------------------------------------------------------ C RETURN END SUBROUTINE ENVIRN C SUBROUTINE CONST(METHOD,VECSTG,REQSTG,LINTPL) C************************************************ C**** NOT A DESIGNATED USER-CALLABLE ROUTINE **** ************************************************* C C Purpose: Set formula definitions and formula characteristics for C selected method. Return storage requirements for the C selected method. C C Input: METHOD C Output: VECSTG, REQSTG, LINTPL C C Common: Initializes: /RKCOM4/ A(*,*), B(*), C(*), BHAT(*), R(*), C E(*), PTR(*), NSTAGE, METHD, INTP, MINTP C /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, C TANANG, RS, RS1, RS2, RS3, RS4, ORDER, C LSTSTG, MAXTRY, NSEC, FSAL C Reads: /RKCOM7/ RNDOFF C Alters: none C C Comments: C ========= C Runge-Kutta formula pairs are described by a set of coefficients C and by the setting of a number of parameters that describe the C characteristics of the pair. The input variable METHD indicates C which of the three pairs available is to be set. In the case of C METHD = 2 additional coefficients are defined that make interpolation C of the results possible and provide an additional error estimator. C VECSTG is the number of columns of workspace required to compute the C stages of a METHD. For interpolation purposes the routine returns via C COMMON the logical variable INTP indicating whether interpolation is C possible and via the call list: C REQSTG - whether the stages are required to form the C interpolant (set .FALSE. if INTP=.FALSE.) C LINTPL - the number of extra columns of storage required for use C with UT (set 0 if INTP=.FALSE.) C C Quantities set in common blocks: C METHD - copy of METHOD C A, B, C, BHAT - coefficients of the selected method C R - extra coefficents for interpolation with METHD = 2 C E - extra coefficients for additional local error estimate C with METHD = 2 C PTR - vector of pointers indicating how individual stages are to C be stored. With it zero coefficients of the formulas can C be exploited to reduce the storage required C NSTAGE - number of stages for the specified METHD C INTP - indicates whether there is an associated interpolant C (depending on the method, the user may have to supply C extra workspace) C MINTP - the degree of the interpolating polynomial, if one exists C FSAL - indicates whether the last stage of a step can be used as C the first stage of the following step C LSTSTG - pointer to location of last stage for use with FSAL=.TRUE. C ORDER - the lower order of the pair of Runge-Kutta formulas that C constitute a METHD C TANANG, C STBRAD - the stability region of the formula used to advance C the integration is approximated by a sector in the left half C complex plane. TANANG is the tangent of the interior angle C of the sector and STBRAD is the radius of the sector. C COST - cost of a successful step in function evaluations C MAXTRY - limit on the number of iterations in the stiffness check. As C set, no more than 24 function evaluations are made in the check. C NSEC - each step of size H in the primary integration corresponds to C NSEC steps of size H/NSEC in the secondary integration when C global error assessment is done. C EXPON - used to adjust the step size; this code implements an error C per step control for which EXPON = 1/(ORDER + 1). C SAFETY - quantity used in selecting the step size C TOOSML - quantity used to determine when a step size is too small for C the precision available C RS, RS1, C RS2, RS3, C RS4 - quantities used in determining the maximum and minimum change C change in step size (set independently of METHD) C C Further comments on SAFETY: C ============================ C The code estimates the largest step size that will yield the specified C accuracy. About half the time this step size would result in a local C error that is a little too large, and the step would be rejected. For C this reason a SAFETY factor is used to reduce the "optimal" value to one C more likely to succeed. Unfortunately, there is some art in choosing this C value. The more expensive a failed step, the smaller one is inclined to C take this factor. However, a small factor means that if the prediction were C good, more accuracy than desired would be obtained and the behavior of the C error would then be irregular. The more stringent the tolerance, the better C the prediction, except near limiting precision. Thus the general range of C tolerances expected influences the choice of SAFETY. C C .. Scalar Arguments .. INTEGER LINTPL, METHOD, VECSTG LOGICAL REQSTG C .. Common Block to hold Formula Definitions .. DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6), & E(7) INTEGER PTR(13), NSTAGE, METHD, MINTP LOGICAL INTP COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD, & MINTP, INTP SAVE /RKCOM4/ C .. Common Block to hold Formula Characterisitcs .. DOUBLE PRECISION TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4 INTEGER ORDER, LSTSTG, MAXTRY, NSEC LOGICAL FSAL COMMON /RKCOM5/ TOOSML, COST, SAFETY, EXPON, STBRAD, TANANG, & RS, RS1, RS2, RS3, RS4, ORDER, LSTSTG, MAXTRY, & NSEC, FSAL SAVE /RKCOM5/ C .. Common Block for Environment Parameters .. DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY INTEGER OUTCH COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY, & OUTCH SAVE /RKCOM7/ C .. Parameters .. DOUBLE PRECISION ONE, ZERO, TWO, FIFTY, FIVEPC PARAMETER (ONE=1.0D+0,ZERO=0.0D+0,TWO=2.0D+0,FIFTY=50.D+0, & FIVEPC=0.05D+0) C .. Local Scalars .. DOUBLE PRECISION CDIFF, DIFF INTEGER I, J C .. Intrinsic Functions .. INTRINSIC ABS, DBLE, INT, MAX, MIN C .. Executable Statements .. C METHD = METHOD C GO TO (20,40,100) METHD C C METHD = 1. C This pair is from "A 3(2) Pair of Runge-Kutta Formulas" by P. Bogacki C and L.F. Shampine, Appl. Math. Lett., 2, pp. 321-325, 1989. The authors C are grateful to P. Bogacki for his assistance in implementing the pair. C 20 CONTINUE NSTAGE = 4 FSAL = .TRUE. ORDER = 2 TANANG = 8.9D0 STBRAD = 2.3D0 SAFETY = 0.8D0 INTP = .TRUE. MINTP = 3 REQSTG = .FALSE. LINTPL = 2 NSEC = 3 C PTR(1) = 0 PTR(2) = 1 PTR(3) = 2 PTR(4) = 3 C A(2,1) = 1.D0/2.D0 A(3,1) = 0.D0 A(3,2) = 3.D0/4.D0 A(4,1) = 2.D0/9.D0 A(4,2) = 1.D0/3.D0 A(4,3) = 4.D0/9.D0 C C The coefficients BHAT(*) refer to the formula used to advance the C integration, here the one of order 3. The coefficients B(*) refer C to the other formula, here the one of order 2. For this pair, BHAT(*) C is not needed since FSAL = .TRUE. C B(1) = 7.D0/24.D0 B(2) = 1.D0/4.D0 B(3) = 1.D0/3.D0 B(4) = 1.D0/8.D0 C C(1) = 0.D0 C(2) = 1.D0/2.D0 C(3) = 3.D0/4.D0 C(4) = 1.D0 C GO TO 120 C C METHD = 2 C This pair is from "An Efficient Runge-Kutta (4,5) Pair" by P. Bogacki C and L.F. Shampine, Rept. 89-20, Math. Dept., Southern Methodist C University, Dallas, Texas, USA, 1989. The authors are grateful to C P. Bogacki for his assistance in implementing the pair. Shampine and C Bogacki subsequently modified the formula to enhance the reliability of C the pair. The original fourth order formula is used in an estimate of C the local error. If the step fails, the computation is broken off. If C the step is acceptable, the first evaluation of the next step is done, C i.e., the pair is implemented as FSAL and the local error of the step C is again estimated with a fourth order formula using the additional data. C The step must succeed with both estimators to be accepted. When the C second estimate is formed, it is used for the subsequent adjustment of C the step size because it is of higher quality. The two fourth order C formulas are well matched to leading order, and only exceptionally do C the estimators disagree -- problems with discontinuous coefficients are C handled more reliably by using two estimators as is global error C estimation. C 40 CONTINUE NSTAGE = 8 FSAL = .TRUE. ORDER = 4 TANANG = 5.2D0 STBRAD = 3.9D0 SAFETY = 0.8D0 INTP = .TRUE. REQSTG = .TRUE. MINTP = 6 LINTPL = 6 NSEC = 2 C PTR(1) = 0 PTR(2) = 1 PTR(3) = 2 PTR(4) = 3 PTR(5) = 4 PTR(6) = 5 PTR(7) = 6 PTR(8) = 7 C A(2,1) = 1.D0/6.D0 A(3,1) = 2.D0/27.D0 A(3,2) = 4.D0/27.D0 A(4,1) = 183.D0/1372.D0 A(4,2) = -162.D0/343.D0 A(4,3) = 1053.D0/1372.D0 A(5,1) = 68.D0/297.D0 A(5,2) = -4.D0/11.D0 A(5,3) = 42.D0/143.D0 A(5,4) = 1960.D0/3861.D0 A(6,1) = 597.D0/22528.D0 A(6,2) = 81.D0/352.D0 A(6,3) = 63099.D0/585728.D0 A(6,4) = 58653.D0/366080.D0 A(6,5) = 4617.D0/20480.D0 A(7,1) = 174197.D0/959244.D0 A(7,2) = -30942.D0/79937.D0 A(7,3) = 8152137.D0/19744439.D0 A(7,4) = 666106.D0/1039181.D0 A(7,5) = -29421.D0/29068.D0 A(7,6) = 482048.D0/414219.D0 A(8,1) = 587.D0/8064.D0 A(8,2) = 0.D0 A(8,3) = 4440339.D0/15491840.D0 A(8,4) = 24353.D0/124800.D0 A(8,5) = 387.D0/44800.D0 A(8,6) = 2152.D0/5985.D0 A(8,7) = 7267.D0/94080.D0 C C The coefficients B(*) refer to the formula of order 4. C B(1) = 2479.D0/34992.D0 B(2) = 0.D0 B(3) = 123.D0/416.D0 B(4) = 612941.D0/3411720.D0 B(5) = 43.D0/1440.D0 B(6) = 2272.D0/6561.D0 B(7) = 79937.D0/1113912.D0 B(8) = 3293.D0/556956.D0 C C The coefficients E(*) refer to an estimate of the local error based on C the first formula of order 4. It is the difference of the fifth order C result, here located in A(8,*), and the fourth order result. By C construction both E(2) and E(7) are zero. C E(1) = -3.D0/1280.D0 E(2) = 0.D0 E(3) = 6561.D0/632320.D0 E(4) = -343.D0/20800.D0 E(5) = 243.D0/12800.D0 E(6) = -1.D0/95.D0 E(7) = 0.D0 C C(1) = 0.D0 C(2) = 1.D0/6.D0 C(3) = 2.D0/9.D0 C(4) = 3.D0/7.D0 C(5) = 2.D0/3.D0 C(6) = 3.D0/4.D0 C(7) = 1.D0 C(8) = 1.D0 C C To do interpolation with this pair, some extra stages have to be computed. C The following additional A(*,*) and C(*) coefficients are for this purpose. C In addition there is an array R(*,*) that plays a role for interpolation C analogous to that of BHAT(*) for the basic step. C C(9) = 1.D0/2.D0 C(10) = 5.D0/6.D0 C(11) = 1.D0/9.D0 C A(9,1) = 455.D0/6144.D0 A(10,1) = -837888343715.D0/13176988637184.D0 A(11,1) = 98719073263.D0/1551965184000.D0 A(9,2) = 0.D0 A(10,2) = 30409415.D0/52955362.D0 A(11,2) = 1307.D0/123552.D0 A(9,3) = 10256301.D0/35409920.D0 A(10,3) = -48321525963.D0/759168069632.D0 A(11,3) = 4632066559387.D0/70181753241600.D0 A(9,4) = 2307361.D0/17971200.D0 A(10,4) = 8530738453321.D0/197654829557760.D0 A(11,4) = 7828594302389.D0/382182512025600.D0 A(9,5) = -387.D0/102400.D0 A(10,5) = 1361640523001.D0/1626788720640.D0 A(11,5) = 40763687.D0/11070259200.D0 A(9,6) = 73.D0/5130.D0 A(10,6) = -13143060689.D0/38604458898.D0 A(11,6) = 34872732407.D0/224610586200.D0 A(9,7) = -7267.D0/215040.D0 A(10,7) = 18700221969.D0/379584034816.D0 A(11,7) = -2561897.D0/30105600.D0 A(9,8) = 1.D0/32.D0 A(10,8) = -5831595.D0/847285792.D0 A(11,8) = 1.D0/10.D0 A(10,9) = -5183640.D0/26477681.D0 A(11,9) = -1.D0/10.D0 A(11,10) = -1403317093.D0/11371610250.D0 C DO 60 I = 1, 11 R(I,1) = 0.D0 60 CONTINUE DO 80 I = 1, 6 R(2,I) = 0.D0 80 CONTINUE R(1,6) = -12134338393.D0/1050809760.D0 R(1,5) = -1620741229.D0/50038560.D0 R(1,4) = -2048058893.D0/59875200.D0 R(1,3) = -87098480009.D0/5254048800.D0 R(1,2) = -11513270273.D0/3502699200.D0 C R(3,6) = -33197340367.D0/1218433216.D0 R(3,5) = -539868024987.D0/6092166080.D0 R(3,4) = -39991188681.D0/374902528.D0 R(3,3) = -69509738227.D0/1218433216.D0 R(3,2) = -29327744613.D0/2436866432.D0 C R(4,6) = -284800997201.D0/19905339168.D0 R(4,5) = -7896875450471.D0/165877826400.D0 R(4,4) = -333945812879.D0/5671036800.D0 R(4,3) = -16209923456237.D0/497633479200.D0 R(4,2) = -2382590741699.D0/331755652800.D0 C R(5,6) = -540919.D0/741312.D0 R(5,5) = -103626067.D0/43243200.D0 R(5,4) = -633779.D0/211200.D0 R(5,3) = -32406787.D0/18532800.D0 R(5,2) = -36591193.D0/86486400.D0 C R(6,6) = 7157998304.D0/374350977.D0 R(6,5) = 30405842464.D0/623918295.D0 R(6,4) = 183022264.D0/5332635.D0 R(6,3) = -3357024032.D0/1871754885.D0 R(6,2) = -611586736.D0/89131185.D0 C R(7,6) = -138073.D0/9408.D0 R(7,5) = -719433.D0/15680.D0 R(7,4) = -1620541.D0/31360.D0 R(7,3) = -385151.D0/15680.D0 R(7,2) = -65403.D0/15680.D0 C R(8,6) = 1245.D0/64.D0 R(8,5) = 3991.D0/64.D0 R(8,4) = 4715.D0/64.D0 R(8,3) = 2501.D0/64.D0 R(8,2) = 149.D0/16.D0 R(8,1) = 1.D0 C R(9,6) = 55.D0/3.D0 R(9,5) = 71.D0 R(9,4) = 103.D0 R(9,3) = 199.D0/3.D0 R(9,2) = 16.0D0 C R(10,6) = -1774004627.D0/75810735.D0 R(10,5) = -1774004627.D0/25270245.D0 R(10,4) = -26477681.D0/359975.D0 R(10,3) = -11411880511.D0/379053675.D0 R(10,2) = -423642896.D0/126351225.D0 C R(11,6) = 35.D0 R(11,5) = 105.D0 R(11,4) = 117.D0 R(11,3) = 59.D0 R(11,2) = 12.D0 C GO TO 120 C C METHD = 3 C This pair is from "High Order Embedded Runge-Kutta Formulae" by P.J. C Prince and J.R. Dormand, J. Comp. Appl. Math.,7, pp. 67-75, 1981. The C authors are grateful to P. Prince and J. Dormand for their assistance in C implementing the pair. C 100 CONTINUE NSTAGE = 13 FSAL = .FALSE. ORDER = 7 TANANG = 11.0D0 STBRAD = 5.2D0 SAFETY = 0.8D0 INTP = .FALSE. REQSTG = .FALSE. MINTP = 0 LINTPL = 0 NSEC = 2 C PTR(1) = 0 PTR(2) = 1 PTR(3) = 2 PTR(4) = 1 PTR(5) = 3 PTR(6) = 2 PTR(7) = 4 PTR(8) = 5 PTR(9) = 6 PTR(10) = 7 PTR(11) = 8 PTR(12) = 9 PTR(13) = 1 C A(2,1) = 5.55555555555555555555555555556D-2 A(3,1) = 2.08333333333333333333333333333D-2 A(3,2) = 6.25D-2 A(4,1) = 3.125D-2 A(4,2) = 0.D0 A(4,3) = 9.375D-2 A(5,1) = 3.125D-1 A(5,2) = 0.D0 A(5,3) = -1.171875D0 A(5,4) = 1.171875D0 A(6,1) = 3.75D-2 A(6,2) = 0.D0 A(6,3) = 0.D0 A(6,4) = 1.875D-1 A(6,5) = 1.5D-1 A(7,1) = 4.79101371111111111111111111111D-2 A(7,2) = 0.D0 A(7,3) = 0.0D0 A(7,4) = 1.12248712777777777777777777778D-1 A(7,5) = -2.55056737777777777777777777778D-2 A(7,6) = 1.28468238888888888888888888889D-2 A(8,1) = 1.6917989787292281181431107136D-2 A(8,2) = 0.D0 A(8,3) = 0.D0 A(8,4) = 3.87848278486043169526545744159D-1 A(8,5) = 3.59773698515003278967008896348D-2 A(8,6) = 1.96970214215666060156715256072D-1 A(8,7) = -1.72713852340501838761392997002D-1 A(9,1) = 6.90957533591923006485645489846D-2 A(9,2) = 0.D0 A(9,3) = 0.D0 A(9,4) = -6.34247976728854151882807874972D-1 A(9,5) = -1.61197575224604080366876923982D-1 A(9,6) = 1.38650309458825255419866950133D-1 A(9,7) = 9.4092861403575626972423968413D-1 A(9,8) = 2.11636326481943981855372117132D-1 A(10,1) = 1.83556996839045385489806023537D-1 A(10,2) = 0.D0 A(10,3) = 0.D0 A(10,4) = -2.46876808431559245274431575997D0 A(10,5) = -2.91286887816300456388002572804D-1 A(10,6) = -2.6473020233117375688439799466D-2 A(10,7) = 2.84783876419280044916451825422D0 A(10,8) = 2.81387331469849792539403641827D-1 A(10,9) = 1.23744899863314657627030212664D-1 A(11,1) = -1.21542481739588805916051052503D0 A(11,2) = 0.D0 A(11,3) = 0.D0 A(11,4) = 1.66726086659457724322804132886D1 A(11,5) = 9.15741828416817960595718650451D-1 A(11,6) = -6.05660580435747094755450554309D0 A(11,7) = -1.60035735941561781118417064101D1 A(11,8) = 1.4849303086297662557545391898D1 A(11,9) = -1.33715757352898493182930413962D1 A(11,10) = 5.13418264817963793317325361166D0 A(12,1) = 2.58860916438264283815730932232D-1 A(12,2) = 0.D0 A(12,3) = 0.D0 A(12,4) = -4.77448578548920511231011750971D0 A(12,5) = -4.3509301377703250944070041181D-1 A(12,6) = -3.04948333207224150956051286631D0 A(12,7) = 5.57792003993609911742367663447D0 A(12,8) = 6.15583158986104009733868912669D0 A(12,9) = -5.06210458673693837007740643391D0 A(12,10) = 2.19392617318067906127491429047D0 A(12,11) = 1.34627998659334941535726237887D-1 A(13,1) = 8.22427599626507477963168204773D-1 A(13,2) = 0.D0 A(13,3) = 0.D0 A(13,4) = -1.16586732572776642839765530355D1 A(13,5) = -7.57622116690936195881116154088D-1 A(13,6) = 7.13973588159581527978269282765D-1 A(13,7) = 1.20757749868900567395661704486D1 A(13,8) = -2.12765911392040265639082085897D0 A(13,9) = 1.99016620704895541832807169835D0 A(13,10) = -2.34286471544040292660294691857D-1 A(13,11) = 1.7589857770794226507310510589D-1 A(13,12) = 0.D0 C C The coefficients BHAT(*) refer to the formula used to advance the C integration, here the one of order 8. The coefficients B(*) refer C to the other formula, here the one of order 7. C BHAT(1) = 4.17474911415302462220859284685D-2 BHAT(2) = 0.D0 BHAT(3) = 0.D0 BHAT(4) = 0.D0 BHAT(5) = 0.D0 BHAT(6) = -5.54523286112393089615218946547D-2 BHAT(7) = 2.39312807201180097046747354249D-1 BHAT(8) = 7.0351066940344302305804641089D-1 BHAT(9) = -7.59759613814460929884487677085D-1 BHAT(10) = 6.60563030922286341461378594838D-1 BHAT(11) = 1.58187482510123335529614838601D-1 BHAT(12) = -2.38109538752862804471863555306D-1 BHAT(13) = 2.5D-1 C B(1) = 2.9553213676353496981964883112D-2 B(2) = 0.D0 B(3) = 0.D0 B(4) = 0.D0 B(5) = 0.D0 B(6) = -8.28606276487797039766805612689D-1 B(7) = 3.11240900051118327929913751627D-1 B(8) = 2.46734519059988698196468570407D0 B(9) = -2.54694165184190873912738007542D0 B(10) = 1.44354858367677524030187495069D0 B(11) = 7.94155958811272872713019541622D-2 B(12) = 4.44444444444444444444444444445D-2 B(13) = 0.D0 C C(1) = 0.D0 C(2) = 5.55555555555555555555555555556D-2 C(3) = 8.33333333333333333333333333334D-2 C(4) = 1.25D-1 C(5) = 3.125D-1 C(6) = 3.75D-1 C(7) = 1.475D-1 C(8) = 4.65D-1 C(9) = 5.64865451382259575398358501426D-1 C(10) = 6.5D-1 C(11) = 9.24656277640504446745013574318D-1 C(12) = 1.D0 C(13) = C(12) C GO TO 120 C C The definitions of all pairs come here for the calculation of C LSTSTG, RS1, RS2, RS3, RS4, COST, MAXTRY, EXPON, TOOSML, and VECSTG. C 120 CONTINUE LSTSTG = PTR(NSTAGE) IF (FSAL) THEN COST = DBLE(NSTAGE-1) ELSE COST = DBLE(NSTAGE) END IF C C MAXTRY - limit on the number of iterations of a computation made in C diagnosing stiffness. There are at most Q = 3 function calls per C iteration. MAXTRY is determined so that Q*MAXTRY <= 5% of the cost of C 50 steps and 1 <= MAXTRY <= 8. This limits the number of calls to FCN C in each diagnosis of stiffness to 24 calls. C MAXTRY = MIN(8,MAX(1,INT(FIVEPC*COST*FIFTY))) C EXPON = ONE/(ORDER+ONE) C C In calculating CDIFF it is assumed that there will be a non-zero C difference |C(I) - C(J)| less than one. If C(I) = C(J) for any I not C equal to J, they should be made precisely equal by assignment. C CDIFF = ONE DO 160 I = 1, NSTAGE - 1 DO 140 J = I + 1, NSTAGE DIFF = ABS(C(I)-C(J)) IF (DIFF.NE.ZERO) CDIFF = MIN(CDIFF,DIFF) 140 CONTINUE 160 CONTINUE TOOSML = RNDOFF/CDIFF C C Determine the number of columns needed in STAGES(1:NEQ,*) (this will be C at most NSTAGE-1 since the first stage is held in a separate array). C The PTR array contains the column positions of the stages. C VECSTG = 0 DO 180 I = 2, NSTAGE VECSTG = MAX(PTR(I),VECSTG) 180 CONTINUE C RS = TWO RS1 = ONE/RS RS2 = RS**2 RS3 = RS**3 RS4 = ONE/RS3 C RETURN END SUBROUTINE CONST C SUBROUTINE FORMI(F,NEQ,NWANT,Y,YP,YOLD,YPOLD,STAGES,CALSTG, & XSTAGE,P) C************************************************ C**** NOT A DESIGNATED USER-CALLABLE ROUTINE **** C************************************************ C C Purpose: Forms an interpolating polynomial for use with C METHDs 1 or 2. C C Input: NEQ, NWANT, T, Y(*), YP(*), HOLD, YOLD(*), YPOLD(*), C STAGES(NEQ,*), CALSTG C Output: P(*), XSTAGE(NEQ) C External: F C C Common: Initializes: none C Reads: /RKCOM4/ A(*,*), C(*), R(*), METHD, MINTP C /RKCOM2/ T, TOLD, HOLD C Alters: /RKCOM2/ NFCN C C Comments: C ========= C The integration has reached T with a step HOLD from TOLD = T-HOLD. C Y(*),YP(*) and YOLD(*),YPOLD(*) approximate the solution and its C derivative at T and TOLD respectively. STAGES(NEQ,*) holds the stages C computed in taking this step. In the case of METHD = 2 it is necessary C to compute some more stages in this subroutine. CALSTG indicates whether C or not the extra stages need to be computed. A(*,*) and C(*) are used in C computing these stages. The extra stages are stored in STAGES(NEQ,*) and C XSTAGE(*). The coefficients of the interpolating polynomials for the first C NWANT components of the solution are returned in the array P(*). The C polynomial is of degree MINTP = 3 for METHD = 1 and of degree MINTP = 6 C for METHD = 2. The vector R(*) is used for workspace when METHD = 2. C C .. Scalar Arguments .. INTEGER NEQ, NWANT LOGICAL CALSTG C .. Array Arguments .. DOUBLE PRECISION P(*), STAGES(NEQ,*), XSTAGE(*), Y(*), YOLD(*), & YP(*), YPOLD(*) C .. Subroutine Arguments .. EXTERNAL F C .. Common Block to hold Problem Status .. DOUBLE PRECISION T, H, TOLD, HOLD INTEGER NFCN, SVNFCN, OKSTP, FLSTP LOGICAL FIRST, LAST COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP, & FIRST, LAST SAVE /RKCOM2/ C .. Common Block to hold Formula Definitions .. DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6), & E(7) INTEGER PTR(13), NSTAGE, METHD, MINTP LOGICAL INTP COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD, & MINTP, INTP SAVE /RKCOM4/ C .. Local Scalars .. DOUBLE PRECISION D1, D2, D3, D4, HYP, HYPOLD INTEGER I, J, K, L C .. Executable Statements .. C IF (METHD.EQ.1) THEN C C METHD = 1. Use the cubic Hermite interpolant that is is fully C specified by the values and slopes at the two ends of the step. C DO 20 L = 1, NWANT D1 = Y(L) - YOLD(L) HYP = HOLD*YP(L) HYPOLD = HOLD*YPOLD(L) D2 = HYP - D1 D3 = D1 - HYPOLD D4 = D2 - D3 P(L) = D2 + D4 P(NWANT+L) = D4 20 CONTINUE C ELSE C C METHD = 2. C IF (CALSTG) THEN C C Compute the extra stages needed for interpolation using the facts that C 1. Stage 1 is YPOLD(*). C 2. Stage i (i>1) is stored in STAGES(1:NEQ,i). C 3. This pair is FSAL, i.e. STAGES(1:NEQ,7)=YP(1:NEQ), which frees C up STAGES(1:NEQ,7) for use by stage 9. C 4. XSTAGE(1:NEQ) is used for stage 10. C 5. The coefficient of stage 2 in the interpolant is always 0, so C STAGES(1:NEQ,1) is used for stage 11. C The vector P(1:NEQ) is used as workspace for computing the stages. C DO 180 I = 9, 11 DO 40 L = 1, NEQ P(L) = A(I,1)*YPOLD(L) 40 CONTINUE DO 140 J = 2, I - 1 IF (J.LE.7) THEN DO 60 L = 1, NEQ P(L) = P(L) + A(I,J)*STAGES(L,J-1) 60 CONTINUE ELSE IF (J.EQ.8) THEN DO 80 L = 1, NEQ P(L) = P(L) + A(I,J)*YP(L) 80 CONTINUE ELSE IF (J.EQ.9) THEN DO 100 L = 1, NEQ P(L) = P(L) + A(I,J)*STAGES(L,7) 100 CONTINUE ELSE IF (J.EQ.10) THEN DO 120 L = 1, NEQ P(L) = P(L) + A(I,J)*XSTAGE(L) 120 CONTINUE END IF 140 CONTINUE DO 160 L = 1, NEQ P(L) = YOLD(L) + HOLD*P(L) 160 CONTINUE IF (I.EQ.9) THEN CALL F(TOLD+C(I)*HOLD,P,STAGES(1,7)) NFCN = NFCN + 1 ELSE IF (I.EQ.10) THEN CALL F(TOLD+C(I)*HOLD,P,XSTAGE) NFCN = NFCN + 1 ELSE CALL F(TOLD+C(I)*HOLD,P,STAGES(1,1)) NFCN = NFCN + 1 END IF 180 CONTINUE END IF C C Form the coefficients of the interpolating polynomial in its shifted C and scaled form. The transformation from the form in which the C polynomial is derived can be somewhat ill-conditioned. The terms C are grouped so as to minimize the errors of the transformation. C C Coefficient of SIGMA**6 K = 4*NWANT DO 200 L = 1, NWANT P(K+L) = R(5,6)*STAGES(L,4) + & ((R(10,6)*XSTAGE(L)+R(8,6)*YP(L))+ & (R(7,6)*STAGES(L,6)+R(6,6)*STAGES(L,5))) + & ((R(4,6)*STAGES(L,3)+R(9,6)*STAGES(L,7))+ & (R(3,6)*STAGES(L,2)+R(11,6)*STAGES(L,1))+ & R(1,6)*YPOLD(L)) 200 CONTINUE C C Coefficient of SIGMA**5 K = 3*NWANT DO 220 L = 1, NWANT P(K+L) = (R(10,5)*XSTAGE(L)+R(9,5)*STAGES(L,7)) + & ((R(7,5)*STAGES(L,6)+R(6,5)*STAGES(L,5))+ & R(5,5)*STAGES(L,4)) + ((R(4,5)*STAGES(L,3)+ & R(8,5)*YP(L))+(R(3,5)*STAGES(L,2)+R(11,5)* & STAGES(L,1))+R(1,5)*YPOLD(L)) 220 CONTINUE C C Coefficient of SIGMA**4 K = 2*NWANT DO 240 L = 1, NWANT P(K+L) = ((R(4,4)*STAGES(L,3)+R(8,4)*YP(L))+ & (R(7,4)*STAGES(L,6)+R(6,4)*STAGES(L,5))+ & R(5,4)*STAGES(L,4)) + ((R(10,4)*XSTAGE(L)+ & R(9,4)*STAGES(L,7))+(R(3,4)*STAGES(L,2)+ & R(11,4)*STAGES(L,1))+R(1,4)*YPOLD(L)) 240 CONTINUE C C Coefficient of SIGMA**3 K = NWANT DO 260 L = 1, NWANT P(K+L) = R(5,3)*STAGES(L,4) + R(6,3)*STAGES(L,5) + & ((R(3,3)*STAGES(L,2)+R(9,3)*STAGES(L,7))+ & (R(10,3)*XSTAGE(L)+R(8,3)*YP(L))+R(1,3)* & YPOLD(L))+((R(4,3)*STAGES(L,3)+R(11,3)* & STAGES(L,1))+R(7,3)*STAGES(L,6)) 260 CONTINUE C C Coefficient of SIGMA**2 C DO 280 L = 1, NWANT P(L) = R(5,2)*STAGES(L,4) + ((R(6,2)*STAGES(L,5)+ & R(8,2)*YP(L))+R(1,2)*YPOLD(L)) + & ((R(3,2)*STAGES(L,2)+R(9,2)*STAGES(L,7))+ & R(10,2)*XSTAGE(L)) + ((R(4,2)*STAGES(L,3)+ & R(11,2)*STAGES(L,1))+R(7,2)*STAGES(L,6)) 280 CONTINUE C C Scale all the coefficients by the step size. C DO 300 L = 1, NWANT*(MINTP-1) P(L) = HOLD*P(L) 300 CONTINUE C END IF C RETURN END SUBROUTINE FORMI C SUBROUTINE EVALI(Y,YP,P,TWANT,REQEST,NWANT,YWANT,YPWANT) C************************************************ C**** NOT A DESIGNATED USER-CALLABLE ROUTINE **** ************************************************* C C Purpose: Evaluation of an interpolating polynomial and/or its C first derivative. C C Input: Y(*), YP(*), P(NWANT,*), TWANT, REQEST, NWANT C Output: YWANT(*), YPWANT(*) C C Common: Initializes: none C Reads: /RKCOM2/ HOLD, T C /RKCOM4/ MINTP C Alters: none C C Comments: C ========= C The interpolant is evaluated at TWANT to approximate the solution, C YWANT, and/or its first derivative there, YPWANT. Only the first C NWANT components of the answer are computed. There are three cases C that are indicated by the first character of REQEST: C REQEST(1:1) = `S' or `s'- compute approximate `S'olution only. C = `D' or `d'- compute approximate first `D'erivative C of the solution only. C = `B' or `b'- compute `B'oth approximate solution and C first derivative. C The coefficents of the polynomial are contained in Y(*), YP(*) and C P(NWANT,*). C C .. Scalar Arguments .. DOUBLE PRECISION TWANT INTEGER NWANT CHARACTER*(*) REQEST C .. Array Arguments .. DOUBLE PRECISION P(NWANT,*), Y(*), YP(*), YPWANT(*), YWANT(*) C .. Common Block to hold Problem Status .. DOUBLE PRECISION T, H, TOLD, HOLD INTEGER NFCN, SVNFCN, OKSTP, FLSTP LOGICAL FIRST, LAST COMMON /RKCOM2/ T, H, TOLD, HOLD, NFCN, SVNFCN, OKSTP, FLSTP, & FIRST, LAST SAVE /RKCOM2/ C .. Common Block to hold Formula Definitions .. DOUBLE PRECISION A(13,13), B(13), C(13), BHAT(13), R(11,6), & E(7) INTEGER PTR(13), NSTAGE, METHD, MINTP LOGICAL INTP COMMON /RKCOM4/ A, B, C, BHAT, R, E, PTR, NSTAGE, METHD, & MINTP, INTP SAVE /RKCOM4/ C .. Local Scalars .. DOUBLE PRECISION SIGMA INTEGER K, L CHARACTER REQST1 C .. Executable Statements .. C C Evaluate the interpolating polynomial of degree MINTP in terms of the C shifted and scaled independent variable SIGMA. C SIGMA = (TWANT-T)/HOLD C REQST1 = REQEST(1:1) IF (REQST1.EQ.'S' .OR. REQST1.EQ.'s' .OR. & REQST1.EQ.'B' .OR. REQST1.EQ.'b') THEN C DO 20 L = 1, NWANT YWANT(L) = P(L,MINTP-1)*SIGMA 20 CONTINUE DO 60 K = MINTP - 2, 1, -1 DO 40 L = 1, NWANT YWANT(L) = (YWANT(L)+P(L,K))*SIGMA 40 CONTINUE 60 CONTINUE DO 80 L = 1, NWANT YWANT(L) = (YWANT(L)+HOLD*YP(L))*SIGMA + Y(L) 80 CONTINUE END IF C C Evaluate the derivative of the interpolating polynomial. C IF (REQST1.EQ.'D' .OR. REQST1.EQ.'d' .OR. & REQST1.EQ.'B' .OR. REQST1.EQ.'b') THEN C C The derivative of the interpolating polynomial with respect to TWANT C is the derivative with respect to S divided by HOLD. C DO 100 L = 1, NWANT YPWANT(L) = MINTP*P(L,MINTP-1)*SIGMA 100 CONTINUE DO 140 K = MINTP - 1, 2, -1 DO 120 L = 1, NWANT YPWANT(L) = (YPWANT(L)+K*P(L,K-1))*SIGMA 120 CONTINUE 140 CONTINUE DO 160 L = 1, NWANT YPWANT(L) = (YPWANT(L)+HOLD*YP(L))/HOLD 160 CONTINUE END IF C RETURN END SUBROUTINE EVALI C SUBROUTINE RKMSG(IER,SRNAME,NREC,FLAG) C************************************************ C**** NOT A DESIGNATED USER-CALLABLE ROUTINE **** C************************************************ C C Purpose: To process error messages and terminate the program C in the event of a "catastrophic" failure. C C Input: IER, SRNAME, NREC C Output: FLAG C C Common: Initializes: none C Reads: /RKCOM7/ OUTCH C /RKCOM8/ MSG, UTASK C /RKCOM9/ REC C Alters: none C C Comments: C ========= C The output variable FLAG is assigned the value of the input variable IER. C C IER = -2 reports a successful call of the subroutine SRNAME and C indicates that special action is to be taken elsewhere C in the suite. FLAG is set and a return is effected. C C IER = 1 reports a successful call of the subroutine SRNAME. FLAG C is set and a return is effected. C C 1 < IER < 911 and MSG = .TRUE.: a message of NREC records contained in C the array REC(*) is written to the standard output channel, C OUTCH. FLAG is set and a return is effected. C C IER = 911 reports a "catastrophic" error was detected in SRNAME. A C message is written to OUTCH regardless of the value of MSG and C normally the execution of the program is terminated. The C execution is not terminated if the error is the result of an C indirect call to CT, RESET, or INTRP through UT (UTASK = .TRUE.). C Termination can be prevented by using the subroutine SOFTFL. C C IER = 912 reports that a "catastrophic" error was detected earlier and C termination was prevented, but the user has failed to take C appropriate remedial action. Execution is terminated. C C .. Scalar Arguments .. INTEGER FLAG, IER, NREC CHARACTER*(*) SRNAME C .. Common Block for Environment Parameters .. DOUBLE PRECISION MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY INTEGER OUTCH COMMON /RKCOM7/ MCHEPS, DWARF, RNDOFF, SQRRMC, CUBRMC, TINY, & OUTCH SAVE /RKCOM7/ C .. Common Block for Integrator Options .. LOGICAL MSG, UTASK COMMON /RKCOM8/ MSG, UTASK SAVE /RKCOM8/ C .. Common Block for Error Message .. CHARACTER*80 REC(10) COMMON /RKCOM9/ REC SAVE /RKCOM9/ C .. Parameters .. INTEGER PLUS1 LOGICAL ASK, TELL PARAMETER (PLUS1=1,ASK=.TRUE.,TELL=.FALSE.) C .. Local Scalars .. INTEGER I LOGICAL BADERR, OK, ON, UTCALL C .. External Subroutines .. EXTERNAL CHKFL, RKSIT, SOFTFL C .. Executable Statements .. C C Check where the call came from - if it is an indirect call from UT, C the run is not STOPped. UTCALL = (SRNAME.EQ.'RESET' .OR. SRNAME.EQ.'CT' .OR. & SRNAME.EQ.'INTRP') .AND. UTASK C C Check if can continue with integrator. OK = (SRNAME.EQ.'CT' .OR. SRNAME.EQ.'UT') .AND. & (IER.EQ.2 .OR. IER.EQ.3 .OR. IER.EQ.4) C C Check if program termination has been overridden. CALL SOFTFL(ASK,ON) C IF ((MSG.AND.IER.GT.PLUS1) .OR. IER.GE.911) THEN WRITE (OUTCH,'(/A)') ' **' WRITE (OUTCH,'(A)') (REC(I),I=1,NREC) IF (IER.GE.911) THEN WRITE (OUTCH,'(A/A,A,A/A/)') &' **', &' ** Catastrophic error detected in ', SRNAME, '.', &' **' IF ((.NOT.(UTCALL.OR.ON).AND.IER.EQ.911) .OR. & IER.EQ.912) THEN WRITE (OUTCH,'(A/A/A)') &' **', &' ** Execution of your program is being terminated.', &' **' STOP END IF ELSE IF (OK) THEN WRITE (OUTCH,'(A/A,A,A,I2,A/A/A)') &' **', &' ** Warning from routine ', SRNAME, ' with flag set ',IER, '.', &' ** You can continue integrating this problem.', &' **' ELSE WRITE (OUTCH,'(A/A,A,A,I2,A/A/A)') &' **', &' ** Warning from routine ', SRNAME, ' with flag set ',IER, '.', &' ** You cannot continue integrating this problem.', &' **' END IF END IF DO 20 I = NREC + 1, 10 REC(I) = ' ' 20 CONTINUE FLAG = IER C C TELL RKSIT the status of the routine associated with SRNAME CALL RKSIT(TELL,SRNAME,FLAG) C C Indicate that a catastrophic error has been detected BADERR = FLAG .GE. 911 CALL CHKFL(TELL,BADERR) C RETURN C END SUBROUTINE RKMSG C SUBROUTINE RKSIT(ASK,SRNAME,STATE) C************************************************ C**** NOT A DESIGNATED USER-CALLABLE ROUTINE **** C************************************************ C C Purpose: To save or enquire about the status of each C subprogram in the suite. C C Input: ASK, SRNAME C Input/output: STATE C C C Comments: C ========= C SRNAME indicates which routine is of interest in the call to RKSIT. C C If ASK=.FALSE., then the value of STATE (which as used is the error C flag) for the routine SRNAME is saved internally. This value of STATE C is usually positive. There are two exceptions: C 1. SRNAME='SETUP' and STATE=-1 indicates a completely new problem, C so all SAVEd states are cleared. C 2. STATE=-2 is used by some routines in the suite to indicate C circumstances which require special action. C C If ASK=.TRUE., then RKSIT first checks to see if there were any C catastrophic errors, that is, a SAVEd state has value 911. This should C happen only when the user has overridden program termination in the event C of catastrophic failures from routines in the package but has failed to C take appropriate action. If this is the case, then RKSIT returns a value C of STATE = 911 which forces a termination of execution inside RKMSG. If C no catastrophic errors are flagged, then STATE returns the saved state C value for the routine specified by SRNAME. C C .. Scalar Arguments .. INTEGER STATE LOGICAL ASK CHARACTER*(*) SRNAME C .. Parameters .. INTEGER STATES, MINUS1 PARAMETER (STATES=7,MINUS1=-1) C .. Local Scalars .. INTEGER I, NAME C .. Local Arrays .. INTEGER SVSTA(STATES) C .. Save statement .. SAVE SVSTA C .. Data statements .. DATA SVSTA/STATES*MINUS1/ C .. Executable Statements .. C IF (SRNAME.EQ.'SETUP') THEN NAME = 1 ELSE IF (SRNAME.EQ.'UT') THEN NAME = 2 ELSE IF (SRNAME.EQ.'STAT') THEN NAME = 3 ELSE IF (SRNAME.EQ.'GLBERR') THEN NAME = 4 ELSE IF (SRNAME.EQ.'CT') THEN NAME = 5 ELSE IF (SRNAME.EQ.'INTRP') THEN NAME = 6 ELSE IF (SRNAME.EQ.'RESET') THEN NAME = 7 ELSE NAME = 0 END IF C C (Re)initialize if SETUP is telling RKSIT to do so. IF (.NOT.ASK .AND. NAME.EQ.1 .AND. STATE.EQ.MINUS1) THEN DO 20 I = 1, STATES SVSTA(I) = MINUS1 20 CONTINUE GO TO 60 END IF C C Check for 911 on exit from a previous call. IF (ASK) THEN DO 40 I = 1, STATES IF (SVSTA(I).EQ.911) THEN STATE = 911 GO TO 60 END IF 40 CONTINUE END IF C IF (ASK) THEN STATE = SVSTA(NAME) ELSE SVSTA(NAME) = STATE END IF C 60 CONTINUE C RETURN END SUBROUTINE RKSIT C SUBROUTINE TRUERR(F,NEQ,Y,TOL,WEIGHT,ZY,ZYP,ZERROR,ZYNEW,ZERRES, & ZSTAGE,IER) C************************************************ C**** NOT A DESIGNATED USER-CALLABLE ROUTINE **** C************************************************ C C Purpose: Compute a running RMS measure of the true (global) error C for a general Runge-Kutta pair. C C C Input: NEQ, Y(*), TOL, WEIGHT(*), C Input/output: ZY(*), ZYP(*), ZERROR(*) C Workspace: ZYNEW(*), ZERRES(*), ZSTAGE(NEQ,*) C Output: IER C External: F C C Common: Initializes: none C Reads: /RKCOM2/ T, HOLD C /RKCOM5/ TOOSML, ORDER, NSEC C /RKCOM7/ TINY C Alters: /RKCOM6/ MAXERR, LOCMAX, GNFCN C C Comments: C ========= C A secondary integration is performed using a fraction of the step size C of the primary integration. ZY(*) and ZYP(*) are the approximate solution C and first derivative of this secondary integration. ZERRES(*) contains the C error estimates for the secondary integration. ZYNEW(*) and ZSTAGE(*,*) are C workspace for taking a step. The error assessment is computed using the C difference of the primary and secondary solutions at the primary C integration points as an estimate of the true error there. The weights C used are those of the error test of the primary integration. This error C assessment is maintained in the vector ZERROR(*). MAXERR and LOCMAX C contain the maximum contribution to the assessment and its location, C respectively. The number of calls to F is counted by GNFCN. C C .. Scalar Arguments .. DOUBLE PRECISION TOL INTEGER IER, NEQ C .. Array Arguments .. DOUBLE PRECISION WEIGHT(*), Y(*), ZERRES(*), ZERROR(*), & ZSTAGE(NEQ,*), ZY(*), ZYNEW(*), ZYP(*) C .. Su