      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
