SUBROUTINE INITAL
C...
C...  H6 - TUBULAR HEAT EXCHANGER
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=101)
      COMMON/T/     TIME,     NSTOP,     NORUN
     1      /Y/    T(NZ)
     2      /F/   TT(NZ)
     3      /S/   TZ(NZ)
     4      /P/        V,        ZL,       TAU,
     5               TEC,       TAC,       T0C
     6      /I/       IP,     NCASE
C...
C...  BOUNDARY CONDITION (AT Z = 0)
C...
C...     NCASE = 1 - CONSTANT
C...
C...     NCASE = 2 - VARIABLE
C...
      NCASE=1
C...
C...  SELECT ONE RUN INSTEAD OF THREE (NORUN MAY NOT INITIALLY BE 1)
      NORUN=1
C...
C...  PARAMETERS
C...
C...     LENGTH
         ZL=100.0D0
C...
C...     VELOCITY
         V=10.0D0
C...
C...     DENSITY
         RHO=1.0D0
C...
C...     SPECIFIC HEAT
         CP=1.0D0
C...
C...     TUBE DIAMETER
         D=2.0D0
C...
C...     HEAT TRANSFER COEFFICIENT
         U=0.1D0
C...
C...  TAU
      TAU=(RHO*CP*D)/(4.0D0*U)
C...
C...  ENTERING TEMPERATURE, ANNULUS TEMPERATURE
      IF(NORUN.EQ.1)THEN
         TEC=100.0D0
         TAC=0.0D0
      ELSE
     +IF(NORUN.EQ.2)THEN
         TEC=0.0D0
         TAC=100.0D0
      ELSE
     +IF(NORUN.EQ.3)THEN
         TEC=100.0D0
         TAC=100.0D0
      END IF
C...
C...  INITIAL CONDITION
      T0C=0.0D0
      DO 1 I=1,NZ
         T(I)=T0C
1     CONTINUE
C...
C...  INITIALIZE COUNTER FOR THE OUTPUT
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=101)
      COMMON/T/     TIME,     NSTOP,     NORUN
     1      /Y/    T(NZ)
     2      /F/   TT(NZ)
     3      /S/   TZ(NZ)
     4      /P/        V,        ZL,       TAU,
     5               TEC,       TAC,       T0C
     6      /I/       IP,     NCASE
C...
C...  BOUNDARY CONDITION
      IF(NCASE.EQ.1)T(1)=TEC
      IF(NCASE.EQ.2)T(1)=TEC*FBC(TIME)
      TT(1)=0.0D0
C...
C...  DERIVATIVE TZ
C...
C...     FIVE POINT BIASED UPWIND
         CALL DSS020(0.0D0,ZL,NZ,T,TZ,V)
C...
C...     VARIABLE ORDER UPWIND
C...     CALL DSS024(0.0D0,ZL,NZ,T,TZ,V)
C...
C...  PDE
      DO 1 I=2,NZ
         TT(I)=-V*TZ(I)-(1.0D0/TAU)*(T(I)-TAC)
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=101)
      COMMON/T/     TIME,     NSTOP,     NORUN
     1      /Y/    T(NZ)
     2      /F/   TT(NZ)
     3      /S/   TZ(NZ)
     4      /P/        V,        ZL,       TAU,
     5               TEC,       TAC,       T0C
     6      /I/       IP,     NCASE
C...
C...  MONITOR OUTPUT (ON UNIT *)
      IF(TIME.LT.0.01D0)WRITE(*,4)NCASE,TEC,TAC,T0C,TAU
4     FORMAT(//,' NCASE = ',I3,//,
     +         ' TEC = ',F7.2,'   TAC = ',F7.2,'   T0C = ',F7.2,
     +       '   TAU = ',F7.2,/)
      WRITE(*,3)NORUN,TIME,T(NZ),EXACT(NCASE)
3     FORMAT(' NORUN = ',I2,'    TIME = ',F6.2,'   T(NZ) = ',F7.2,
     +       '    TA = ',F7.2)
C...
C...  PRINT A HEADING FOR THE NUMERICAL AND EXACT SOLUTIONS
      IP=IP+1
      IF(IP.EQ.1)THEN
         WRITE(NO,1)
1        FORMAT(' TUBULAR HEAT EXCHANGER',//,
     1          9X,'t',5X,'T(zl,t)',4X,'Ta(zl,t)',8X,'Diff')
      END IF
C...
C...  PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS, AND THE DIFFERENCE
      TA=EXACT(NCASE)
      DIFF=T(NZ)-TA
      WRITE(NO,2)TIME,T(NZ),TA,DIFF
2     FORMAT(F10.2,3F12.2)
C...
C...  MATLAB PLOTTING OF THE SOLUTION
      CALL PLOTM
      RETURN
      END
      SUBROUTINE PLOTM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=101)
      COMMON/T/     TIME,     NSTOP,     NORUN     1      /Y/    T(NZ)
     2      /F/   TT(NZ)
     3      /S/   TZ(NZ)
     4      /P/        V,        ZL,       TAU,
     5               TEC,       TAC,       T0C
     6      /I/       IP,     NCASE
C...
C...  OPEN A FILE FOR MATLAB PLOTTING
C...
C...  OPEN PLOTTING FILE AFTER THREE RUNS
C...  IF((NORUN.EQ.1).AND.(IP.EQ.1))THEN
C...
C...  OPEN PLOTTING FILE AFTER ONE RUN (NORUN MAY NOT EQUAL 1)
      IF(IP.EQ.1)THEN
         OPEN(1,FILE='h6.out')
      END IF
C...
C...  WRITE THE SOLUTION TO THE FILE FOR MATLAB PLOTTING
      WRITE(1,1)TIME,T(NZ),EXACT(NCASE)
1     FORMAT(3F10.4)
      RETURN
      END
      DOUBLE PRECISION FUNCTION EXACT(NCASE)
C...
C...  FUNCTION EXACT COMPUTES THE EXACT SOLUTION TO THE PDE FOR
C...  A TUBULAR HEAT EXCHANGER
C...
C...     NCASE = 1 - CONSTANT BOUNDARY CONDITION TEMPERATURE
C...
C...     NCASE = 2 - VARIABLE BOUNDARY CONDITION TEMPERATURE
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NZ=101)
      COMMON/T/     TIME,     NSTOP,     NORUN
     1      /Y/   TD(NZ)
     2      /F/   TT(NZ)
     3      /S/   TZ(NZ)
     4      /P/        V,        ZL,       TAU,
     5               TEC,       TAC,       T0C
     6      /I/       IP,      NDUM
C...
C...  PRECOMPUTE SOME EXPONENTIALS
      EXPZ=DEXP(-ZL/(TAU*V))
      EXPT=DEXP(-TIME/TAU)
C...
C...  SOLUTION WITH T LT Z/V
      IF((TIME-ZL/V).LE.0.0D0)THEN
         EXACT=TAC*(1.0D0-EXPT)+
     +         T0C*EXPT
C...
C...  SOLUTION WITH T GT Z/V
      ELSE
     +IF((TIME-ZL/V).GT.0.0D0)THEN
         EXPTD=DEXP(-(TIME-ZL/V)/TAU)
C...
C...     COMPLETE SOLUTION ACCORDING TO BOUNDARY CONDITION
C...
C...     CONSTANT BOUNDARY CONDITION
         IF(NCASE.EQ.1)THEN
            EXACT=TEC*EXPZ+
     +            TAC*((1.0D0-EXPT)-(EXPZ-EXPT))+
     +            T0C*(EXPT-EXPZ*EXPTD)
C...
C...     VARIABLE BOUNDARY CONDITION
         ELSE
     +   IF(NCASE.EQ.2)THEN
            EXACT=TEC*EXPZ*FBC(TIME-ZL/V)+
     +            TAC*((1.0D0-EXPT)-(EXPZ-EXPT))+
     +            T0C*(EXPT-EXPZ*EXPTD)
         END IF
      END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION FBC(TIME)
C...
C...  FUNCTION FBC DEFINES THE VARIABLE (TIME DEPENDENT) BOUNDARY
C...  CONDITION
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  CHARACTERISTIC TIME FOR THE EXPONENTIAL FUNCTION
      TAUBC=1.0D0
C...
C...  EXPONENTIAL INCREASE OF THE BOUNDARY TEMPERATURE WITH TIME
      IF(TIME.LE.0.0D0)FBC=0.0D0
      IF(TIME.GT.0.0D0)FBC=1.0D0-DEXP(-TIME/TAUBC)
      RETURN
      END