      SUBROUTINE INITAL
C...
C...  DYNAMIC MODEL OF A FOUR-PASS SHELL AND TUBE HEAT EXCHANGER WITH
C...  MAPPING OF THE MODEL JACOBIAN MATRIX
C...
C...  A HEAT EXCHANGER IS A DISTRIBUTED SYSTEM IN WHICH CONVECTION AND
C...  TRANSFER OF HEAT ALONG THE EXCHANGER ARE THE PRINCIPAL MODES OF
C...  OPERATION.  THUS THE SPATIAL DISTRIBUTION OF TEMPERATURE WITHIN
C...  THE EXCHANGER MUST BE TAKEN INTO ACCOUNT IN A REALISTIC ANALYSIS.
C...  IF THE TRANSIENT OPERATION OF AN EXCHANGER IS ALSO TO BE CON-
C...  SIDERED, THE MATHEMATICAL MODEL MUST NECESSARILY INVOLVE SPACE
C...  AND TIME, I.E., THE MODEL IS EXPRESSED IN TERMS OF PDES.  ALSO,
C...  MIXING IN THE HEADERS OF THE EXCHANGER (THE ENTRANCE CHAMBER FOR
C...  EACH TUBE PASS) IS MODELED BY ODES WHICH SERVE AS BOUNDARY CON-
C...  DITIONS FOR THE PDES.  THE ANALYSIS OF EXPERIMENTAL DATA FOR HEAT
C...  EXCHANGERS HAS INDICATED THAT THE MIXING IN THE HEADERS CANNOT BE
C...  NEGLECTED AND THUS A MIXED ODE/PDE MODEL IS REQUIRED.  TO ILLUS-
C...  TRATE SOME OF THESE CONCEPTS, CONSIDER THE FOLLOWING SYSTEM OF
C...  NINE PDES AND FIVE PDES FOR A FOUR-PASS, SHELL AND TUBE HEAT
C...  EXCHANGER (SUBSCRIPTS T AND X DENOTE PARTIAL DERIVATIVES WITH
C...  RESPECT TO T AND X, RESPECTIVELY)
C...
C...                   T1  = -C1*T1  + C3*(T5 - T1)                  (1)
C...                     T         X
C...
C...                   T2  =  C1*T2  + C3*(T6 - T2)                  (2)
C...                     T         X
C...
C...                   T3  = -C1*T3  + C3*(T7 - T3)                  (3)
C...                     T         X
C...
C...                   T4  =  C1*T4  + C3*(T8 - T4)                  (4)
C...                     T         X
C...
C...                T5  = C4*(T1 - T5) + C5*(T9 - T5)                (5)
C...                  T
C...
C...                T6  = C4*(T2 - T6) + C5*(T9 - T6)                (6)
C...                  T
C...
C...                T7  = C4*(T3 - T7) + C5*(T9 - T7)                (7)
C...                  T
C...
C...                T8  = C4*(T4 - T8) + C5*(T9 - T8)                (8)
C...                  T
C...
C...           T9  = -C2*T9  + C6*((T5 - T9) + (T6 - T9) +
C...             T         X                                         (9)
C...                               (T7 - T9) + (T8 - T9))
C...
C...  THE NINE DEPENDENT VARIABLES, T1(X,T), T2(X,T), ..., T9(X,T) ARE
C...  TO BE COMPUTED.  C1, C2. ..., C6 ARE GIVEN NUMERICAL CONSTANTS
C...  WHICH REFLECT THE THERMAL CAPACITANCES OF THE LIQUIDS FLOWING
C...  THROUGH THE EXCHANGER, THE THERMAL CAPACITANCE OF THE METAL IN
C...  THE EXCHANGER TUBES AND THE HEAT TRANSFER CHARACTERISTICS.  EQUA-
C...  TIONS (1) TO (4) AND (9) EACH REQUIRE A BOUNDARY CONDITION IN X
C...  AND ALL NINE EQUATIONS REQUIRE AN INITIAL CONDITION IN T.  THESE
C...  WILL BE TAKEN AS
C...
C...                 DT1(0,T)/DT = B1*(TI - T1(0,T))                (10)
C...
C...              DT2(L,T)/DT = B2*(T1(L,T) - T2(L,T))              (11)
C...
C...              DT3(0,T)/DT = B2*(T2(0,T) - T3(0,T))              (12)
C...
C...              DT4(L,T)/DT = B2*(T3(L,T) - T4(L,T))              (13)
C...
C...                           T9(0,T) = TSI                        (14)
C...
C...             T1(X,0) = T2(X,0) = T3(X,0) = T4(X,0) =
C...
C...             T5(X,0) = T6(X,0) = T7(X,0) = T8(X,0) =            (15)
C...
C...                            T9(X,0) = 0
C...
C...  FINALLY, ONE ORDINARY DIFFERENTIAL EQUATION WILL BE ADDED TO MODEL
C...  THE MIXING IN THE HEADER AT THE EXIT OF THE FOURTH PASS
C...
C...                    DT0/DT = B1*(T4(0,T) - T0)                  (16)
C...
C...                             T0(0) = 0                          (17)
C...
C...  EQUATIONS (1) TO (17) CONSTITUTE THE COMPLETE SYSTEM (NINE PDES
C...  AND FIVE ODES).  L, B1 AND B2 ARE THE EXCHANGER LENGTH AND MIXING
C...  TIME CONSTANTS FOR THE HEADERS, RESPECTIVELY.  TI AND TSI ARE THE
C...  ENTERING TEMPERATURES OF THE TUBE SIDE AND SHELL SIDE FLUIDS,
C...  RESPECTIVELY.
C...
C...  THE SOLUTION OF EQUATIONS (1) TO (17) IS BASED ON THE APPROXI-
C...  MATION OF THE DERIVATIVES IN X IN EQUATIONS (1) TO (4) AND (9)
C...  OVER A GRID OF 11 POINTS.  AN APPROXIMATING ODE IS WRITTEN FOR
C...  EACH GRID POINT FOR EACH PDE.  THUS THERE WILL BE A TOTAL OF 11*9
C...  = 99 ODES APPROXIMATING EQUATIONS (1) TO (9).  FINALLY, EQUATION
C...  (16) BRINGS THE TOTAL ODES TO 100.
C...
C...  IN ORDER TO GAIN A PICTURE OF THE 100-ODE SYSTEM STRUCTURE, THE
C...  FOLLOWING PROGRAM ALSO CALLS SUBROUTINES JMAP AND EIGEN TO MAP THE
C...  100-ODE SYSTEM JACOBIAN MATRIX AND COMPUTE THE TEMPORAL EIGEN-
C...  VALUES, RESPECTIVELY.  SUBROUTINE EIGEN IN TURN CALLS EISPACK
C...  ROUTINES BALANC, ELMHES, ELTRAN, HQR2 AND BALBAK.
C...
C...  THE USUAL FUNCTION OF SUBROUTINE INITAL, EVALUATION OF THE INITIAL
C...  CONDITION VECTOR FOR THE DEPENDENT VARIABLE VECTOR IN COMMON/Y/,
C...  IS ACCOMPLISHED BY A BLOCK DATA ROUTINE.  THEREFORE, SUBROUTINE
C...  INITAL IS ESSENTIALLY A DUMMY ROUTINE, AND DOES NOT REQUIRE THE
C...  USUAL COMMON/T/, /Y/ AND /F/.
      RETURN
      END
      BLOCK DATA
C...
C...  THIS BLOCK DATA ROUTINE SETS THE INITIAL CONDITION VECTOR FOR THE
C...  DEPENDENT VARIABLE VECTOR IN /Y/, AND THE PROBLEM PAREMETERS IN
C...  /C/
      COMMON/T/
     1       TIME,     NSTOP,     NORUN
     2      /Y/
     3     T1(11),    T2(11),    T3(11),    T4(11),    T5(11),
     4     T6(11),    T7(11),    T8(11),    T9(11),        T0
     5      /F/
     6    T1T(11),   T2T(11),   T3T(11),   T4T(11),   T5T(11),
     7    T6T(11),   T7T(11),   T8T(11),   T9T(11),       T0T
     8      /S/
     9    T1X(11),   T2X(11),   T3X(11),   T4X(11),   T9X(11)
      COMMON/C/
     1         XL,        C1,        C2,        C3,        C4,
     2         C5,        C6,        B1,        B2,        TI,
     3        TSI,        NG,         N
C...
C...  INITIALIZE THE PROBLEM SYSTEM PARAMETERS
      DATA
     1         XL,        C1,        C2,        C3,        C4,
     2         C5,        C6,        B1,        B2,        TI,
     3        TSI,        NG,         N/
     4      12.E0,    1.98E0,     0.3E0,    0.13E0,   0.476E0,
     5    0.434E0, 0.00963E0,   0.588E0,   0.294E0,     10.E0,
     6       0.E0,        11,       100/
C...
C...  INITIALIZE THE PROBLEM SYSTEM DEPENDENT VARIABLES
      DATA T1,T2,T3,T4,T5,T6,T7,T8,T9,T0/100*0.E+00/
      END
      SUBROUTINE DERV
      COMMON/T/
     1       TIME,     NSTOP,     NORUN
     2      /Y/
     3     T1(11),    T2(11),    T3(11),    T4(11),    T5(11),
     4     T6(11),    T7(11),    T8(11),    T9(11),        T0
     5      /F/
     6    T1T(11),   T2T(11),   T3T(11),   T4T(11),   T5T(11),
     7    T6T(11),   T7T(11),   T8T(11),   T9T(11),       T0T
     8      /S/
     9    T1X(11),   T2X(11),   T3X(11),   T4X(11),   T9X(11)
      COMMON/C/
     1         XL,        C1,        C2,        C3,        C4,
     2         C5,        C6,        B1,        B2,        TI,
     3        TSI,        NG,         N
C...
C...  SET BOUNDARY CONDITION (14)
      T9(1)=TSI
C...
C...  COMPUTE THE FIRST DERIVATIVES WITH RESPECT TO X BY FIVE-POINT
C...  CENTERED DIFFERENCES
      CALL DSS020(0.E+00,XL,NG,T1,T1X, 1.E+00)
      CALL DSS020(0.E+00,XL,NG,T2,T2X,-1.E+00)
      CALL DSS020(0.E+00,XL,NG,T3,T3X, 1.E+00)
      CALL DSS020(0.E+00,XL,NG,T4,T4X,-1.E+00)
      CALL DSS020(0.E+00,XL,NG,T9,T9X, 1.E+00)
C...
C...  ASSEMBLE THE PDES, EQUATIONS (1) TO (9)
      DO 1 I=1,NG
      T1T(I)=-C1*T1X(I) + C3*(T5(I)-T1(I))
      T2T(I)= C1*T2X(I) + C3*(T6(I)-T2(I))
      T3T(I)=-C1*T3X(I) + C3*(T7(I)-T3(I))
      T4T(I)= C1*T4X(I) + C3*(T8(I)-T4(I))
      T5T(I)=             C4*(T1(I)-T5(I))
     1                  + C5*(T9(I)-T5(I))
      T6T(I)=             C4*(T2(I)-T6(I))
     1                  + C5*(T9(I)-T6(I))
      T7T(I)=             C4*(T3(I)-T7(I))
     1                  + C5*(T9(I)-T7(I))
      T8T(I)=             C4*(T4(I)-T8(I))
     1                  + C5*(T9(I)-T8(I))
      T9T(I)=-C2*T9X(I) - 4.E+00*C6*T9(I)
     1                  + C6*(T5(I)+T6(I)+T7(I)+T8(I))
1     CONTINUE
C...
C...  EVALUATE THE ODE DERIVATIVES FROM EQUATIONS (10) AND (13) AND
C...  EQUATION (16).  NOTE THAT THESE DERIVATIVES MUST BE COMPUTED
C...  AFTER THOSE FOR THE PDES
      T1T( 1)=B1*(    TI-T1( 1))
      T2T(NG)=B2*(T1(NG)-T2(NG))
      T3T( 1)=B2*(T2( 1)-T3( 1))
      T4T(NG)=B2*(T3(NG)-T4(NG))
      T9T( 1)=0.E+00
      T0T    =B1*(T4( 1)-    T0)
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/
     1       TIME,     NSTOP,     NORUN
     2      /Y/
     3     T1(11),    T2(11),    T3(11),    T4(11),    T5(11),
     4     T6(11),    T7(11),    T8(11),    T9(11),        T0
     5      /F/
     6    T1T(11),   T2T(11),   T3T(11),   T4T(11),   T5T(11),
     7    T6T(11),   T7T(11),   T8T(11),   T9T(11),       T0T
     8      /S/
     9    T1X(11),   T2X(11),   T3X(11),   T4X(11),   T9X(11)
      COMMON/C/
     1         XL,        C1,        C2,        C3,        C4,
     2         C5,        C6,        B1,        B2,        TI,
     3        TSI,        NG,         N
C...
C...  THIS SECTION OF SUBROUTINE PRINT
C...
C...     (1)  CALLS SUBROUTINE JMAP TO MAP THE JACOBIAN MATRIX OF
C...          100-ODE SYSTEM.
C...
C...     (2)  CALLS A SERIES OF EISPACK ROUTINES VIA SUBROUTINE EIGEN
C...          TO COMPUTE THE TEMPORAL EIGENVALUES OF THE ODE SYSTEM,
C...          AND OPTIONALLY, THE ASSOCIATED EIGENVECTORS.
C...
C...  ABSOLUTE DIMENSIONING OF THE ARRAYS REQUIRED BY SUBROUTINE JMAP
C...  (A, Y, YOLD, F, FOLD) AND SUBROUTINE EIGEN (WI, WR, Z, RW, IW)
      REAL
     1       A(100,100), WR(100), WI(100), Z(100,100), RW(100),
     2                             Y(100),  YOLD(100),
     3                             F(100),  FOLD(100)
      INTEGER                     IW(100)
C...
C...  THE FOLLOWING EQUIVALENCE FACILITATES TRANSFER OF THE DEPENDENT
C...  VARIABLE VECTOR IN /Y/ TO ARRAY Y, AND THE DERIVATIVE VECTOR IN
C...  IN /F/ TO ARRAY F FOR USE IN THE SUBSEQUENT CALL TO SUBROUTINE
C...  JMAP
      EQUIVALENCE (T1(1),Y(1)),(T1T(1),F(1))
C...
C...  MAP THE JACOBIAN MATRIX OF THE ODE SYSTEM DEFINED IN SUBROUTINE
C...  DERV, AND COMPUTE ITS TEMPORAL EIGENVALUES
      IF(TIME.GT.1.E+00)GO TO 5
      CALL JMAP(N,A,Y,YOLD,F,FOLD)
C...
C...  SUBROUTINE EIGEN (PART OF DSS/2) CALLS A SERIES OF EISPACK ROU-
C...  TINES TO COMPUTE THE TEMPORAL EIGENVALUES, AND OPTIONALLY, THE
C...  EIGENVECTORS OF THE 100-ODE SYSTEM JACOBIAN MATRIX
      CALL EIGEN(N,A,WR,WI,Z,RW,IW)
C...
C...  PRINT A HEADING FOR THE NUMERICAL SOLUTION
      WRITE(NO,3)
3     FORMAT(1H1,9X,50HMIXED ORDINARY/PARTIAL DIFFERENTIAL EQUATION MODE
     1L,/)
      WRITE(NO,1)C1,C2,C3,C4
1     FORMAT(
     13X,4HC1 = ,E10.3,3X,4HC2 = ,E10.3,3X,4HC3 = ,E10.3,
     23X,4HC4 = ,E10.3)
      WRITE(NO,2)C5,C6,B1,B2
2     FORMAT(
     13X,4HC5 = ,E10.3,3X,4HC6 = ,E10.3,3X,4HB1 = ,E10.3,
     23X,4HB2 = ,E10.3,//)
      WRITE(NO,4)
4     FORMAT(
     1 6X,4HTIME,7X,7HT1(X=0),7X,8HT4(X=12),8X,2HT0,10X,
     2 8HT9(X=12))
C...
C...  PRINT THE SOLUTION
5     WRITE(NO,6)TIME,T1(1),T4(NG),T0,T9(NG)
6     FORMAT(E12.3,4E14.5)
      RETURN
      END
