PROGRAM PRO5P1 C... C... THIS MAIN PROGRAM ESSENTIALLY CALLS SUBROUTINE RKF45 FOR THE C... SOLUTION OF A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS (ODES). C... C... THE MODEL INITIAL CONDITIONS ARE SET IN SUBROUTINE INITAL, AND C... THE MODEL DERIVATIVES ARE PROGRAMMED IN SUBROUTINE DERV. THE C... NUMERICAL SOLUTION IS PRINTED AND PLOTTED IN SUBROUTINE PRINT. C... C... SUBROUTINES INITAL, DERV AND PRINT. THE FOLLOWING CODING IS FOR C... 1000 ORDINARY DIFFERENTIAL EQUATIONS (ODES). IF MORE ODES ARE TO C... BE INTEGRATED, ALL OF THE 1000*S SHOULD BE CHANGED TO THE REQUIRED C... NUMBER IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y(1000) 2 /F/ F(1000) C... C... THE NUMBER OF DIFFERENTIAL EQUATIONS IS IN COMMON/N/ FOR USE IN C... SUBROUTINE FCN COMMON/N/ NEQN C... C... COMMON AREA TO PROVIDE THE INPUT/OUTPUT UNIT NUMBERS TO OTHER C... SUBROUTINES COMMON/IO/ NI, NO C... C... ABSOLUTE DIMENSIONING OF THE ARRAYS REQUIRED BY RKF45 DIMENSION YV(1000), WORK(6003), IWORK(5) C... C... EXTERNAL THE DERIVATIVE ROUTINE CALLED BY RKF45 EXTERNAL FCN C... C... ARRAY FOR THE TITLE (FIRST LINE OF DATA), CHARACTERS END OF RUNS CHARACTER TITLE(20)*4, ENDRUN(3)*4 C... C... DEFINE THE CHARACTERS END OF RUNS DATA ENDRUN/'END ','OF R','UNS '/ C... C... DEFINE THE INPUT/OUTPUT UNIT NUMBERS NI=5 NO=6 C... C... DEFINE INPUT/OUTPUT FILES OPEN(NI,FILE='DATA' ,STATUS='OLD') OPEN(NO,FILE='OUTPUT',STATUS='NEW') C... C... INITIALIZE THE RUN COUNTER NORUN=0 C... C... BEGIN A RUN 1 NORUN=NORUN+1 C... C... INITIALIZE THE RUN TERMINATION VARIABLE NSTOP=0 C... C... READ THE FIRST LINE OF DATA READ(NI,1000,END=999)(TITLE(I),I=1,20) C... C... TEST FOR END OF RUNS IN THE DATA DO 2 I=1,3 IF(TITLE(I).NE.ENDRUN(I))GO TO 3 2 CONTINUE C... C... AN END OF RUNS HAS BEEN READ, SO TERMINATE EXECUTION 999 STOP C... C... READ THE SECOND LINE OF DATA 3 READ(NI,1001,END=999)T0,TF,TP C... C... READ THE THIRD LINE OF DATA READ(NI,1002,END=999)NEQN,TOL C... C... PRINT A DATA SUMMARY WRITE(NO,1003)NORUN,(TITLE(I),I=1,20),T0,TF,TP,NEQN,TOL C... C... INITIALIZE TIME T=T0 C... C... SET THE INITIAL CONDITIONS CALL INITAL C... C... SET THE INITIAL CONDITIONS FOR SUBROUTINE RKF45 DO 5 I=1,NEQN YV(I)=Y(I) 5 CONTINUE C... C... SET THE PARAMETERS FOR SUBROUTINE RKF45 C... C... INITIAL TIME TV=T0 C... C... ERROR TOLERANCE ABSERR=TOL RELERR=TOL C... C... INITIALIZATION FOR MULTI-STEP MODE (RKF45 NOT CONSTRAINED C... BY OUTPUT INTERVAL TP) IFLAG=1 C... C... INITIALIZATION FOR ONE-STEP MODE (RKF45 CONSTRAINED BY C... OUTPUT INTERVAL TP) C... IFLAG=-1 C... C... FINAL TIME TEND=T0 C... C... CALL SUBROUTINE RKF45 TO START THE SOLUTION FROM THE INITIAL C... CONDITION (IFLAG = 1 OR -1) OR COMPUTE THE SOLUTION TO THE NEXT C... PRINT POINT (IFLAG = 2 OR -2) 4 CALL RKF45(FCN,NEQN,YV,TV,TEND,RELERR,ABSERR,IFLAG,WORK,IWORK) C... C... PRINT THE INITIAL CONDITION (IFLAG = 1 OR -1) OR THE SOLUTION AT C... THE NEXT OUTPUT POINT (IFLAG = 2 OR -2) T=TV DO 6 I=1,NEQN Y(I)=YV(I) 6 CONTINUE CALL PRINT(NI,NO) C... C... TEST FOR AN ERROR CONDITION IF((IFLAG.NE.2).AND.(IFLAG.NE.-2))THEN C... C... PRINT A MESSAGE INDICATING AN ERROR CONDITION WRITE(NO,1004)IFLAG C... C... GO ON TO THE NEXT RUN GO TO 1 END IF C... C... CHECK FOR A RUN TERMINATION IF(NSTOP.NE.0)GO TO 1 C... C... CHECK FOR THE END OF THE RUN TEND=TV+TP IF(TV.LT.(TF-0.5D0*TP))GO TO 4 C... C... THE CURRENT RUN IS COMPLETE, SO GO ON TO THE NEXT RUN GO TO 1 C... C... ***************************************************************** C... C... FORMATS C... 1000 FORMAT(20A4) 1001 FORMAT(3D10.0) 1002 FORMAT(I5,20X,D10.0) 1003 FORMAT(1H1, 1 ' RUN NO. - ',I3,2X,20A4,//, 2 ' INITIAL T - ',D10.3,//, 3 ' FINAL T - ',D10.3,//, 4 ' PRINT T - ',D10.3,//, 5 ' NUMBER OF DIFFERENTIAL EQUATIONS - ',I3,//, 6 ' INTEGRATION ALGORITHM - RKF45 ',//, 7 ' MAXIMUM INTEGRATION ERROR - ',D10.3,//, 8 1H1) 1004 FORMAT(1H ,//,' IFLAG = ',I3,//, 1 ' INDICATING AN INTEGRATION ERROR, SO THE CURRENT RUN' ,/, 2 ' IS TERMINATED. PLEASE REFER TO THE DOCUMENTATION FOR' ,/, 3 ' SUBROUTINE',//,25X,'RKF45',//, 4 ' FOR AN EXPLANATION OF THESE ERROR INDICATORS' ) END SUBROUTINE FCN(TV,YV,YDOT) C... C... SUBROUTINE FCN IS AN INTERFACE ROUTINE BETWEEN SUBROUTINES RKF45 C... AND DERV C... C... ODE COMMON AREA IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y(1) 2 /F/ F(1) C... C... THE NUMBER OF DIFFERENTIAL EQUATIONS IS AVAILABLE THROUGH COMMON C... /N/ COMMON/N/ NEQN C... C... ABSOLUTE DIMENSION THE DEPENDENT, DERIVATIVE VECTORS DIMENSION YV(1000), YDOT(1000) C... C... TRANSFER THE INDEPENDENT VARIABLE, DEPENDENT VARIABLE VECTOR C... FOR USE IN SUBROUTINE DERV T=TV DO 1 I=1,NEQN Y(I)=YV(I) 1 CONTINUE C... C... EVALUATE THE DERIVATIVE VECTOR CALL DERV C... C... TRANSFER THE DERIVATIVE VECTOR FOR USE BY SUBROUTINE RKF45 DO 2 I=1,NEQN YDOT(I)=F(I) 2 CONTINUE RETURN END