*DECK PRO3P10
      SUBROUTINE INITAL
C...
C...  SUBROUTINE INITAL IS CALLED BY MAIN PROGRAM PRO1P3 TO DEFINE THE
C...  INITIAL CONDITION IN THE NUMOL SOLUTION OF EQUATIONS (3.73) TO
C...  (3.75)
C...
      COMMON/T/   TIME,  NSTOP,  NORUN
     1      /Y/  T(51)
     2      /F/ TT(51)
     3      /P/      L,      V,      A,      W,     PI,      S,    SSE,
     4            NFCN,      N,     IP
C...
C...  TYPE SELECTED VARIABLES AS REAL
      REAL L
C...
C...  LENGTH, VELOCITY
      L=1.0
      V=1.0
C...
C...  SELECT A PARTICULAR FUNCTION AND SET ITS PARAMETERS
C...
C...     (1)  NFCN = 1 - GAUSSIAN STEP FUNCTION
C...
C...     (2)  NFCN = 2 - TRUNCATED RAMP FUNCTION
C...
C...     (3)  NFCN = 3 - COSINE STEP FUNCTION
C...
C...     (4)  NFCN = 4 - COSINE PULSE FUNCTION
C...
C...  FIGURES 3.8 AND 3.9
      IF((NORUN.EQ.1).OR.(NORUN.EQ.2))NFCN=2
C...
C...  FIGURES 3.10 AND 3.11
      IF((NORUN.EQ.3).OR.(NORUN.EQ.4))NFCN=4
      GO TO(1,2,3,4),NFCN
C...
C...  GAUSSIAN STEP FUNCTION
1     A=50.
      GO TO 10
C...
C...  TRUNCATED RAMP FUNCTION
2     S=5.
      GO TO 10
C...
C...  COSINE STEP FUNCTION
3     W=5.
      PI=4.*ATAN(1.)
      GO TO 10
C...
C...  COSINE PULSE FUNCTION
4     W=5.
      PI=4.*ATAN(1.)
      GO TO 10
C...
C...  GRID SPACING
10    N=51
C...
C...  INITIAL CONDITION (3.74)
      DO 11 I=1,N
      T(I)=0.
11    CONTINUE
C...
C...  INITIALIZE INDEX FOR PLOTTED POINTS, SUM OF SQUARES OF ERRORS,
C...  DERIVATIVES
      IP=0
      SSE=0.
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
C...
C...  SUBROUTINE DERV IS CALLED BY INTEGRATOR RKF45 TO DEFINE THE
C...  ODES IN THE NUMOL SOLUTION OF EQUATIONS (3.73) TO (3.75) USING
C...  SECOND ORDER APPROXIMATIONS (DSS002 CALLED WHEN NORUN = 1, 3) OR
C...  FOURTH ORDER APPROXIMATIONS (DSS004 CALLED WHEN NORUN = 2, 4)
C...
      COMMON/T/   TIME,  NSTOP,  NORUN
     1      /Y/  T(51)
     2      /F/ TT(51)
     3      /P/      L,      V,      A,      W,     PI,      S,    SSE,
     4            NFCN,      N,     IP
C...
C...  TYPE SELECTED VARIABLES AS REAL
      REAL L, TX(51)
C...
C...  LEFT BOUNDARY CONDITION
      GO TO(1,2,3,4),NFCN
C...
C...  GAUSSIAN STEP FUNCTION
1     IF(TIME.LT.0.)T(1)=0.
      AT2=-A*TIME**2
      IF((TIME.GE.0.).AND.(ABS(AT2).GT.100.))T(1)=1.
      IF((TIME.GE.0.).AND.(ABS(AT2).LE.100.))T(1)=1.-EXP(AT2)
      GO TO 10
C...
C...  TRUNCATED RAMP FUNCTION
2     IF(TIME.LT.0.)T(1)=0.
      IF(TIME.GT.(1./S))T(1)=1.
      IF((TIME.GE.0.).AND.(TIME.LE.(1./S)))T(1)=S*TIME
      GO TO 10
C...
C...  COSINE STEP FUNCTION
3     IF(TIME.LT.0.)T(1)=0.
      IF((W*TIME).GT.(PI/2.))T(1)=1.
      IF(((W*TIME).GE.(0.0)).AND.((W*TIME).LE.(PI/2.)))
     1  T(1)=1.-COS(W*TIME)
      GO TO 10
C...
C...  COSINE PULSE FUNCTION
4     IF(TIME.LT.0.)T(1)=0.
      IF((W*TIME).GT.(PI))T(1)=0.
      IF(((W*TIME).GE.(0.0)).AND.((W*TIME).LE.(PI/2.)))
     1  T(1)=1.-COS(W*TIME)
      IF(((W*TIME).GT.(PI/2.)).AND.((W*TIME).LE.(PI)))
     1  T(1)=1.+COS(W*TIME)
      GO TO 10
10    TT(1)=0.
C...
C...  COMPUTE THE SPATIAL DERIVATIVE TX
C...
C...  FIGURES 3.8 AND 3.10
      IF((NORUN.EQ.1).OR.(NORUN.EQ.3))CALL DSS012(0.,L,N,T,TX,V)
C...
C...  FIGURES 3.9 AND 3.11
      IF((NORUN.EQ.2).OR.(NORUN.EQ.4))CALL DSS020(0.,L,N,T,TX,V)
C...
C...  EQUATION (3.73)
      DO 11 I=2,N
      TT(I)=-V*TX(I)
11    CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
C...
C...  SUBROUTINE PRINT IS CALLED BY MAIN PROGRAM PRO1P3 TO PRINT THE
C...  NUMOL SOLUTION OF EQUATIONS (3.73) TO (3.75)
C...
      COMMON/T/   TIME,  NSTOP,  NORUN
     1      /Y/  T(51)
     2      /F/ TT(51)
     3      /P/      L,      V,      A,      W,     PI,      S,    SSE,
     4            NFCN,      N,     IP
C...
C...  TYPE SELECTED VARIABLES REAL
      REAL L
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IF(IP.EQ.0)THEN
      WRITE(NO,103)
103   FORMAT('   ADVECTION EQUATION',//,
     1    3X,'   TIME      T(0,t)      T(1,t)     TE(1,t)       DIFF')
      END IF
C...
C...  COMPUTE THE EXACT SOLUTION FOR PRINTING
      TXV=TIME-L/V
      GO TO(1,2,3,4),NFCN
C...
C...  GAUSSIAN STEP FUNCTION
1     IF(TXV.LT.0.)TE=0.
      AT2=-A*TXV**2
      IF((TXV.GE.0.).AND.(ABS(AT2).GT.100.))TE=1.
      IF((TXV.GE.0.).AND.(ABS(AT2).LE.100.))TE=1.-EXP(AT2)
      GO TO 10
C...
C...  TRUNCATED RAMP FUNCTION
2     IF(TXV.LT.0.)TE=0.
      IF(TXV.GT.(1./S))TE=1.
      IF((TXV.GE.0.).AND.(TXV.LE.(1./S)))TE=S*TXV
      GO TO 10
C...
C...  COSINE STEP FUNCTION
3     IF(TXV.LT.0.)TE=0.
      IF((W*TXV).GT.(PI/2.))TE=1.
      IF(((W*TXV).GE.(0.0)).AND.((W*TXV).LE.(PI/2.)))TE=1.-COS(W*TXV)
      GO TO 10
C...
C...  COSINE PULSE FUNCTION
4     IF(TXV.LT.0.)TE=0.
      IF((W*TXV).GT.(PI))TE=0.
      IF(((W*TXV).GE.(0.0)).AND.((W*TXV).LE.(PI/2.)))TE=1.-COS(W*TXV)
      IF(((W*TXV).GT.(PI/2.)).AND.((W*TXV).LE.(PI))) TE=1.+COS(W*TXV)
      GO TO 10
C...
C...  COMPUTE THE FIGURE OF MERIT, SSE
10    DIFF=T(N)-TE
      SSE=SSE+DIFF**2
C...
C...  PRINT THE NUMOL AND EXACT SOLUTIONS
      WRITE(NO,11)TIME,T(1),T(N),TE,DIFF
11    FORMAT(F10.1,4F12.4)
      IP=IP+1
C...
C...  PRINT THE SUM OF SQUARES OF THE ERRORS AT THE END OF THE RUN
      IF(IP.LT.21)RETURN
      WRITE(NO,104)SSE
104   FORMAT(/,' SSE = ',F7.5,//)
      IP=0
      RETURN
      END
