      SUBROUTINE INITAL
C...
C...  H5 - DYNAMIC GRAETZ PROBLEMS WITH CONSTANT WALL TEMPERATURE
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,        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 EIGHT CASES)
      NCASE=8
C...
C...  PROBLEM PARAMETERS
      PE=60.0D0
      R0= 1.0D0
      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)=0.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) WITH UZZ,
C...                 UZ(R,ZL,T) = 0
         IF(NCASE.EQ.1)THEN
            CALL DERV1
C...
C...     NCASE = 2 - EXPLICIT FINITE DIFFERENCES WITH UZZ,
C...                 UT(R,ZL,T) = -PE*V(R)*UZ(R,ZL,T)
         ELSE
     +   IF(NCASE.EQ.2)THEN
            CALL DERV2
C...
C...     NCASE = 3 - EXPLICIT FINITE DIFFERENCES WITH UZZ,
C...                 UT(R,ZL,T) = -PE*V(R)*UZ(R,ZL,T) +
C...                 URR(R,ZL,T) + (1/R)*UR(R,ZL,T)
         ELSE
     +   IF(NCASE.EQ.3)THEN
            CALL DERV3
C...
C...     NCASE = 4 - EXPLICIT FINITE DIFFERENCES WITHOUT UZZ
         ELSE
     +   IF(NCASE.EQ.4)THEN
            CALL DERV4
C...
C...     NCASE = 5 - FINITE DIFFERENCES FROM DSS034 WITH UZZ,
C...                 UZ(R,ZL,T) = 0
         ELSE
     +   IF(NCASE.EQ.5)THEN
            CALL DERV5
C...
C...     NCASE = 6 - FINITE DIFFERENCES FROM DSS034 WITH UZZ,
C...                 UT(R,ZL,T) = -PE*V(R)*UZ(R,ZL,T)
         ELSE
     +   IF(NCASE.EQ.6)THEN
            CALL DERV6
C...
C...     NCASE = 7 - FINITE DIFFERENCES FROM DSS034 WITH UZZ,
C...                 UT(R,ZL,T) = -PE*V(R)*UZ(R,ZL,T) +
C...                 URR(R,ZL,T) + (1/R)*UR(R,ZL,T)
         ELSE
     +   IF(NCASE.EQ.7)THEN
            CALL DERV7
C...
C...     NCASE = 8 - FINITE DIFFERENCES FROM DSS034 WITHOUT UZZ
         ELSE
     +   IF(NCASE.EQ.8)THEN
            CALL DERV8
         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,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  Z = 0
      DO 1 IR=1,NR
          U(IR,1)=0.0D0
         UT(IR,1)=0.0D0
1     CONTINUE
C...
C...  R = R0, Z NE 0
      DO 2 JZ=2,NZ
          U(NR,JZ)=1.0D0
         UT(NR,JZ)=0.0D0
2     CONTINUE
C...
C...  R = 0, Z NE 0, ZL
      DO 3 JZ=2,NZ-1
         UT(1,JZ)=-PE*V(1)*(U(1,JZ)-U(1,JZ-1))/DZ
     +            +4.0D0*(U(2,JZ)-U(1,JZ))/DRS
     +            +(U(1,JZ+1)-2.0D0*U(1,JZ)+U(1,JZ-1))/DZS
3     CONTINUE
C...
C...  R = 0, Z = ZL
         UT(1,NZ)=-PE*V(1)*0.0D0
     +            +4.0D0*(U(2,NZ)  -U(1,NZ))/DRS
     +            +2.0D0*(U(1,NZ-1)-U(1,NZ))/DZS
C...
C...  Z = ZL, R NE 0, R0
      DO 4 IR=2,NR-1
         UT(IR,NZ)=-PE*V(IR)*0.0D0
     +             +(U(IR+1,NZ)-2.0D0*U(IR,NZ)+U(IR-1,NZ))/DRS
     +             +(1.0D0/R(IR))*(U(IR+1,NZ)-U(IR-1,NZ))/(2.0D0*DR)
     +             +2.0D0*(U(IR,NZ-1)-U(IR,NZ))/DZS
4     CONTINUE
C...
C...  R NE 0, R0, Z NE 0, ZL
      DO 5 JZ=2,NZ-1
      DO 5 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)
     +             +(U(IR,JZ+1)-2.0D0*U(IR,JZ)+U(IR,JZ-1))/DZS
5     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,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  Z = 0
      DO 1 IR=1,NR
          U(IR,1)=0.0D0
         UT(IR,1)=0.0D0
1     CONTINUE
C...
C...  R = R0, Z NE 0
      DO 2 JZ=2,NZ
          U(NR,JZ)=1.0D0
         UT(NR,JZ)=0.0D0
2     CONTINUE
C...
C...  R = 0, Z NE 0, ZL
      DO 3 JZ=2,NZ-1
         UT(1,JZ)=-PE*V(1)*(U(1,JZ)-U(1,JZ-1))/DZ
     +            +4.0D0*(U(2,JZ)-U(1,JZ))/DRS
     +            +(U(1,JZ+1)-2.0D0*U(1,JZ)+U(1,JZ-1))/DZS
3     CONTINUE
C...
C...  Z = ZL, R NE R0
      DO 4 IR=1,NR-1
         UT(IR,NZ)=-PE*V(IR)*(U(IR,NZ)-U(IR,NZ-1))/DZ
4     CONTINUE
C...
C...  R NE 0, R0, Z NE 0, ZL
      DO 5 JZ=2,NZ-1
      DO 5 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)
     +             +(U(IR,JZ+1)-2.0D0*U(IR,JZ)+U(IR,JZ-1))/DZS
5     CONTINUE
      RETURN
      END
      SUBROUTINE DERV3
      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,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  Z = 0
      DO 1 IR=1,NR
          U(IR,1)=0.0D0
         UT(IR,1)=0.0D0
1     CONTINUE
C...
C...  R = R0, Z NE 0
      DO 2 JZ=2,NZ
          U(NR,JZ)=1.0D0
         UT(NR,JZ)=0.0D0
2     CONTINUE
C...
C...  R = 0, Z NE 0, ZL
      DO 3 JZ=2,NZ-1
         UT(1,JZ)=-PE*V(1)*(U(1,JZ)-U(1,JZ-1))/DZ
     +            +4.0D0*(U(2,JZ)-U(1,JZ))/DRS
     +            +(U(1,JZ+1)-2.0D0*U(1,JZ)+U(1,JZ-1))/DZS
3     CONTINUE
C...
C...  R = 0, Z = ZL
         UT(1,NZ)=-PE*V(1)*(U(1,NZ)-U(1,NZ-1))/DZ
     +            +4.0D0*(U(2,NZ)  -U(1,NZ))/DRS
C...
C...  Z = ZL, R NE 0, R0
      DO 4 IR=2,NR-1
         UT(IR,NZ)=-PE*V(IR)*(U(IR,NZ)-U(IR,NZ-1))/DZ
     +             +(U(IR+1,NZ)-2.0D0*U(IR,NZ)+U(IR-1,NZ))/DRS
     +             +(1.0D0/R(IR))*(U(IR+1,NZ)-U(IR-1,NZ))/(2.0D0*DR)
4     CONTINUE
C...
C...  R NE 0, R0, Z NE 0, ZL
      DO 5 JZ=2,NZ-1
      DO 5 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)
     +             +(U(IR,JZ+1)-2.0D0*U(IR,JZ)+U(IR,JZ-1))/DZS
5     CONTINUE
      RETURN
      END
      SUBROUTINE DERV4
      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,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  Z = 0
      DO 1 IR=1,NR
          U(IR,1)=0.0D0
         UT(IR,1)=0.0D0
1     CONTINUE
C...
C...  R = R0, Z NE 0
      DO 2 JZ=2,NZ
          U(NR,JZ)=1.0D0
         UT(NR,JZ)=0.0D0
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 5 JZ=2,NZ
      DO 5 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)
5     CONTINUE
      RETURN
      END
      SUBROUTINE DERV5
      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,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  ARRAYS FOR UR, URR, UZ, UZZ
      DIMENSION UR(NR,NZ), URR(NR,NZ), UZ(NR,NZ), UZZ(NR,NZ)
C...
C...  Z = 0
      DO 1 IR=1,NR
          U(IR,1)=0.0D0
         UT(IR,1)=0.0D0
1     CONTINUE
C...
C...  R = R0, Z NE 0
      DO 2 JZ=2,NZ
          U(NR,JZ)=1.0D0
         UT(NR,JZ)=0.0D0
2     CONTINUE
C...
C...  UR
      CALL DSS034(0.0D0,R0,NR,NZ,1,U,UR,0.0D0)
C...
C...  R = 0
      DO 3 JZ=1,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 CONDUCTION)
      CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,0.0D0)
C...
C...  Z = ZL
      DO 4 IR=1,NR
         UZ(IR,NZ)=0.0D0
4     CONTINUE
C...
C...  UZZ (FOR CONDUCTION)
      CALL DSS034(0.0D0,ZL,NR,NZ,2,UZ,UZZ,0.0D0)
C...
C...  UZ (FOR CONVECTION)
      CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,1.0D0)
C...
C...  Z = ZL
      DO 5 IR=1,NR
         UZ(IR,NZ)=0.0D0
5     CONTINUE
C...
C...  R NE R0, Z NE 0
      DO 6 JZ=2,NZ
      DO 6 IR=1,NR-1
C...
C...     R = 0
         IF(IR.EQ.1)THEN
            UT(1,JZ)=-PE*V(1)*UZ(1,JZ)
     +               +2.0D0*URR(1,JZ)
     +               +UZZ(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)
     +                +UZZ(IR,JZ)
         END IF
6     CONTINUE
      RETURN
      END
      SUBROUTINE DERV6
      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,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  ARRAYS FOR UR, URR, UZ, UZZ
      DIMENSION UR(NR,NZ), URR(NR,NZ), UZ(NR,NZ), UZZ(NR,NZ)
C...
C...  Z = 0
      DO 1 IR=1,NR
          U(IR,1)=0.0D0
         UT(IR,1)=0.0D0
1     CONTINUE
C...
C...  R = R0, Z NE 0
      DO 2 JZ=2,NZ
          U(NR,JZ)=1.0D0
         UT(NR,JZ)=0.0D0
2     CONTINUE
C...
C...  UR
      CALL DSS034(0.0D0,R0,NR,NZ,1,U,UR,0.0D0)
C...
C...  R = 0
      DO 3 JZ=1,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 CONDUCTION)
      CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,0.0D0)
C...
C...  UZZ (FOR CONDUCTION)
      CALL DSS034(0.0D0,ZL,NR,NZ,2,UZ,UZZ,0.0D0)
C...
C...  UZ (FOR CONVECTION)
      CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,1.0D0)
C...
C...  R NE R0, Z NE 0, ZL
      DO 4 JZ=2,NZ-1
      DO 4 IR=1,NR-1
C...
C...     R = 0
         IF(IR.EQ.1)THEN
            UT(1,JZ)=-PE*V(1)*UZ(1,JZ)
     +               +2.0D0*URR(1,JZ)
     +               +UZZ(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)
     +                +UZZ(IR,JZ)
         END IF
4     CONTINUE
C...
C...  R NE R0, Z = ZL
      DO 5 IR=1,NR-1
         UT(IR,NZ)=-PE*V(IR)*UZ(IR,NZ)
5     CONTINUE
      RETURN
      END
      SUBROUTINE DERV7
      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,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  ARRAYS FOR UR, URR, UZ, UZZ
      DIMENSION UR(NR,NZ), URR(NR,NZ), UZ(NR,NZ), UZZ(NR,NZ)
C...
C...  Z = 0
      DO 1 IR=1,NR
          U(IR,1)=0.0D0
         UT(IR,1)=0.0D0
1     CONTINUE
C...
C...  R = R0, Z NE 0
      DO 2 JZ=2,NZ
          U(NR,JZ)=1.0D0
         UT(NR,JZ)=0.0D0
2     CONTINUE
C...
C...  UR
      CALL DSS034(0.0D0,R0,NR,NZ,1,U,UR,0.0D0)
C...
C...  R = 0
      DO 3 JZ=1,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 CONDUCTION)
      CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,0.0D0)
C...
C...  UZZ (FOR CONDUCTION)
      CALL DSS034(0.0D0,ZL,NR,NZ,2,UZ,UZZ,0.0D0)
C...
C...  UZ (FOR CONVECTION)
      CALL DSS034(0.0D0,ZL,NR,NZ,2,U,UZ,1.0D0)
C...
C...  R NE R0, Z NE 0, ZL
      DO 4 JZ=2,NZ-1
      DO 4 IR=1,NR-1
C...
C...     R = 0
         IF(IR.EQ.1)THEN
            UT(1,JZ)=-PE*V(1)*UZ(1,JZ)
     +               +2.0D0*URR(1,JZ)
     +               +UZZ(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)
     +                +UZZ(IR,JZ)
         END IF
4     CONTINUE
C...
C...  R NE R0, Z = ZL
      DO 5 IR=1,NR-1
C...
C...     R = 0
         IF(IR.EQ.1)THEN
            UT(1,NZ)=-PE*V(1)*UZ(1,NZ)
     +               +2.0D0*URR(1,NZ)
C...
C...     R NE 0
         ELSE
     +   IF(IR.NE.1)THEN
            UT(IR,NZ)=-PE*V(IR)*UZ(IR,NZ)
     +                +URR(IR,NZ)+(1.0D0/R(IR))*UR(IR,NZ)
         END IF
5     CONTINUE
      RETURN
      END
      SUBROUTINE DERV8
      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,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
C...
C...  ARRAYS FOR UR, URR, UZ, UZZ
      DIMENSION UR(NR,NZ), URR(NR,NZ), UZ(NR,NZ), UZZ(NR,NZ)
C...
C...  Z = 0
      DO 1 IR=1,NR
          U(IR,1)=0.0D0
         UT(IR,1)=0.0D0
1     CONTINUE
C...
C...  R = R0, Z NE 0
      DO 2 JZ=2,NZ
          U(NR,JZ)=1.0D0
         UT(NR,JZ)=0.0D0
2     CONTINUE
C...
C...  UR
      CALL DSS034(0.0D0,R0,NR,NZ,1,U,UR,0.0D0)
C...
C...  R = 0
      DO 3 JZ=1,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...  R NE R0
      DO 4 JZ=2,NZ
      DO 4 IR=1,NR-1
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,        R0,        ZL,
     +                  DR,       DRS,        DZ,       DZS
     +      /I/         IP,     NCASE
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)',/,(11F6.3))
         WRITE(NO,11)(( UT(IR,JZ),IR=1,NR),JZ=1,NZ)
11       FORMAT(//,'  ut(r,z,t)',/,(11F6.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 EVERY FIFTH CALL TO PRINT
      IP=IP+1
      IF(((IP-1)/5*5).EQ.(IP-1))THEN
         WRITE(NO,1)T
1        FORMAT(//,' t = ',F6.2,//,5X,'u(r,z,t)',//,
     +    4X,'z',4X,'r = 0.0',3X,'r = 0.2',3X,'r = 0.4',
     +           3X,'r = 0.6',3X,'r = 0.8',3X,'r = 1.0')
         DO 3 JZ=1,NZ
            NRI=2
            WRITE(NO,2)Z(JZ),(U(IR,JZ),IR=1,NR,NRI)
2           FORMAT(F5.1,6F10.4)
3        CONTINUE
      END IF
C...
C...  OPEN FILE FOR MATLAB PLOTTING
      IF((NORUN.EQ.1).AND.(T.LT.0.001D0))THEN
         OPEN(1,FILE='h5a.out')
         OPEN(2,FILE='h5b.out')
      END IF
C...
C...  WRITE NUMERICAL SOLUTION AT R = 0 FOR MATLAB PLOTTING
      DO 4 JZ=1,NZ
         IF(NORUN.EQ.1)WRITE(1,5)Z(JZ),U(1,JZ)
         IF(NORUN.EQ.2)WRITE(2,5)Z(JZ),U(1,JZ)
5        FORMAT(2F10.4)
4     CONTINUE
      RETURN
      END
