*DECK PRO1P2
      PROGRAM PRO1P2
C...
C...  MAIN PROGRAM FOR THE NUMOL INTEGRATION OF EQUATIONS (1.28) TO
C...  (1.31) USING DIFFERENTIATION IN SPACE SUBROUTINE DSS002 AND ODE
C...  INTEGRATOR RKF45
C...
C...  SET THE NUMBER OF ODES
      PARAMETER (N=51)
C...
C...  ABSOLUTE DIMENSION THE ARRAYS (WORK AND IWORK ARE REAL AND INTEGER
C...  WORK ARRAYS REQUIRED BY SUBROUTINE RKF45)
      REAL T(N), TT(N), TE(N), DIFF(N), WORK(1000)
      INTEGER IWORK(5)
C...
C...  EXTERNAL THE DERIVATIVE SUBROUTINE CALLED BY RKF45
      EXTERNAL FCN
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...  SET THE INITIAL, FINAL AND PRINT INTERVAL TIMES
      TIME=0.
      TF=0.5
      TP=0.1
C...
C...  SET THE PARAMETERS FOR SUBROUTINE RKF45
      RELERR=1.0E-04
      ABSERR=1.0E-06
      IFLAG=1
      TOUT=TP
C...
C...  CALL SUBROUTINE FCN TO SET THE INITIAL CONDITION
      CALL FCN(TIME,T,TT)
C...
C...  CALCULATE THE EXACT SOLUTION, AND THE DIFFERENCE BETWEEN THE
C...  NUMOL AND EXACT SOLUTIONS
1     DO 3 I=1,N,10
      TE(I)=EXACT(I,TIME)
      DIFF(I)=T(I)-TE(I)
3     CONTINUE
C...
C...  PRINT THE INITIAL CONDITION (IFLAG = 1) OR THE SOLUTION AT THE
C...  NEXT OUTPUT POINT (IFLAG = 2)
      WRITE(NO,2)TIME,(T(I),I=1,N,10),
     1               (TE(I),I=1,N,10),
     2             (DIFF(I),I=1,N,10)
2     FORMAT(' TIME = ',F6.2,/,15X,'  X=0',5X,'X=0.2',5X,'X=0.4',5X,
     1                             'X=0.6',5X,'X=0.8',5X,'  X=1',/,
     2                         '    T(X,T)',6F10.6,/,
     3                         '   TE(X,T)',6F10.6,/,
     4                         ' DIFF(X,T)',6F10.6,/)
C...
C...  CHECK FOR THE FINAL VALUE OF TIME
      IF(TIME.GT.(TF-0.5*TP))THEN
C...
C...     THE SOLUTION IS COMPLETE SO TERMINATE EXECUTION
         STOP
      ELSE
C...
C...     SET THE NEXT OUTPUT TIME AND CONTINUE THE INTEGRATION
         TOUT=TIME+TP
      END IF
C...
C...  CALL SUBROUTINE RKF45 TO START THE SOLUTION FROM THE INITIAL
C...  CONDITION (IFLAG = 1) OR COMPUTE THE SOLUTION TO THE NEXT OUTPUT
C...  POINT (IFLAG = 2) AT TOUT
      CALL RKF45(FCN,N,T,TIME,TOUT,RELERR,ABSERR,IFLAG,WORK,IWORK)
C...
C...  TEST FOR AN ERROR CONDITION
      IF(IFLAG.NE.2)THEN
C...
C...     PRINT A MESSAGE INDICATING AN ERROR CONDITION AND TERMINATE
C...     EXECUTION
         WRITE(NO,1000)IFLAG
         STOP
      ELSE
C...
C...     CONTINUE THE INTEGRATION
         GO TO 1
      END IF
C...
C...  *****************************************************************
C...
C...  FORMATS
C...
1000  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(TIME,T,TT)
C...
C...  SUBROUTINE FCN IS CALLED BY INTEGRATOR RKF45 TO DEFINE THE
C...  INITIAL CONDITION AND ODE DERIVATIVES IN THE NUMOL SOLUTION OF
C...  EQUATIONS (1.28) TO (1.31)
C...
C...  DIMENSION THE ARRAYS REQUIRED IN FCN
      PARAMETER (N=51)
      REAL T(N), TT(N), X(N), TX(N), TXX(N)
C...
C...  TYPE SELECTED VARIABLES AS REAL
      REAL L
C...
C...  CALCULATE PI FOR USE IN DO LOOP 1
      PI=ACOS(-1.0)
C...
C...  LENGTH
      L=1.0
C...
C...  GRID SPACING
      DX=L/FLOAT(N-1)
C...
C...  INITIAL CONDITION (1.22)
      IF(TIME.EQ.0.)THEN
         DO 1 I=1,N
         X(I)=FLOAT(I-1)*DX
         T(I)=COS(PI*X(I)/L)
1        CONTINUE
      END IF
C...
C...  DERIVATIVE TX
      CALL DSS002(0.,L,N,T,TX)
C...
C...  BOUNDARY CONDITION (1.30)
      TX(1)=0.
C...
C...  BOUNDARY CONDITION (1.31)
      TX(N)=0.
C...
C...  DERIVATIVE TXX
      CALL DSS002(0.,L,N,TX,TXX)
C...
C...  EQUATION (1.28)
      DO 2 I=1,N
      TT(I)=TXX(I)
2     CONTINUE
      RETURN
      END
      REAL FUNCTION EXACT(I,T)
C...
C...  FUNCTION EXACT COMPUTES THE EXACT SOLUTION TO EQUATIONS (1.28)
C...  TO (1.31), I.E., EQUATION (1.32), AT GRID INDEX I AND TIME T
C...
C...  NUMBER OF GRID POINTS
      PARAMETER (N=51)
C...
C...  TYPE SELECTED VARIABLES AS REAL
      REAL L
C...
C...  PI FOR USE IN THE EXACT SOLUTION
      PI=ACOS(-1.)
C...
C...  LENGTH
      L=1.
C...
C...  GRID SPACING
      DX=L/FLOAT(N-1)
C...
C...  X FOR WHICH THE EXACT SOLUTION IS TO BE COMPUTED
      X=FLOAT(I-1)*DX
C...
C...  EXACT SOLUTION AT X AND T
      EXACT=EXP((-PI**2/L**2)*T)*COS(PI*X/L)
      RETURN
      END
