      SUBROUTINE INITAL
C...
C...  M2 - UNSTEADY POISEUILLE FLOW IN A CIRCULAR TUBE
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NR=31)
      COMMON/T/      TAU,    NSTOP,    NORUN
     +      /Y/  PHI(NR)
     +      /F/ PHIT(NR)
     +      /C/       DR,      DRS,    R(NR)
     +      /I/       NG,       IP
C...
C...  SET NUMBER OF RADIAL GRID POINTS
      IF(NORUN.EQ.1)NG=11
      IF(NORUN.EQ.2)NG=21
      IF(NORUN.EQ.3)NG=31
C...
C...  SPATIAL GRID
      DR=1.0D0/DFLOAT(NG-1)
      DRS=DR**2
      DO 1 I=1,NG
         R(I)=DFLOAT(I-1)*DR
1     CONTINUE
C...
C...  INITIAL CONDITION
      DO 2 I=1,NG
         PHI(I)=0.0D0
2     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NR=31)
      COMMON/T/      TAU,    NSTOP,    NORUN
     +      /Y/  PHI(NR)
     +      /F/ PHIT(NR)
     +      /C/       DR,      DRS,    R(NR)
     +      /I/       NG,       IP
C...
C...  STEP THROUGH NG GRID POINTS
      DO 1 I=1,NG
C...
C...     R = 0
         IF(I.EQ.1)THEN
            PHIT(1)=4.0D0+4.0D0*(PHI(2)-PHI(1))/DRS
C...
C...     R EQ 1
         ELSE IF(I.EQ.NG)THEN
             PHI(NG)=0.0D0
            PHIT(NG)=0.0D0
C...
C...     R NE 0 OR 1
         ELSE
            PHIT(I)=4.0D0+(PHI(I+1)-2.0D0*PHI(I)+PHI(I-1))/DRS
     +             +(1.0D0/R(I))*(PHI(I+1)-PHI(I-1))/(2.0D0*DR)
         END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NR=31)
      COMMON/T/      TAU,    NSTOP,    NORUN
     +      /Y/  PHI(NR)
     +      /F/ PHIT(NR)
     +      /C/       DR,      DRS,    R(NR)
     +      /I/       NG,       IP
C...
C...  MONITOR SOLUTION
      WRITE(*,*)TAU,PHI(1)
C...
C...  OPEN A FILE FOR MATLAB PLOTTING OF THE SOLUTION FOR THE 31
C...  POINT GRID (NORUN = 31)
      IF((TAU.LT.0.001D0).AND.(NORUN.EQ.3))OPEN(1,FILE='m2.out')
C...
C...  WRITE SOLUTION
      IF(NORUN.EQ.1)NGI=1
      IF(NORUN.EQ.2)NGI=2
      IF(NORUN.EQ.3)NGI=3
      WRITE(NO,1)TAU,(PHI(I),I=1,NG,NGI)
1     FORMAT(/,' TAU = ',F6.2,/,'   R',
     + 3X,'  0',3X,'0.1',3X,'0.2',3X,'0.3',3X,'0.4',3X,'0.5',
     + 3X,'0.6',3X,'0.7',3X,'0.8',3X,'0.9',3X,'1.0',/,' PHI',
     + 11F6.3,/)
C...
C...  WRITE CENTERLINE VELOCITY TO MORE FIGURES FOR COMPARISON
C...  OF THE SOLUTIONS
      WRITE(NO,2)PHI(1)
2     FORMAT(' PHI(0) = ',F8.5,//)
C...
C...  WRITE THE SOLUTION FOR MATLAB PLOTTING
      IP=IP+1
      IF((NORUN.EQ.3).AND.
     +   ((IP.EQ. 1).OR.(IP.EQ. 2).OR.(IP.EQ. 3).OR.
     +    (IP.EQ. 4).OR.(IP.EQ. 5).OR.(IP.EQ. 7).OR.
     +    (IP.EQ. 9).OR.(IP.EQ.11).OR.(IP.EQ.21)))THEN
      DO 3 I=1,NG
         WRITE(1,4)R(I),PHI(I)
4        FORMAT(F10.3,F10.5)
3        CONTINUE
      END IF
      RETURN
      END
