*DECK PRO5P1
      SUBROUTINE INITAL
C...
C...  DYNAMICS OF A HUMIDIFICATION COLUMN
C...
C...  A HUMIDIFICATION COLUMN (A PACKED COLUMN WITH WATER FLOWING DOWN
C...  AND AIR FLOWING UP) IS DIAGRAMMED BELOW
C...
C...                     L,TL(ZL,T)   V,Y(ZL,T),
C...                                   TG(ZL,T)
C...                          .           +
C...                          .           .
C...                          .           .
C...                          +           .
C...                          .  .  .  .  . Z = ZL
C...                        .              .
C...                       .                .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      . L,TL       V,Y,TG.
C...                      .  .            +  .
C...                      .  .            .  .
C...                      .  +            .  .
C...                      .................... Z+DZ
C...                      .      .           .
C...                      . H20  .           .
C...                      .PHASE . ----+ H20 .
C...                      .      .    AIR    .
C...                      .      .   PHASE   .
C...                      ....................  Z
C...                      .  .            +  .
C...                      .  .            .  .
C...                      .  +            .  .
C...                      . H20          AIR .
C...                      .STREAM      STREAM.
C...                      .(DOWN)       (UP) .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                       .                .
C...                        .              .
C...                          .  .  .  .  . Z = 0
C...                          .           +
C...                          .           .
C...                          .           .
C...                          +           .
C...                      L,TL(0,T)   V,Y(0,T),
C...                                   TG(0,T)
C...
C...  THE VARIABLES AND PARAMETERS FOR THIS SYSTEM ARE
C...
C...     TG      GAS TEMPERATURE (UNITS OF THE MODEL VARIABLES ARE
C...             DEFINED BELOW WITH THE NUMERICAL VALUES OF THE MODEL
C...             PARAMETERS AND AUXILIARY CONDITIONS)
C...
C...     EV      INTERNAL ENERGY OF THE GAS STREAM
C...             = (CVA*TG + Y*(CVV*TG + DHVAP))
C...
C...     EP      ENTHALPY OF THE GAS STRED
C...             = (CPA*TG + Y*(CPV*TG + DHVAP))
C...
C...     TL      LIQUID TEMPERATURE
C...
C...     Y       MOLE RATIO OF H20 IN THE GAS
C...
C...     YS      MOLE RATIO OF H20 IN THE GAS WHICH WOULD BE IN
C...             EQUILIBRIUM WITH THE LIQUID AT TEMPERATURE TL
C...
C...     T       TIME
C...
C...     Z       AXIAL POSITION ALONG THE COLUMN
C...
C...     ZL      LENGTH OF THE COLUMN
C...
C...     V       DRY AIR MOLAR FLOW RATE
C...
C...     G       DRY AIR MOLAR HOLDUP
C...
C...     S       COLUMN CROSS SECTIONAL AREA
C...
C...     KY      MASS TRANSFER COEFFICIENT
C...
C...     AV      HEAT AND MASS TRANSFER AREAS PER UNIT VOLUME OF
C...             COLUMN
C...
C...     H       HEAT TRANSFER COEFFICIENT
C...
C...     CVV     SPECIFIC HEAT OF WATER VAPOR AT CONSTANT VOLUME
C...
C...     CPV     SPECIFIC HEAT OF WATER VAPOR AT CONSTANT PRESSURE
C...
C...     CVA     SPECIFIC HEAT OF DRY AIR AT CONSTANT VOLUME
C...
C...     CPA     SPECIFIC HEAT OF DRY AIR AT CONSTANT PRESSURE
C...
C...     CL      SPECIFIC HEAT OF WATER
C...
C...     L       LIQUID (H20) MOLAR FLOW RATE
C...
C...     H       LIQUID (H20) MOLAR HOLDUP
C...
C...     DHVAP   HEAT OF VAPORIZATION OF WATER
C...
C...  THE DEPENDENT VARIABLES OF PARTICULAR INTEREST ARE Y(Z,T), TG(Z,T)
C...  TL(Z,T).  IN ORDER TO COMPUTE THESE VARIABLES, BASIC CONSERVATION
C...  PRINCIPLES ARE APPLIED TO THE AIR AND WATER STREAMS WITHIN THE
C...  COLUMN OVER THE INCREMENTAL SECTION OF LENGTH DZ (SEE THE DIAGRAM)
C...  AND THEN THE LIMIT OF DZ ---+ 0 GIVES THE MODEL DIFFERENTIAL EQUA-
C...  TIONS.  SINCE THE MODEL IS DYNAMIC, TIME, T, AND AXIAL POSITION,
C...  Z, APPEAR IN THE EQUATIONS AS INDEPENDENT VARIABLES AND THEREFORE
C...  PARTIAL DIFFERENTIAL EQUATIONS (PDES) MUST BE USED.
C...
C...  THE PDES ARE DERIVED AS FOLLOWS.  A H20 BALANCE ON THE GAS PHASE
C...  GIVES
C...
C...     (Y*G*S*DZ)  = V*Y  - V*Y     + KY*AV*(YS - Y)*S*DZ     (1A)
C...               T      Z      Z+DZ
C...
C...  WHERE THE SUBSCRIPT T ON THE LHS DENOTES A PARTIAL DERIVATIVE WITH
C...  RESPECT TO TIME, T.
C...
C...  THE PHYSICAL INTERPRETATION AND UNITS OF EACH OF THE TERMS IN
C...  EQUATION (1A) ARE INDICATED BELOW
C...
C...     (Y*G*S*DZ)
C...               T
C...
C...        = (MOLS H20/MOL DRY AIR)*(MOLS DRY AIR/CM**3)*(CM**2)*
C...          CM*(1/HR) = MOLS H20/HR ACCUMULATING IN THE GAS PHASE
C...          WITHIN THE INCREMENTAL SECTION OF LENGTH DZ (SEE THE
C...          DIAGRAM BELOW)
C...
C...     V*Y
C...        Z
C...
C...        = (MOLS DRY AIR/HR)*(MOLS H20/MOL DRY AIR) = MOLS H20/HR
C...          FLOWING INTO THE INCREMENTAL SECTION AT Z
C...
C...     V*Y
C...        Z+DZ
C...
C...        = (MOLS DRY AIR/HR)*(MOLS H20/MOL DAY AIR) = MOLS H20/HR
C..           FLOWING OUT OF THE INCREMENTAL SECTION AT Z+DZ
C...
C...     KY*AV*(YS - Y)*S*DZ
C...
C...        = (MOLS DRY AIR/HR-CM**2)*(CM**2/CM**3)*(MOLS H20/MOL DRY
C...           AIR)*(CM**2)*(CM) = MOLS H20/HR TRANSFERRED TO OR FROM
C...           THE AIR STREAM IN THE INCREMENTAL SECTION (THE DIRECTION
C...           OF TRANSFER DEPENDS ON THE SIGN OF (YS - Y))
C...
C...  IN THE LIMIT AS DZ ---+ 0, EQUATION (1A) BECOMES
C...
C...     Y   = -(V/(G*S))*Y  + (KY*AV/G)*(YS - Y)                (1)
C...      T                Z
C...
C...  AN ENERGY BALANCE ON THE GAS PHASE GIVES
C...
C...     ((CVA*TG + Y*(CVV*TG + DHVAP))*G*S*DZ)  =
C...                                           T
C...
C...     (CPA*TG + Y*(CPV*TG + DHVAP))*V  -
C...                                    Z                       (2A)
C...
C...     (CPA*TG + Y*(CPV*TG + DHVAP))*V     -
C...                                    Z+DZ
C...
C...     H*AV*(TL - TG)*S*DZ + KY*AV*(YS - Y)*S*DZ*(CVV*TG + DHVAP)
C...
C...  THE PHYSICAL INTERPRETATION AND UNITS OF EACH OF THE TERMS IN
C...  EQUATION (2A) ARE INDICATED BELOW
C...
C...     ((CVA*TG + Y*(CVV*TG + DHVAP))G*S*DZ)
C...                                          T
C...
C...        = (CAL/MOL DRY AIR-C)*(C) + (MOLS H20/MOL DRY AIR)*
C...          (CAL/MOL H20-C)*(C) + CAL/MOL H20))*(MOLS DRY AIR/
C...          CM**3)*(CM**2)*CM*(1/HR) = CAL/HR ACCUMULATING IN THE
C...          GAS PHASE WITHIN THE INCREMENTAL SECTION
C...
C...     (CPA*TG + Y*(CPV*TG + DHVAP))*V
C...                                    Z
C...
C...        = (CAL/MOL DRY AIR-C)*(C) + (MOLS H20/MOL DRY AIR)*
C...          (CAL/MOL H20-C)*(C) + CAL/MOL H20))*(MOLS DRY AIR/HR)
C...        = CAL/HR FLOWING INTO THE INCREMENTAL SECTION AT Z
C...
C...     (CPA*TG + Y*(CPV*TG + DHVAP))*V
C...                                    Z+DZ
C...
C...        = (CAL/MOL DRY AIR-C)*(C) + (MOLS H20/MOL DRY AIR)*
C...          (CAL/MOL H20-C)*(C) + CAL/MOL H20))*(MOLS DRY AIR/HR)
C...          = CAL/HR FLOWING OUT OF THE INCREMENTAL SECTION AT Z+DZ
C...
C...     H*AV*(TL - TG)*S*DZ
C...
C...        = (CAL/HR-CM**2-C)*(CM**2/CM**3)*(C)*(CM**2)*CM = CAL/HR
C...          TRANSFERRED BETWEEN THE H20 AND GAS PHASES DUE TO TEMP-
C...          ERATURE DIFFERENCES (THE DIRECTION OF TRANSFER DEPENDS
C...          ON THE SIGN OF (TL - TG))
C...
C...     KY*AV*(YS - Y)*S*DZ*(CVV*TG + DHVAP)
C...
C...        = (MOLS DRY AIR/HR-CM**2)*(CM**2/CM**3)*(MOLS H20/MOL DRY
C...          AIR)*(CM**2)*(CM)*((CAL/MOL H20-C)*C + CAL/MOL H20) =
C...          CAL/HR TRANSFERRED BETWEEN THE H20 AND GAS PHASES DUE TO
C...          DIFFERENCES IN HUMIDITY
C...
C...  IN THE LIMIT AS DZ ---+ 0, EQUATION (2A) BECOMES
C...
C...     EV  = -(V/(G*S))*EP  + (H*AV/G)*(TL - TG)
C...       T                Z                                    (2)
C...
C...                          + (KY*AV/G)*(YS - Y)*(CVV*TG + DHVAP)
C...
C...  NOTE THE ENERGY GROUPS, EV (ENERGY AT CONSTANT VOLUME) AND EP
C...  (ENERGY AT CONSTANT PRESSURE) HAVE BEEN USED TO SIMPLIFY THE
C...  TERMS IN EQUATION (2).  EV AND EP ARE DEFINED IN THE TABLE OF
C...  NOMENCLATURE GIVEN EARLIER.
C...
C...  AN ENERGY BALANCE ON THE LIQUID GIVES
C...
C...     (CL*TL*H*S*DZ)  = (CL*TL*L)     - (CL*TL*L)
C...                   T            Z+DZ            Z           (3A)
C...
C...     -H*AV*(TL - TG)*S*DZ - KY(AV*(YS - Y)*S*DZ*(CVV*TG + HDVAP)
C...
C...  THE PHYSICAL INTERPRETATION AND UNITS OF EACH OF THE TERMS IN
C...  EQUATION (3A) ARE GIVEN BELOW
C...
C...     (CL*TL*H*S*DZ)
C...                   T
C...
C...        (CAL/MOL H20-C)*(C)*(MOLS H20/CM**3)*(CM**2)*CM*(1/HR)
C...        = CAL/HR ACCUMULATING IN THE H20 PHASE WITHIN THE
C...        INCREMENTAL SECTION OF LENGTH DZ
C...
C...     (CL*TL*L)
C...              Z+DZ
C...
C...        = (CAL/MOL H20-C)*(C)*(MOLS H20/HR) = CAL/HR FLOWING INTO
C...          THE INCREMENTAL SECTION AT Z+DZ
C...
C...     (CL*TL*L)
C...              Z
C...
C...        = (CAL/MOL H20-C)*(C)*(MOLS H20/HR) = CAL/HR FLOWING OUT OF
C...          THE INCREMENTAL SECTION AT Z
C...
C...     H*AV*(TL - TG)*S*DZ
C...
C...        = (CAL/HR-CM**2-C)*(CM**2/CM**3)*(C)*(CM**2)*CM = CAL/HR
C...          TRANSFERRED BETWEEN THE H20 AND GAS PHASES DUE TO TEMP-
C...          ERATURE DIFFERENCES (THE DIRECTION OF TRANSFER DEPENDS
C...          ON THE SIGN OF (TL - TG))
C...
C...     KY*AV*(YS - Y)*S*DZ*(CVV*TG + DHVAP)
C...
C...        = (MOLS DRY AIR/HR-CM**2)*(CM**2/CM**3)*(MOLS H20/MOL DRY
C...          AIR)*(CM**2)*(CM)*((CAL/MOL H20-C)*C + CAL/MOL H20) =
C...          CAL/HR TRANSFERRED BETWEEN THE H20 AND GAS PHASES DUE TO
C...          DIFFERENCES IN HUMIDITY
C...
C...  IN THE LIMIT AS DZ ---+ 0, EQUATION (3A) BECOMES
C...
C...     TL  =  (L/(H*S))*TL  - (H*AV/(CL*H))*(TL - TG)
C...       T                Z                                    (3)
C...
C...            - (KY*AV/(CL*H))*(YS - Y)*(CVV*TG + DHVAP)
C...
C...  AN EXTENSION OF THE PRECEDING SYSTEM IS ALSO SIMULATED IN THE
C...  FOLLOWING CODING.  SPECIFICALLY, THE EXITING WATER TEMPERATURE
C...  TL(0,T), IS SENSED, COMPARED WITH A SET POINT VALUE, TLSET,
C...  AND THE DIFFERENCE (ERROR), E = TL(0,T) - TLSET, IS USED TO
C...  ADJUST THE ENTERING AIR FLOW RATE, V, TO MOVE TL(0,T) TO TLSET.
C...  THE ADJUSTMENT OF V IS ACHIEVED WITH A PROPORTIONAL-INTEGRAL (PI)
C...  CONTROLLER.  THE FOLLOWING DIAGRAM ILLUSTRATES THE CONTROLLER
C...  AT THE BOTTOM OF THE HUMIDIFICATION COLUMN
C...
C...                     L,TL(ZL,T)   V,Y(ZL,T),
C...                                   TG(ZL,T)
C...                          .           +
C...                          .           .
C...                          .           .
C...                          +           .
C...                          .  .  .  .  . Z = ZL
C...                        .              .
C...                       .                .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      . L,TL       V,Y,TG.
C...                      .  .            +  .
C...                      .  .            .  .
C...                      .  +            .  .
C...                      .................... Z+DZ
C...                      .      .           .
C...                      . H20  .           .
C...                      .PHASE . ----+ H20 .
C...                      .      .    AIR    .
C...                      .      .   PHASE   .
C...                      ....................  Z
C...                      .  .            +  .
C...                      .  .            .  .
C...                      .  +            .  .
C...                      . H20          AIR .
C...                      .STREAM      STREAM.
C...                      .(DOWN)       (UP) .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                      .                  .
C...                       .                .
C...                        .              .
C...                          .  .  .  .  . Z = 0
C...                          .           +             PI
C...                          .           .         CONTROLLER
C...                          .           .          .........
C...            TEMPERATURE .....         .          .       .
C...              SENSOR    .   .......... ........... KC,TI .+--- TLSET
C...                        .....         .          .       .
C...                          .           .          .........
C...                          .           .              .
C...                          .          ...   ..        .
C...                          .           ...... .........
C...                          .          ...   ..
C...                          .           . CONTROL VALVE
C...                          +           .
C...                      L,TL(0,T)   V,Y(0,T),
C...                                   TG(0,T)
C...
C...  THE PURPOSE OF THE SIMULATION IS ESSENTIALLY TO DETERMINE THE
C...  VALUES OF THE CONTROLLER GAIN, KC, AND INTEGRAL TIME, TI, SO
C...  THAT THE EXITING LIQUID TEMPERATURE, TL(0,T), IS CONTROLLED CLOSE
C...  TO THE SET POINT, TLSET.
C...
C...  THE CONTROLLER EQUATIONS ARE
C...
C...     V = CV*SQRT(DP)*X                                         (4)
C...
C...     X = XSS + KC*(E + (1/TI)INT(E)DT)
C...
C...     E = TL(0,T) - TLSET
C...
C...     X       CONTROL VALVE STEM POSITION
C...
C...     XSS     X AT STEADY STATE
C...
C...     KC      CONTROLLER GAIN
C...
C...     TI      CONTROLLER INTEGRAL TIME
C...
C...     CVDP    PRODUCT OF THE CONTROL VALVE CONSTANT AND SQUARE
C...             ROOT OF THE PRESSURE DROP ACROSS THE VALVE
C...
C...     TLSET   CONTROLLER SET POINT
C...
C...  THE MODEL PARAMETERS AND AUXILIARY CONDITIONS IN ENGLISH UNITS ARE
C...  LISTED BELOW, FOLLOWED BY EQUIVALENT VALUES IN METRIC UNITS.  THE
C...  PROGRAM CAN BE EXECUTED WITH EITHER SET BY SELECTING THE VALUE OF 
C...  THEVARIABLE METRIC IN THE PROGRAM (METRIC = 0 FOR ENGLISH UNITS,
C...  METRIC = 1 FOR METRIC UNITS)
C...
C...  ******************************************************************
C...
C...   ENGLISH UNITS
C...
C...   H*AV = 50 BTU/(HR-FT**3-F)
C...
C...   KY*AV = 4.1 LB MOLS DRY AIR/(HR-FT**3)
C...
C...   ZL = 8 FT, CL = 18 BTU/(LB MOL-F), DHVAP = 17000 BTU/(LB MOL)
C...
C...   CPV = 7.3 BTU/(LB MOL-F), CPA = 7.3 BTU/(LB MOL-F)
C...
C...   CVV = CPV - R = 7.3 - 1.987 = 5.3 BTU/(LB MOL-F)
C...
C...   CVA = CPA - R = 7.3 - 1.987 = 5.3 BTU/(LB MOL-F)
C...
C...   L = 50 LB MOLS/HR, S = 1 FT**2
C...
C...   H = 0.5 LB MOLS/FT**3, G = 0.01 LB MOLS/FT**3
C...
C...   P = 10**(7.96681 - 3002.4/(378.4 + TL)) MM HG (TL IS IN F)
C...
C...   YS = P/(760 - P) LB MOLS H20/LB MOL DRY AIR
C...
C...   TG(0,T) = 110 F, TL(ZL,T) = 110 F, Y(0,T) = 0.01 LB MOLS H20/
C...                                              LB MOL DRY AIR)
C...   TG(Z,0) = TL(Z,0) = 110 F, Y(Z,0) = 0.01
C...
C...   KC = 0.02 (1/F), TI = 0.1 HR, CV*(DP**0.5) = 50 LB MOLS/HR,
C...
C...   XSS = 0.5, TLSET = 90 F
C...
C...  ******************************************************************
C...
C...   METRIC UNITS
C...
C...   H*AV = 0.8015 CAL/(HR-CM**3-C)
C...
C...   KY*AV = 0.06573 GM MOLS DRY AIR/(HR-CM**3)
C...
C...   ZL = 243.8 CM, CL = 18 CAL/(GM MOL-C), DHVAP = 9443.6 CAL/GM MOL
C...
C...   CPV = 7.3 CAL/(GM MOL-C), CPA = 7.3 CAL/(GM MOL-C)
C...
C...   CVV = CPV - R = 7.3 - 1.987 = 5.3 CAL/(GM MOL-C)
C...
C...   CVA = CPA - R = 7.3 - 1.987 = 5.3 CAL/(GM MOL-C)
C...
C...   L = 22700 GM MOLS/HR, S = 929.03 CM**2
C...
C...   H = 0.008015 GM MOLS/CM**3, G = 0.0001603 GM MOLS/CM**3
C...
C...   P = 10**(7.96681 - 3002.4/(378.4 + 1.8*TL + 32)) MM HG
C...                                                    (TL IS IN C)
C...
C...   YS = P/(760 - P) GM MOLS H20/GM MOL DRY AIR
C...
C...   TG(0,T) = 43.33 C, TL(ZL,T) = 43.33 C, Y(0,T) = 0.01 GM MOL H20/
C...                                                  GM MOL DRY AIR)
C...
C...   TG(Z,0) = TL(Z,0) = 43.33 C, Y(Z,0) = 0.01
C...
C...   KC = 0.036 (1/C), TI = 0.1 HR, CV*(DP**0.5) = 22700 GM MOLS/HR
C...
C...   XSS = 0.5, TLSET = 32.22 C
C...
C...  ******************************************************************
C...
C...  THE CONVERSION FACTORS USED TO GO FROM ONE SET OF UNITS TO THE
C...  OTHER ARE
C...
C...     1 BTU = 252.2 CAL
C...
C...     1 BTU/LB MOL-F = 1 CAL/GM MOL-C
C...
C...     1 LB = 454 GM
C...
C...     1 LB MOL = 454 GM MOL
C...
C...     1 F = 1/1.8 C
C...
C...     1 FT = 30.48 CM
C...
C...     1 FT**3 = 28316.8 CM**3
C...
C...     TF = 1.8*TC + 32
C...
C...  THE EXITING GAS COMPOSITION AND TEMPERATURE, Y(ZL,T) AND TG(ZL,T),
C...  THE EXITING LIQUID TEMPERATURE, TL(0,T), AND THE GAS FLOW RATE,
C...  G, ARE OF PARTICULAR INTEREST IN THIS SIMULATION.
C...
C...  COMMON AREA TO ESTABLISH LINKAGE WITH OTHER SUBROUTINES
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  Y(11), EV(11), TL(11),     EI
     2      /F/ YT(11),EVT(11),TLT(11),      E
     3      /S/ YZ(11),EPZ(11),TLZ(11)
     4      /C/      V,      L,      G,      H,   HLAV,   KYAV,
     5              CL,    CPV,    CPA,    CVV,    CVA,     ZL,
     6           DHVAP,      S,     P1,     P2,     P3,     P4,
     7              P5,     P6,     P7, TG(11), EP(11), YS(11),
     8              KC,     TI,   CVDP,    XSS,  TLSET,      X,
     9              IP,      N, METRIC
      REAL           L,   KYAV,     KC
C...
C...  SELECT ENGLISH (METRIC = 0) OR METRIC (METRIC = 1) UNITS
      IF(NORUN.EQ.1)METRIC=0
      IF(NORUN.EQ.2)METRIC=1
C...
C...  SET THE MODEL PARAMETERS IN THE UNITS SELECTED
C...
C...  ENGLISH UNITS
      IF(METRIC.EQ.0)THEN
      N=11
      L=50.
      G=0.01
      H=0.5
      HLAV=50.
      KYAV=4.1
      CL=18.
      CPV=7.3
      CPA=7.3
      CVV=5.3
      CVA=5.3
      ZL=8.
      DHVAP=17000.
      S=1.
C...
C...  SET THE CONTROLLER PARAMETERS
      KC=0.02
      TI=0.1
      CVDP=50.
      XSS=0.5
      TLSET=90.
C...
C...  METRIC UNITS
      ELSE
      N=11
      L=22700.
      G=0.0001603
      H=0.008015
      HLAV=0.8015
      KYAV=0.06573
      CL=18.
      CPV=7.3
      CPA=7.3
      CVV=5.3
      CVA=5.3
      ZL=243.8
      DHVAP=9443.6
      S=929.03
C...
C...  SET THE CONTROLLER PARAMETERS
      KC=0.036
      TI=0.1
      CVDP=22700.
      XSS=0.5
      TLSET=32.22
      END IF
C...
C...  COMPUTE SOME CONSTANTS FOR USE IN OTHER SUBROUTINES
      P2=KYAV/G
      P3=HLAV/G
      P4=L/(H*S)
      P5=HLAV/(CL*H)
      P6=KYAV/(CL*H)
C...
C...  INITIAL CONDITIONS
      DO 1 I=1,N
      Y(I)=0.01
      IF(METRIC.EQ.0)THEN
      TL(I)=110.
      TG(I)=110.
      ELSE
      TL(I)=43.33
      TG(I)=43.33
      END IF
      EV(I)=CVA*TG(I)+Y(I)*(CVV*TG(I)+DHVAP)
1     CONTINUE
      EI=1.0E-05
C...
C...  INITIALIZE THE MODEL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  Y(11), EV(11), TL(11),     EI
     2      /F/ YT(11),EVT(11),TLT(11),      E
     3      /S/ YZ(11),EPZ(11),TLZ(11)
     4      /C/      V,      L,      G,      H,   HLAV,   KYAV,
     5              CL,    CPV,    CPA,    CVV,    CVA,     ZL,
     6           DHVAP,      S,     P1,     P2,     P3,     P4,
     7              P5,     P6,     P7, TG(11), EP(11), YS(11),
     8              KC,     TI,   CVDP,    XSS,  TLSET,      X,
     9              IP,      N, METRIC
      REAL           L,   KYAV,     KC
C...
C...  CONTROLLER EQUATIONS
      E=TL(1)-TLSET
      X=XSS+KC*(E+(1./TI)*EI)
      IF(X.LT.0.)X=0.
      IF(X.GT.1.)X=1.
      V=CVDP*X
      P1=V/(G*S)
C...
C...  MODEL ALGEBRA
      DO 1 I=1,N
      TG(I)=(EV(I)-Y(I)*DHVAP)/(CVA+Y(I)*CVV)
      EP(I)=CPA*TG(I)+Y(I)*(CPV*TG(I)+DHVAP)
      IF(METRIC.EQ.0)THEN
      P=10.**(7.96681-3002.4/(378.4+TL(I)))
      ELSE
      P=10.**(7.96681-3002.4/(378.4+1.8*TL(I)+32.))
      END IF
      YS(I)=P/(760.-P)
1     CONTINUE
C...
C...  BOUNDARY CONDITIONS
      Y(1)=0.01
      IF(METRIC.EQ.0)THEN
      TL(N)=110.
      TG(1)=110.
      ELSE
      TL(N)=43.33
      TG(1)=43.33
      END IF
      EP(1)=CPA*TG(1)+Y(1)*(CPV*TG(1)+DHVAP)
C...
C...  CONSTRAINTS FOR THE DEPENDENT VARIABLES AND ASSOCIATED
C...  VARIABLES
      DO 3 I=1,N
      IF( Y(I).LT. 0.) Y(I)=0.
      IF(YS(I).LT. 0.)YS(I)=0.
3     CONTINUE
C...
C...  SPATIAL DERIVATIVES
      CALL DSS020(0.,ZL,N, Y, YZ, 1.)
      CALL DSS020(0.,ZL,N,EP,EPZ, 1.)
      CALL DSS020(0.,ZL,N,TL,TLZ,-1.)
C...
C...  PARTIAL DIFFERENTIAL EQUATIONS
      DO 2 I=1,N
      YT(I)=-P1*YZ(I)+P2*(YS(I)-Y(I))
      P7=CVV*TG(I)+DHVAP
      EVT(I)=-P1*EPZ(I)+P3*(TL(I)-TG(I))+P2*(YS(I)-Y(I))*P7
      TLT(I)= P4*TLZ(I)-P5*(TL(I)-TG(I))-P6*(YS(I)-Y(I))*P7
2     CONTINUE
C...
C...  ZERO THE TIME DERIVATIVES TO MAINTAIN THE BOUNDARY CONDITIONS
      YT(1)=0.
      EVT(1)=0.
      TLT(N)=0.
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  Y(11), EV(11), TL(11),     EI
     2      /F/ YT(11),EVT(11),TLT(11),      E
     3      /S/ YZ(11),EPZ(11),TLZ(11)
     4      /C/      V,      L,      G,      H,   HLAV,   KYAV,
     5              CL,    CPV,    CPA,    CVV,    CVA,     ZL,
     6           DHVAP,      S,     P1,     P2,     P3,     P4,
     7              P5,     P6,     P7, TG(11), EP(11), YS(11),
     8              KC,     TI,   CVDP,    XSS,  TLSET,      X,
     9              IP,      N, METRIC
      REAL           L,   KYAV,     KC
C...
C...  DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION, INDIVIDUAL TERMS
C...  IN THE PDES
      DIMENSION       TLP(101),       TGP(101),        YP(101),
     1                 VP(101),        XP(101),       ERP(101),
     2                 TP(101),        TEP(11),     TGLP(2,11)
C...
C...  DIMENSION THE ARRAYS FOR SUBROUTINE JMAP AND EIGEN
      DIMENSION       A(34,34),         WR(34),         WI(34),
     1                Z(34,34),         RW(34),         IW(34),
     2                                  SV(34),      SVOLD(34),
     3                                   F(34),       FOLD(34)
C...
C...  EQUIVALENCE THE STATE VARIABLE VECTOR, S(34), TO THE DEPPENDENT
C...  VARIABLE VECTORS IN COMMON/Y/
      EQUIVALENCE (Y(1),SV(1)),   (YT(1),F(1))
C...
C...  PRINT A DETAILED SOLUTION (THE INDIVIDUAL TERMS IN THE PDES AND
C...  THE DEPENDENT VARIABLES) FOR IP = 1, 51, 101 CALLS TO SUBROUTINE
C...  PRINT.  FIRST, UPDATE ALL OF THE MODEL CALCULATIONS BY A CALL
C...  TO SUBROUTINE DERV (THE INITIAL-VALUE INTEGRATION THROUGH COMMON
C...  /Y/ AND /F/ IS DONE CORRECTLY WITHOUT THIS UPDATE, BUT SOME OF THE
C...  MODEL VARIABLES MAY BE ONE INTEGRATION TIME STEP BEHIND IN TIME,
C...  T, WHEN THEY ARE PRINTED OUT IN SUBROUTINE PRINT IF THE FOLLOWING
C...  CALL TO SUBROUTINE DERV IS NOT EXECUTED)
      IP=IP+1
      IF((IP.NE.1).AND.(IP.NE.51).AND.(IP.NE.101))GO TO 100
      CALL DERV
C...
C...  ***************************************************************
C...  GAS-PHASE MATERIAL BALANCE
C...  ************************************************************(1)
C...
C...     -(V/(G*S))*Y
C...                 Z
         DO 1 I=1,N
         TEP(I)=-P1*EPZ(I)
1        CONTINUE
         N2=1
         WRITE(NO,2)T,(TEP(I),I=1,N,N2)
2        FORMAT( 1H1,5H T = ,F7.3,//,
     1          25H -(V/(G*S))*YZ           ,5E11.3,/,13X,6E11.3,//)
C...
C...  ************************************************************(2)
C...
C...     (KY*AV/G)*(YS - Y)
         DO 3 I=1,N
         TEP(I)=P3*(TL(I)-TG(I))
3        CONTINUE
         WRITE(NO,4)(TEP(I),I=1,N,N2)
4        FORMAT(25H (KA*AV/G)*(YS - Y)      ,5E11.3,/,13X,6E11.3,//)
C...
C...  ************************************************************(3)
C...
C...     Y
C...      T
         WRITE(NO,5)(YT(I),I=1,N,N2)
5        FORMAT(25H YT                      ,5E11.3,/,13X,6E11.3,//)
C...
C...  ************************************************************(4)
C...
C...     Y
         WRITE(NO,6)(Y(I),I=1,N,N2)
6        FORMAT(25H Y                       ,5E11.3,/,13X,6E11.3,//)
C...
C...  ***************************************************************
C...  EQUILIBRIUM RELATIONSHIP
C...  ************************************************************(5)
C...
C...     YS
         WRITE(NO,7)(YS(I),I=1,N,N2)
7        FORMAT(25H YS                      ,5E11.3,/,13X,6E11.3,//)
C...
C...  ***************************************************************
C...  GAS-PHASE ENERGY BALANCE
C...  ************************************************************(6)
C...
C...     -(V/(G*S))*EP  = -(V/(G*S))*(CPA*TG + Y*(CPV*TG + DHVAP))
C...                  Z                                           Z
         DO 8 I=1,N
         TEP(I)=-P1*EPZ(I)
8        CONTINUE
         WRITE(NO,9)(TEP(I),I=1,N,N2)
9        FORMAT(25H -(V/(G*S))*EPZ          ,5E11.3,/,13X,6E11.3,//)
C...
C...  ************************************************************(7)
C...
C...     (H*AV/G)*(TL - TG)
         DO 10 I=1,N
         TEP(I)=P3*(TL(I)-TG(I))
10       CONTINUE
         WRITE(NO,11)(TEP(I),I=1,N,N2)
11       FORMAT(25H -(H*AV/G)*(TL - TG)     ,5E11.3,/,13X,6E11.3,//)
C...
C...  ************************************************************(8)
C...
C...     (KY*AV/G)*(YS - Y)*(CVV*TG + DHVAP)
         DO 12 I=1,N
         TEP(I)=-P2*(YS(I)-Y(I))*(CVV*TG(I)+DHVAP)
12       CONTINUE
         WRITE(NO,13)(TEP(I),I=1,N,N2)
13       FORMAT(25H (KY*AV/G)*(YS - Y)*(...),5E11.3,/,13X,6E11.3,//)
C...
C...  ************************************************************(9)
C...
C...     EV  = (CVA*TG + Y*(CVV*TG + DHVAP))
C...       T                                T
         WRITE(NO,14)(EVT(I),I=1,N,N2)
14       FORMAT(25H EVT                     ,5E11.3,/,13X,6E11.3,//)
C...  ************************************************************(10)
C...
C...     TG
         WRITE(NO,15)(TG(I),I=1,N,N2)
15       FORMAT(25H TG                      ,5E11.3,/,13X,6E11.3,//)
C...
C...  ****************************************************************
C...  LIQUID-PHASE ENERGY BALANCE
C...  ************************************************************(11)
C...
C...     -(L/(H*S))*TL
C...                  Z
         DO 16 I=1,N
         TEP(I)=P4*TLZ(I)
16       CONTINUE
         WRITE(NO,17)(TEP(I),I=1,N,N2)
17       FORMAT(25H -(L/(H*S))*TLZ          ,5E11.3,/,13X,6E11.3,//)
C...
C...  ************************************************************(12)
C...
C...     -(H*AV/(CL*H))*(TL - TG)
         DO 18 I=1,N
         TEP(I)=-P5*(TL(I)-TG(I))
18       CONTINUE
         WRITE(NO,19)(TEP(I),I=1,N,N2)
19       FORMAT(25H -(H*AV/(CL*H))*(TL-TG)  ,5E11.3,/,13X,6E11.3,//)
C...
C...  ************************************************************(13)
C...
C...     -(KY*AV/(CL*H))*(YS - Y)*(CVV*TG + DHVAP)
         DO 20 I=1,N
         TEP(I)=-P6*(YS(I)-Y(I))*(CVV*TG(I)+DHVAP)
20       CONTINUE
         WRITE(NO,21)(TEP(I),I=1,N,N2)
21       FORMAT(25H -(KY*AV/(CL*H))*(YS-Y)..,5E11.3,/,13X,6E11.3,//)
C...
C...  ************************************************************(14)
C...
C...     TL
C...       T
         WRITE(NO,22)(TLT(I),I=1,N,N2)
22       FORMAT(25H TLT                     ,5E11.3,/,13X,6E11.3,//)
C...
C...  ************************************************************(15)
C...
C...     TL
         WRITE(NO,23)(TL(I),I=1,N,N2)
23       FORMAT(25H TL                      ,5E11.3,/,13X,6E11.3,//)
C...
C...  ****************************************************************
C...
C...  MAP THE JACOBIAN MATRIX OF THE 34 ODES DEFINED IN COMMON/Y/ AND
C...  /F/
      NODE=34
      CALL JMAP(NODE,A,SV,SVOLD,F,FOLD)
C...
C...  COMPUTE THE TEMPORAL EIGENVALUES OF THE 34 ODES DEFINED IN COMMON
C...  /Y/ AND /F/
      CALL EIGEN(NODE,A,WR,WI,Z,RW,IW)
C...
C...  PLOT TG AND TL VS Z
      DO 24 I=1,N
      TGLP(1,I)=TG(I)
      TGLP(2,I)=TL(I)
      TEP(I)=ZL*FLOAT(I-1)/FLOAT(N-1)
24    CONTINUE
      CALL TPLOTS(2,N,TEP,TGLP)
      WRITE(NO,25)T
25    FORMAT(1H ,/,36H 1 - TG(Z,T), 2 - TL(Z,T) VS Z, T = ,F7.3)
C...
C...  STORE THE NUMERICAL SOLUTION FOR PLOTTING (VS T)
100   TLP(IP)=TL(1)
      TGP(IP)=TG(N)
       YP(IP)= Y(N)
       VP(IP)= V
       XP(IP)= X
      ERP(IP)= E
       TP(IP)= T
C...
C...  PLOT THE SOLUTION
      IF(IP.LT.101)RETURN
      CALL TPLOTS(1,IP,TP,TLP)
      WRITE(NO,101)
101   FORMAT(1H ,/,30H TL(0,T) VS T                 )
      CALL TPLOTS(1,IP,TP,TGP)
      WRITE(NO,102)
102   FORMAT(1H ,/,30H TG(8,T) VS T                 )
      CALL TPLOTS(1,IP,TP, YP)
      WRITE(NO,103)
103   FORMAT(1H ,/,30H Y(8,T) VS T                  )
      CALL TPLOTS(1,IP,TP, VP)
      WRITE(NO,104)
104   FORMAT(1H ,/,30H V(T) VS T                    )
      CALL TPLOTS(1,IP,TP, XP)
      WRITE(NO,105)
105   FORMAT(1H ,/,30H X(T) VS T                    )
      CALL TPLOTS(1,IP,TP,ERP)
      WRITE(NO,106)
106   FORMAT(1H ,/,30H E(T) VS T                    )
      RETURN
      END

