      SUBROUTINE INITAL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NR=21)
      COMMON/T/       Z,   NSTOP,   NORUN
     +      /Y/   U(NR)
     +      /F/  UZ(NR)
     +      /S/  UR(NR), URR(NR),   V(NR),      PE,
     +               R0,      DR,     DRS,   R(NR),
     +               ZL
     +      /I/      IP
C...
C...
C...  PROBLEM PARAMETERS
C...
C...     PECLET NUMBER
         IF(NORUN.EQ.1)PE= 60.0D0
         IF(NORUN.EQ.2)PE=120.0D0
         IF(NORUN.EQ.3)PE=180.0D0
C...
C...     RADIUS
         R0=1.0D0
C...
C...     LENGTH
         ZL=30.0D0
C...
C...  RADIAL GRID, VELOCITY PROFILE
      DR=R0/DFLOAT(NR-1)
      DRS=DR**2
      DO 1 I=1,NR
         R(I)=DR*DFLOAT(I-1)
         V(I)=(1.0D0-R(I)**2)
1     CONTINUE
C...
C...  INITIAL CONDITION
      DO 2 I=1,NR
         U(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=21)
      COMMON/T/       Z,   NSTOP,   NORUN
     +      /Y/   U(NR)
     +      /F/  UZ(NR)
     +      /S/  UR(NR), URR(NR),   V(NR),      PE,
     +               R0,      DR,     DRS,   R(NR),
     +               ZL
     +      /I/      IP
C...
C...  BOUNDARY CONDITION AT R = R0 = 1
       U(NR)=1.0D0
      UZ(NR)=0.0D0
      NU=1
C...
C...  UR
      CALL DSS004(0.0D0,R0,NR,U,UR)
C...
C...  BOUNDARY CONDITION AT R = 0
      UR(1)=0.0D0
      NL=2
C...
C...  URR
      CALL DSS044(0.0D0,R0,NR,U,UR,URR,NL,NU)
C...
C...  PDE
      DO 1 I=1,NR-1
C...
C...     R = 0
         IF(I.EQ.1)THEN
         UZ(I)=1.0D0/(PE*2.0D0*V(I))*2.0D0*URR(I)
C...
C...     R NE 0
         ELSE IF(I.NE.1)THEN
         UZ(I)=1.0D0/(PE*2.0D0*V(I))*(URR(I)+(1.0D0/R(I))*UR(I))
         END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NR=21)
      COMMON/T/       Z,   NSTOP,   NORUN
     +      /Y/   U(NR)
     +      /F/  UZ(NR)
     +      /S/  UR(NR), URR(NR),   V(NR),      PE,
     +               R0,      DR,     DRS,   R(NR),
     +               ZL
     +      /I/      IP
C...
C...  MONITOR THE OUTPUT
      WRITE(*,*)NORUN,Z
C...
C...  OPEN A MATLAB FILE FOR PLOTTING
      IP=IP+1
      IF(IP.EQ.1)THEN
         OPEN(1,FILE='h5a.out')
         OPEN(2,FILE='h5b.out')
      END IF
C...
C...  WRITE THE NUMERICAL AND AND ANALYTICAL SOLUTIONS AT SELECTED Z
      IF(((IP-1)*(IP-51)*(IP-101)).EQ.0)THEN
C...
C...     NUMERICAL SOLUTION FOR 0 LE R LE R0
         NRI=4
         WRITE(NO,1)Z,(U(I),I=1,NR,NRI),
     +               (UZ(I),I=1,NR,NRI)
1        FORMAT(//,' z = ',F4.1,/,
     +    12X,'r=0',4X,'r=0.2',4X,'r=0.4',
     +     4X,'r=0.6',4X,'r=0.8',4X,'r=1.0',
     +    /,'  u(r,z)',6F8.5,
     +    /,' uz(r,z)',6F8.5,/)
C...
C...     SERIES SOLUTION AT R = 0
         CALL SERIES(PE,Z,UA)
         WRITE(NO,5)PE,Z,UA
5        FORMAT(/,' Series solution, Ua(0,z)',//,
     +         2X,'    Pe = ',F5.0  ,/,
     +         2X,'     z = ',F5.1  ,/,
     +         2X,'    Ua = ',F8.5 ,//)
      END IF
C...
C...  WRITE THE SOLUTION FOR MATLAB PLOTTING
      WRITE(1,2)Z,(U(I),I=1,NR)
2     FORMAT(22F8.4)
C...
C...  COMPUTE THE SERIES SOLUTION TO BE SUPERIMPOSED ON THE NUMERICAL
C...  SOLUTION (AT THE END OF THE THIRD SOLUTION)
      IF((IP.EQ.101).AND.(NORUN.EQ.3))THEN
         DO 3 NCASE=1,3
            IF(NCASE.EQ.1)PE= 60.0D0
            IF(NCASE.EQ.2)PE=120.0D0
            IF(NCASE.EQ.3)PE=180.0D0
         DO 4 I=1,101
            ZA=ZL*DFLOAT(I-1)/DFLOAT(101-1)
            CALL SERIES(PE,ZA,UA)
            WRITE(2,2)ZA,UA
4        CONTINUE
3        CONTINUE
      END IF
      RETURN
      END
      SUBROUTINE SERIES(PE,Z,U)
C...
C...  The Graetz problem stated by Jakob (1949) is
C...
C...     2*vavg*v(r)*ua  = D*(ua   + (1/r)*ua )                   (1)
C...                   z        rr           r
C...
C...     v(r) = 1 - r**(2)
C...
C...  where ua = (T - Ts)/(To - Ts)
C...
C...  For the numerical solution, u = (T - To)/(Ts - To).  u can be
C...  expressed in terms of ua.  Thus,
C...
C...  T = Ts + ua*(To - Ts)
C...
C...  u = (Ts + ua*(To - Ts) - To)/(Ts - To) = 1 - ua
C...
C...  The analytical (series) solution given by Jakob, eq. (22-45),
C...  is
C...
C...  ua = 1.477*exp(-3.658*(1/Pe)*z)*R0(r)
C...
C...     - 0.810*exp(22.178*(1/Pe)*z)*R1(r)
C...
C...     + 0.385*exp(-53.05*(1/Pe)*z)*R2(r)
C...
C...  where z = z/ro, r = r/ro, (the LHS z and r are dimensionaless),
C...  R0(0) = R1(0) = R2(0) = 1 (from Jakob, Table,22-1, page 455),
C...  Pe = ro*vavg/D (D = k/(rho*Cp))
C...
C...  Jakob, M., Heat Transfer, John Wiley and Sons, Inc. New York,
C...  Vol. 1, 1949
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  RATIO OF Z/PE USED IN SERIES SOLUTION
      ZPE=Z/PE
C...
C...  UA(0,Z)
      TERM1= 1.477D0*DEXP( -3.658D0*ZPE)
      TERM2=-0.810D0*DEXP(-22.178D0*ZPE)
      TERM3= 0.385D0*DEXP( -53.05D0*ZPE)
      UA=TERM1+TERM2+TERM3
C...
C...  FOR COMPARISON WITH THE NUMERICAL SOLUTION, U = 1 - UA
      U=1.0D0-UA
C...
C...  END OF CALCULATION
      RETURN
      END
