      SUBROUTINE INITAL
C...
C...  H4 - DYNAMIC GRAETZ PROBLEM WITH CONSTANT WALL HEAT FLUX
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=11,NZ=31)
      COMMON/T/          T,     NSTOP,     NORUN
     +      /Y/   U(NR,NZ)
     +      /F/  UT(NR,NZ)
     +      /G/      R(NR),     Z(NZ),     V(NR)
     +      /C/         PE,        TN,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  SELECT APPROXIMATION OF SPATIAL DERIVATIVES (SEE THE COMMENTS
C...  IN SUBROUTINE DERV FOR AN EXPLANATION OF THE TWO CASES)
      NCASE=2
C...
C...  PROBLEM PARAMETERS
C...
C...     PECLET NUMBER
         PE=60.0D0
C...
C...     QW*R0/K = TN (NORMALIZATION TEMPERATURE)
         TN=1.0D0
C...
C...     OUTSIDE RADIUS
         R0=1.0D0
C...
C...     AXIAL LENGTH
         IF(NORUN.EQ.1)ZL=30.0D0
         IF(NORUN.EQ.2)ZL=90.0D0
C...
C...  RADIAL GRID, VELOCITY PROFILE
      DR=R0/DFLOAT(NR-1)
      DO 1 IR=1,NR
         R(IR)=DFLOAT(IR-1)*DR
         V(IR)=2.0D0*(1.0D0-R(IR)**2)
1     CONTINUE
      DRS=DR**2
C...
C...  AXIAL GRID
      DZ=ZL/DFLOAT(NZ-1)
      DO 2 JZ=1,NZ
         Z(JZ)=DFLOAT(JZ-1)*DZ
2     CONTINUE
      DZS=DZ**2
C...
C...  INITIAL CONDITION
      DO 3 IR=1,NR
      DO 3 JZ=1,NZ
         U(IR,JZ)=1.0D0
3     CONTINUE
C...
C...  GIVE THE TIME DERIVATIVES A VALUE WHICH CAN THEN BE USED TO
C...  TEST FOR UNDEFINED DERIVATIVES BY SETTING THEM TO PI, THEN
C...  CALLING DERV, AND CHECKING IN SUBROUTINE PRINT IF THEY HAVE
C...  ALL BEEN SET TO VALUES OTHER THAN PI.  THIS IS A SIMPLE PRO-
C...  CEDURE TO ENSURE THAT ALL NR X NZ DERIVATIVES HAVE BEEN SET
C...  IN DERV
      PI=3.1415927D0
      DO 4 IR=1,NR
      DO 4 JZ=1,NZ
         UT(IR,JZ)=PI
4     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/I/         IP,     NCASE
C...
C...  APPROXIMATIONS OF SPATIAL DERIVATIVES
C...
C...     NCASE = 1 - EXPLICIT FINITE DIFFERENCES (FDM)
         IF(NCASE.EQ.1)THEN
            CALL DERV1
C...
C...     NCASE = 2 - FINITE DIFFERENCES FROM DSS034
         ELSE
     +   IF(NCASE.EQ.2)THEN
            CALL DERV2
         END IF
      RETURN
      END
      SUBROUTINE DERV1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=11,NZ=31)
      COMMON/T/          T,     NSTOP,     NORUN
     +      /Y/   U(NR,NZ)
     +      /F/  UT(NR,NZ)
     +      /G/      R(NR),     Z(NZ),     V(NR)
     +      /C/         PE,        TN,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  Z = 0
      DO 1 IR=1,NR
          U(IR,1)=1.0D0
         UT(IR,1)=0.0D0
1     CONTINUE
C...
C...  R = R0, Z NE 0
      DO 2 JZ=2,NZ
         UNRP1=U(NR-1,JZ)+(2.0D0*DR)
         UT(NR,JZ)=-PE*0.0D0
     +            +(UNRP1-2.0D0*U(NR,JZ)+U(NR-1,JZ))/DRS
     +            +(UNRP1-U(NR-1,JZ))/(2.0D0*DR)
2     CONTINUE
C...
C...  R = 0, Z NE 0
      DO 3 JZ=2,NZ
         UT(1,JZ)=-PE*V(1)*(U(1,JZ)-U(1,JZ-1))/DZ
     +            +4.0D0*(U(2,JZ)-U(1,JZ))/DRS
3     CONTINUE
C...
C...  R NE 0, R0, Z NE 0
      DO 4 JZ=2,NZ
      DO 4 IR=2,NR-1
         UT(IR,JZ)=-PE*V(IR)*(U(IR,JZ)-U(IR,JZ-1))/DZ
     +             +(U(IR+1,JZ)-2.0D0*U(IR,JZ)+U(IR-1,JZ))/DRS
     +             +(1.0D0/R(IR))*(U(IR+1,JZ)-U(IR-1,JZ))/(2.0D0*DR)
4     CONTINUE
      RETURN
      END
      SUBROUTINE DERV2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=11,NZ=31)
      COMMON/T/          T,     NSTOP,     NORUN
     +      /Y/   U(NR,NZ)
     +      /F/  UT(NR,NZ)
     +      /G/      R(NR),     Z(NZ),     V(NR)
     +      /C/         PE,        TN,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  ARRAYS FOR UR, URR, UZ
      DIMENSION UR(NR,NZ), URR(NR,NZ), UZ(NR,NZ)
C...
C...  Z = 0
      DO 1 IR=1,NR
          U(IR,1)=1.0D0
         UT(IR,1)=0.0D0
1     CONTINUE
C...
C...  UR
      CALL DSS034(0.0D0,R0,NR,NZ,1,U,UR,0.0D0)
C...
C...  R = R0, Z NE 0
      DO 2 JZ=2,NZ
         UR(NR,JZ)=1.0D0
2     CONTINUE
C...
C...  R = 0, Z NE 0
      DO 3 JZ=2,NZ
         UR(1,JZ)=0.0D0
3     CONTINUE
C...
C...  URR
      CALL DSS034(0.0D0,R0,NR,NZ,1,UR,URR,0.0D0)
C...
C...  UZ (FOR CONVECTION)
      CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,1.0D0)
C...
C...  Z NE 0
      DO 4 JZ=2,NZ
      DO 4 IR=1,NR
C...
C...     R = 0
         IF(IR.EQ.1)THEN
            UT(1,JZ)=-PE*V(1)*UZ(1,JZ)
     +               +2.0D0*URR(1,JZ)
C...
C...     R NE 0
         ELSE
     +   IF(IR.NE.1)THEN
            UT(IR,JZ)=-PE*V(IR)*UZ(IR,JZ)
     +                +URR(IR,JZ)+(1.0D0/R(IR))*UR(IR,JZ)
         END IF
4     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=11,NZ=31)
      COMMON/T/          T,     NSTOP,     NORUN
     +      /Y/   U(NR,NZ)
     +      /F/  UT(NR,NZ)
     +      /G/      R(NR),     Z(NZ),     V(NR)
     +      /C/         PE,        TN,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  ARRAYS FOR THE ANALYTICAL SOLUTION
      DIMENSION UZ(NR), UA(NR)
C...
C...  MONITOR SOLUTION
      WRITE(*,*)NCASE,NORUN,T
C...
C...  PRINT THE INITIAL CONDITION AND INITIAL TEMPORAL DERIVATIVE
      IF(IP.EQ.0)THEN
         WRITE(NO,10)T,((U(IR,JZ),IR=1,NR),JZ=1,NZ)
10       FORMAT(' t = ',F6.2,//,'   u(r,z,t)',/,(11F7.3))
         WRITE(NO,11)(( UT(IR,JZ),IR=1,NR),JZ=1,NZ)
11       FORMAT(//,'  ut(r,z,t)',/,(11F7.1))
C...
C...  CHECK FOR UNDEFINED DERIVATIVES
      DO 12 IR=1,NR
      DO 12 JZ=1,NZ
         IF((UT(IR,JZ).GT.3.14D0).AND.(UT(IR,JZ).LT.3.15D0))THEN
         WRITE (NO,13)IR,JZ,U(IR,JZ),UT(IR,JZ)
13       FORMAT(' IR = ',I3,'   JZ = ',I3,
     +   '  U(IR,JZ) = ',F7.3,'   UT(IR,JZ) = ',F7.1)
         END IF
12    CONTINUE
      END IF
C...
C...  PRINT THE NUMERICAL SOLUTION
      WRITE(NO,1)T
1     FORMAT(//,' t = ',F6.2,//,29X,'u(r,z,t)',//,
     + 3X,'z',4X,'r=0.0',4X,'r=0.2',4X,'r=0.4',
     +        4X,'r=0.6',4X,'r=0.8',4X,'r=1.0',/)
      DO 3 JZ=1,NZ
         NRI=2
         WRITE(NO,2)Z(JZ),(U(IR,JZ),IR=1,NR,NRI)
2        FORMAT(F4.1,6F9.4)
3     CONTINUE
C...
C...  AT THE END OF THE RUN, CALCULATE AND PRINT THE ANALYTICAL,
C...  STEADY STATE TEMPERATURE
      IP=IP+1
      IF(IP.EQ.11)THEN
         DO 14 IR=1,NR
         UZ(IR)=(U(IR,NZ)-U(IR,NZ-1))/DZ
         UA(IR)=U(NR,NZ)+2.0D0*PE*TN*UZ(IR)*
     +   ((1.0D0/ 4.0D0)*(R(IR)**2)
     +   -(1.0D0/16.0D0)*(R(IR)**4)-3.0D0/16.0D0)
14       CONTINUE
         WRITE(NO,15)(UZ(IR),IR=1,NR,NRI),(UA(IR),IR=1,NR,NRI)
15       FORMAT(/,10X,'Axial gradient; analytical, steady state,',
     +          /,16X,'thermally developed solution',//,
     +          '  uz',6F9.4,/,'  ua',6F9.4,/)
      END IF
C...
C...  STORE THE SOLUTION AND PLOT THE SOLUTION AT THE END OF THE
C...  SECOND RUN
C...
C...     PLOT RADIAL PROFILES
         CALL PLOTDR
C...
C...     PLOT CENTERLINE AXIAL PROFILES
         CALL PLOTDZ
      RETURN
      END
      SUBROUTINE PLOTDR
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=11,NZ=31)
      COMMON/T/          T,     NSTOP,     NORUN
     +      /Y/   U(NR,NZ)
     +      /F/  UT(NR,NZ)
     +      /G/      R(NR),     Z(NZ),     V(NR)
     +      /C/         PE,        TN,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  OPEN A FILE FOR MATLAB PLOTTING AT THE BEGINNING OF THE
C...  FIRST RUN
      IF((IP.EQ.1).AND.(NORUN.EQ.1))THEN
         OPEN(1,FILE='h4r.out')
      END IF
C...
C...  WRITE THE TEMPERATURES, U(R,ZL,T), FOR SUBSEQUENT MATLAB
C...  PLOTTING
      DO 2 IR=1,NR
         WRITE(1,3)R(IR),U(IR,NZ)
3        FORMAT(F10.3,F10.5)
2     CONTINUE
      RETURN
      END
      SUBROUTINE PLOTDZ
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NR=11,NZ=31)
      COMMON/T/          T,     NSTOP,     NORUN
     +      /Y/   U(NR,NZ)
     +      /F/  UT(NR,NZ)
     +      /G/      R(NR),     Z(NZ),     V(NR)
     +      /C/         PE,        TN,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  OPEN A FILE FOR MATLAB PLOTTING AT THE BEGINNING OF THE
C...  FIRST RUN
      IF((IP.EQ.1).AND.(NORUN.EQ.1))THEN
         OPEN(2,FILE='h4z.out')
      END IF
C...
C...  WRITE THE CENTERLINE TEMPERATURE, U(0,Z,T), FOR SUBSEQUENT
C...  MATLAB PLOTTING
      DO 2 JZ=1,NZ
         WRITE(2,3)Z(JZ),U(1,JZ)
3        FORMAT(F10.3,F10.5)
2     CONTINUE
      RETURN
      END
