PROGRAM MAINODE C... C... PROGRAM MAINODE CALLS: (1) SUBROUTINE INITAL TO DEFINE THE ODE C... INITAL CONDITIONS, (2) SUBROUTINE RKF45 TO INTEGRATE THE ODES, C... AND (3) SUBROUTINE PRINT TO PRINT THE SOLUTION. C... C... THE FOLLOWING CODING IS FOR 450 ODES. IF MORE ODES ARE TO BE INTE- C... GRATED, ALL OF THE 450'S SHOULD BE CHANGED TO THE REQUIRED NUMBER IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y(450) 2 /F/ F(450) 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(450), WORK(3000), 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... OPEN INPUT AND 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,ERROR C... C... PRINT A DATA SUMMARY WRITE(NO,1003)NORUN,(TITLE(I),I=1,20), 1 T0,TF,TP, 2 NEQN,ERROR C... C... INITIALIZE TIME T=T0 C... C... SET THE INITIAL CONDITIONS CALL INITAL C... C... SET THE INITIAL DERIVATIVES (FOR POSSIBLE PRINTING) CALL DERV C... C... PRINT THE INITIAL CONDITIONS CALL PRINT(NI,NO) C... C... SET THE INITIAL CONDITIONS FOR SUBROUTINE RKF45 TV=T0 DO 5 I=1,NEQN YV(I)=Y(I) 5 CONTINUE C... C... SET THE PARAMETERS FOR SUBROUTINE RKF45 RELERR=ERROR ABSERR=ERROR IFLAG=1 TOUT=T0+TP C... C... CALL SUBROUTINE RKF45 TO START THE SOLUTION FROM THE INITIAL C... CONDITION (IFLAG = 1) OR COMPUTE THE SOLUTION TO THE NEXT PRINT C... POINT (IFLAG = 2) 4 CALL RKF45(FCN,NEQN,YV,TV,TOUT,RELERR,ABSERR,IFLAG,WORK,IWORK) C... C... PRINT THE SOLUTION AT THE NEXT OUTPUT POINT T=TV DO 6 I=1,NEQN Y(I)=YV(I) 6 CONTINUE CALL DERV CALL PRINT(NI,NO) C... C... TEST FOR AN ERROR CONDITION IF(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 TOUT=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 ' MAXIMUM INTEGRATION ERROR - ',D10.3,//, 7 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... NOTE THAT THE SIZE OF ARRAYS Y AND F IN THE FOLLOWING COMMON AREA C... IS ACTUALLY SET BY THE CORRESPONDING COMMON STATEMENT IN MAIN C... PROGRAM MAINODE 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 VARIABLE, DERIVATIVE VECTORS DIMENSION YV(450), YDOT(450) 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 DOUBLE PRECISION FUNCTION DFLOAT(I) DFLOAT=DBLE(FLOAT(I)) RETURN END