PROGRAM M2EX
C...
C...  PROGRAM M2EX COMPUTES AN EXPLICIT FINITE DIFFERENCE SOLUTION TO
C...  THE UNSTEADY POISEUILLE FLOW PDE IN A CIRCULAR TUBE.  THIS IS
C...  EQUIVALENT TO AN EXPLICIT EULER INTEGRATION
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
C...
C...  COMMON AREA FOR INPUT/OUTPUT UNIT NUMBERS
      COMMON/IO/ IO, NO
C...
C...  OPEN OUTPUT FILE
      NO=1
      OPEN(NO,FILE='OUTPUT',STATUS='UNKNOWN')
C...
C...  SET INITIAL CONDITIONS
      TAU=0.0D0
      CALL INITAL
C...
C...  OUTPUT INTERVAL FOR SOLUTION
      TAUP=0.05D0
C...
C...  UPDATE THE DERIVATIVE VECTOR WITH BY A CALL TO DERV AND PRINT
C...  SOLUTION (INITIAL CONDITION THE FIRST CALL TO PRINT)
1     CALL DERV
      CALL PRINT(NI,NO)
C...
C...  CHECK FOR END OF SOLUTION
      IF(TAU.GT.0.99D0)STOP
C...
C...  TAKE NSTEPS EULER STEPS, EACH OF LENGTH DTAU
      NSTEPS=200
      DTAU=TAUP/DFLOAT(NSTEPS)
      DO 2 NS=1,NSTEPS
C...
C...     DERIVATIVE VECTOR
         CALL DERV
C...
C...     EULER STEP
         DO 3 I=1,NG
         PHI(I)=PHI(I)+PHIT(I)*DTAU
3        CONTINUE
         TAU=TAU+DTAU
C...
C...     NEXT EULER STEP
2     CONTINUE
C...
C...  NSTEPS COMPLETED.  PRINT SOLUTION
      GO TO 1
      END
      SUBROUTINE INITAL
      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
C...
C...  SET NUMBER OF RADIAL GRIDS
      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
      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
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
C...
C...  MONITOR SOLUTION
      WRITE(*,*)TAU
C...
C...  WRITE SOLUTION
      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,//)
      RETURN
      END
      DOUBLE PRECISION FUNCTION DFLOAT(I)
C...
C...  FUNCTION DFLOAT CONVERTS A SINGLE PRECISION INTEGER INTO A DOUBLE
C...  PRECISION FLOATING POINT NUMBER
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DFLOAT=DBLE(FLOAT(I))
      RETURN
      END