PROGRAM GRAETZ
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...
C...  The output from this program is
C...
C     Series solution =     60.000    30.000    0.7628
C               Terms =     0.2372    0.0000    0.0000
C
C     Series solution =     60.000    90.000    0.9939
C               Terms =     0.0061    0.0000    0.0000
C
C     Series solution =    120.000    30.000    0.4113
C               Terms =     0.5919   -0.0032    0.0000
C
C     Series solution =    120.000    90.000    0.9050
C               Terms =     0.0950    0.0000    0.0000
C
C     Series solution =    180.000    30.000    0.2172
C               Terms =     0.8028   -0.0201    0.0001
C
C     Series solution =    180.000    90.000    0.7628
C               Terms =     0.2372    0.0000    0.0000
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  OUTPUT FILE
      OPEN(1,FILE='output')
C...
C...  PECLET NUMBERS
      DO 3 NPE=1,3
C...
C...    PE
        IF(NPE.EQ.1)PE= 60.0D0
        IF(NPE.EQ.2)PE=120.0D0
        IF(NPE.EQ.3)PE=180.0D0
C...
C...  TWO LENGTHS
      DO 2 NLEN=1,2
C...
C...     Z
         IF(NLEN.EQ.1)Z=30.D0
         IF(NLEN.EQ.2)Z=90.D0
         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...  WRITE SOLUTION
      WRITE(1,1)PE,Z,U,TERM1,TERM2,TERM3
1     FORMAT(' Series solution = ',2F10.3,F10.4,/,
     +   10X,' Terms = ',3F10.4,/)
C...
C...  NEXT LENGTH
2     CONTINUE
C...
C...  NEXT PE
3     CONTINUE
C...
C...  END OF CALCULATION
      STOP
      END