*DECK PRC06
      SUBROUTINE INITAL
C...
C...  PI CONTROL OF A DOUBLE EFFECT EVAPORATOR
C...
C...  CONSIDER THE FOLLOWING DOUBLE EFFECT EVAPORATOR
C...
C...                          WS1                    WS2
C...                      ............           ..........+
C...                      .          .           .
C...                   ......        .        ......
C...                  .      .       .       .      .
C...                 .        .      .      .        .
C...                .          .     .     .          .
C...                .          .     .     .          .
C...                .  EVAP 1  .     .     .  EVAP 2  .
C...                .          .     .     .          .
C...                .          .     .     .          .
C...                ............     .     ............
C...                .       M1 .     .     .       M2 .
C...                .          .     .     .          .
C...      WS0 ............  T1 .     ............  T2 .+....... WF
C...                .   .      +...        .   .      .
C...                .    .  A1 .  .        .    .  A2 .
C...                .   .      .  .        .   .      .
C...                .    .  X1 .  .        .    .  X2 .
C...          ...........      .  .  ...........      .
C...          .     .          .  .  .     .          .
C...          +      .        .   .  +      .        .
C...                  .      .    .          .      .
C...                   ......     .           ......
C...                      .       .              .
C...                      .       .              .
C...                      .       .              .
C...                      .       .              .
C...                      .       .              .
C...                      .       .              .
C...               ...    .       .     ...      .
C...          +..... ......       ....... ........
C...               ...                  ...
C...                .                    .
C...               ...                  ...
C...            (PUMP FOR            (PUMP FOR
C...               W1)                  W2)
C...  WHERE
C...
C...     M1,M2       MASS HOLDUP IN EVAPORATORS 1 AND 2, RESPECTIVELY
C...
C...     T1,T2       SOLUTION TEMPERATURES IN EVAPORATORS 1 AND 2
C...
C...     H1,H2       SOLUTION ENTHALPIES IN EVAPORATORS 1 AND 2
C...
C...     X1,X2       SOLUTE MASS FRACTIONS IN EVAPORATORS 1 AND 2
C...
C...     WS1,WS2     STEAM MASS FLOW RATES FROM EVAPORATORS 1 AND 2
C...
C...     W1,W2       SOLUTION MASS FLOW RATES FROM EVAPORATORS 1 AND 2
C...
C...     WF          FEED MASS FLOW RATE TO EVAPORATOR 2
C...
C...     WS0         STEAM MASS FLOW RATE THROUGH THE COIL IN EVAPORATOR
C...                 1 (WS1 IS THE STEAM MASS FLOW RATE THROUGH THE COIL
C...                 IN EVAPORATOR 2)
C...
C...     TS0,TS1     STEAM TEMPERATURES IN THE COILS IN EVAPORATORS 1
C...                 AND 2
C...
C...     HS1,HS2     ENTHALPY OF THE VAPOR LEAVING EVAPORATORS 1 AND 2
C...
C...     H1,H2       ENTHALPY OF THE SOLUTION IN EVAPORATORS 1 AND 2
C...
C...     HVS0,       HEAT OF VAPORIZATION OF STEAM IN THE COILS OF
C...     HVS1        EVAPORATORS 1 AND 2
C...
C...     HF          ENTHALPY OF THE FEED
C...
C...     TF,TS0      TEMPERATURE OF THE FEED AND ENTERING STEAM
C...
C...     XF          MASS FRACTION OF THE FEED
C...
C...     A1,A2       HEAT TRANSFER AREAS OF THE COILS IN EVAPORATOR 1
C...                 AND 2
C...
C...     U1,U2       HEAT TRANSFER COEFFICIENTS FOR THE COILS IN
C...                 EVAPORATORS 1 AND 2
C...
C...     X1SET       SET POINT FOR X1
C...
C...     M1SET,      SET POINTS FOR M1, M2
C...     M2SET
C...
C...     T2SET       SET POINT FOR T2
C...
C...     EX1 = X1 - X1SET
C...
C...     EM1 = M1 - M1SET
C...
C...     EM2 = M2 - M2SET
C...
C...     ET2 = T2 - T2SET
C...
C...     WFSS        WF AT STEADY STATE
C...
C...     W1SS,       W1 AND W2 AT STEADY STATE
C...     W2SS
C...
C...     WS2SS       WS2 AT STEADY STATE
C...
C...     KCF,KC1     GAINS FOR THE FEED FLOW RATE CONTROLLER, M1 AND M2
C...     KC2,KCS     LIQUID LEVEL CONTROLLERS, STEAM FLOW RATE (FROM THE
C...                 SECOND EVAPORATOR) CONTROLLER
C...
C...     TIF,TI1     INTEGRAL TIMES FOR THE FEED FLOW RATE CONTROLLER,
C...     TI2,TIS     M1 AND M2 LIQUID LEVEL CONTROLLERS, STEAM FLOW RATE
C...                 (FROM THE SECOND EVAPORATOR) CONTROLLER
C...
C...  THE MODEL FOR THIS SYSTEM CONSISTS OF THE FOLLOWING EQUATIONS
C...
C...     (A)  TOTAL MASS BALANCE AROUND EACH EVAPORATOR
C...
C...          DM1/DT = W2 - W1 - WS1, M1(0)                       (1)
C...
C...          DM2/DT = WF - W2 - WS2, M2(0)                       (2)
C...
C...     (B)  COMPONENT BALANCE AROUND EACH EVAPORATOR
C...
C...          D(M1*X1)/DT = W2*X2 - W1*X1, M1(0)*X1(0)            (3)
C...
C...          D(M2*X2)/DT = WF*XF - W2*X2, M2(0)*X2(0)            (4)
C...
C...     (C)  ENERGY BALANCE ON THE COIL IN EACH EVAPORATOR
C...          (NEGLECTING VAPOR DYNAMICS)
C...
C...          WS1*HVS1 = U2*A2*(TS1 - T2)                         (5)
C...
C...          WS0*HVS0 = U1*A1*(TS0 - T1)                         (6)
C...
C...     (D)  ENERGY BALANCE FOR EACH EVAPORATOR
C...
C...          D(M1*H1)/DT = W2*H2 - W1*H1 - WS1*HS1 + U1*A1*(TS0 - T1)
C...                                                              (7)
C...                        M1(0)*H1(0)
C...
C...          D(M2*H2)/DT = WF*HF - W2*H2 - WS2*HS2 + U2*A2*(TS1 - T2)
C...                                                              (8)
C...                        M2(0)*H2(0)
C...
C...     (E)  AUXILIARY EQUATIONS
C...
C...          TS1 = T1, TS2 = T2 (ASSUMING NO BOILING         (9)(10)
C...                                  POINT RISE)
C...
C...          LIQUID ENTHALPY
C...
C...             H1 = CP*(T1 - TR), H2 = CP*(T2 - TR)
C...
C...             OR FOR WATER (CP = 1) WITH TR = 0
C...
C...             H1 = T1, H2 = T2                             (11)(12)
C...
C...          VAPOR ENTHALPY
C...
C...             HS1 = 0.42*TS1 + 590.0
C...  
C...             HS2 = 0.42*TS2 + 590.0
C...
C...          HEAT OF VAPORIZATION
C...
C...             HVS0 = 606.6 - 0.57*TS0                        (15)
C...
C...             HVS1 = 606.6 - 0.57*TS1                        (16)
C...
C...     (F)  CONTROLLER EQUATIONS
C...
C...                                          T
C...          WF = WFSS + KCF*(EX1 + (1/TIF)*INT EX1*DT)         (17)
C...                                          0
C...
C...                                          T
C...          W1 = W1SS + KC1*(EM1 + (1/TI1)*INT EM1*DT)         (18)
C...                                          0
C...
C...                                          T
C...          W2 = W2SS + KC2*(EM2 + (1/TI2)*INT EM2*DT)         (19)
C...                                          0
C...
C...                                            T
C...          WS2 = WS2SS + KCS*(ET2 + (1/TIS)*INT ET2*DT)       (20)
C...                                            0
C...
C...          WHERE THE ERRORS ARE
C...
C...             EX1 = X1 - X1SET
C...
C...             EM1 = M1 - M1SET
C...
C...             EM2 = M2 - M2SET
C...
C...             ET2 = T2 - T2SET
C...
C...          NOTE THAT THE FEED FLOW RATE, WF, THE VAPOR RATE FROM THE
C...          SECOND EVAPORATOR, WS2, AND THE MASS HOLDUPS IN THE TWO
C...          EVAPORATORS, M1 AND M2, ARE CONTROLLED TO BRING THE SYSTEM
C...          TO A STEADY STATE CONDITION (BECAUSE OF THE INTEGRAL
C...          ACTION IN THE PI CONTROLLERS).  THE CONTROLLERS ARE NOT
C...          SHOWN IN THE PRECEDING DIAGRAM, BUT THEIR OPERATION IS
C...          CLEAR FROM EQUATIONS (17) TO (20).
C...
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/     M1,     M2,   M1X1,   M2X2,   M1H1,   M2H2,
     2          EX1INT, EM1INT, EM2INT, ET2INT
     3      /F/    M1T,    M2T,  M1X1T,  M2X2T,  M1H1T,  M2H2T,
     4             EX1,    EM1,    EM2,    ET2
     5      /V/     W1,     W2,     WF,    WS0,    WS1,    WS2,
     6      XF,     X1,     X2,     T1,     T2,    TS0,    TS1,
     7     TS2,     TF,     HF,     H1,     H2,    HS1,    HS2,
     8    HVS0,   HVS1,     U1,     U2,     A1,     A2,   WFSS,
     9    W1SS,   W2SS,  WS2SS,    KCF,    KC1,    KC2,    KCS,
     A     TIF,    TI1,    TI2,    TIS,  X1SET,  M1SET,  M2SET,
     B   T2SET
      REAL          M1,     M2,   M1X1,   M2X2,   M1H1,   M2H2,
     1             M1T,    M2T,  M1X1T,  M2X2T,  M1H1T,  M2H2T,
     2             KCF,    KC1,    KC2,    KCS,  M1SET,  M2SET
C...
C...  MODEL CONSTANTS
C...
C...     HEAT TRANSFER AREAS (M**2)
         A1=200.0
         A2=200.0
C...
C...     HEAT TRANSFER COEFFICIENTS (KCAL/SEC-M**2-C)
         U1=5000.0/3600.0
         U2=7000.0/3600.0
C...
C...     ENTERING STEAM TEMPERATURE, FEED TEMPERATURE (K)
         TS0=164.0
          TF= 38.0
C...
C...     FEED COMPOSITION (KG SOLUTE/KG SOLUTION)
         XF=0.02
C...
C...     STEADY STATE FLOWS (KG SOLUTION/SEC)
         WFSS=40.0
         W1SS=30.0
         W2SS=30.0
         WS2SS=10.0
C...
C...     CONTROLLER SET POINTS
C...        (KG SOLUTE/KG SOLUTION)
            X1SET=0.250
C...
C...        (KG SOLUTION)
            M1SET=250.0
            M2SET=250.0
C...
C...        (K)
            T2SET=38.0
C...
C...     CONTROLLER GAINS (KG/SEC-KG SOLUTE/KG SOLUTION,
C...        (KG SOLUTION/SEC-KG SOLUTE/KG SOLUTION)
            KCF=100.0
C...
C...        (KG SOLUTION/SEC-KG SOLUTION)
            KC1=1.0
            KC2=1.0
C...
C...        (KG STEAM/SEC-K)
            KCS=1.0
C...
C...     CONTROLLER INTEGRAL TIMES (SEC)
         TIF=250.0
         TI1=250.0
         TI2=250.0
         TIS=250.0
C...
C...  INITIAL CONDITIONS
C...
C...     (KG SOLUTION)
         M1=250.0
         M2=250.0
C...
C...     (KG SOLUTE)
         M1X1=250.0*0.1
         M2X2=250.0*0.1
C...
C...     (KCAL)
         M1H1=250.0*35.0
         M2H2=250.0*35.0
C...
C...     INTEGRATED ERROR
C...        (KG SOLUTE-SEC/KG SOLUTION)
            EX1INT=0.
C...
C...        (KG SOLUTION-SEC)
            EM1INT=0.
            EM2INT=0.
C...
C...        (K-SEC)
            ET2INT=0.
C...
C...  COMPUTE THE INITIAL DERIVATIVES
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/     M1,     M2,   M1X1,   M2X2,   M1H1,   M2H2,
     2          EX1INT, EM1INT, EM2INT, ET2INT
     3      /F/    M1T,    M2T,  M1X1T,  M2X2T,  M1H1T,  M2H2T,
     4             EX1,    EM1,    EM2,    ET2
     5      /V/     W1,     W2,     WF,    WS0,    WS1,    WS2,
     6      XF,     X1,     X2,     T1,     T2,    TS0,    TS1,
     7     TS2,     TF,     HF,     H1,     H2,    HS1,    HS2,
     8    HVS0,   HVS1,     U1,     U2,     A1,     A2,   WFSS,
     9    W1SS,   W2SS,  WS2SS,    KCF,    KC1,    KC2,    KCS,
     A     TIF,    TI1,    TI2,    TIS,  X1SET,  M1SET,  M2SET,
     B   T2SET
      REAL          M1,     M2,   M1X1,   M2X2,   M1H1,   M2H2,
     1             M1T,    M2T,  M1X1T,  M2X2T,  M1H1T,  M2H2T,
     2             KCF,    KC1,    KC2,    KCS,  M1SET,  M2SET
C...
C...  MODEL ALGEBRA
C...
C...     MASS FRACTION IN EACH EVAPORATOR
         X1=M1X1/M1
         X2=M2X2/M2
C...
C...     THE MASS FRACTION IN EACH EVAPORATOR MUST BE NONNEGATIVE
         IF(X1.LT.0.)X1=0.
         IF(X2.LT.0.)X2=0.
C...
C...     LIQUID TEMPERATURE IN EACH EVAPORATOR
         H1=M1H1/M1
         H2=M2H2/M2
         T1=H1
         T2=H2
C...
C...     TEMPERATURE OF VAPOR FROM EACH EVAPORATOR
         TS1=T1
         TS2=T2
C...
C...     HEAT OF VAPORIZATION OF THE STEAM IN EACH EVAPORATOR COIL
         HVS0=606.6-0.57*TS0
         HVS1=606.6-0.57*TS1
C...
C...     VAPOR FLOW IN THE COIL OF EACH EVAPORATOR
         WS0=(U1*A1/HVS0)*(TS0-T1)
         WS1=(U2*A2/HVS1)*(TS1-T2)
C...
C...     CONTROLLER ERRORS
         EX1=X1-X1SET
         EM1=M1-M1SET
         EM2=M2-M2SET
         ET2=T2-T2SET
C...
C...     CONTROLLER EQUATIONS
         WF=WFSS+KCF*(EX1+(1./TIF)*EX1INT)
         W1=W1SS+KC1*(EM1+(1./TI1)*EM1INT)
         W2=W2SS+KC2*(EM2+(1./TI2)*EM2INT)
         WS2=WS2SS+KCS*(ET2+(1./TIS)*ET2INT)
C...
C...     THE MANIPULATED VARIABLES MUST BE NONNEGATIVE
         IF(WF.LT.0.)WF=0.
         IF(W1.LT.0.)W1=0.
         IF(W2.LT.0.)W2=0.
         IF(WS2.LT.0.)WS2=0.
C...
C...     VAPOR ENTHALPIES
         HS1=0.42*TS1+590.0
         HS2=0.42*TS2+590.0
C...
C...     FEED ENTHALPY
         HF=TF
C...
C...  MODEL DIFFERENTIAL EQUATIONS
C...
C...     TOTAL MASS BALANCES
         M1T=W2-W1-WS1
         M2T=WF-W2-WS2
C...
C...     COMPONENT MASS BALANCES
         M1X1T=W2*X2-W1*X1
         M2X2T=WF*XF-W2*X2
C...
C...     ENERGY BALANCES
         M1H1T=W2*H2-W1*H1-WS1*HS1+U1*A1*(TS0-T1)
         M2H2T=WF*HF-W2*H2-WS2*HS2+U2*A2*(TS1-T2)
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/     M1,     M2,   M1X1,   M2X2,   M1H1,   M2H2,
     2          EX1INT, EM1INT, EM2INT, ET2INT
     3      /F/    M1T,    M2T,  M1X1T,  M2X2T,  M1H1T,  M2H2T,
     4             EX1,    EM1,    EM2,    ET2
     5      /V/     W1,     W2,     WF,    WS0,    WS1,    WS2,
     6      XF,     X1,     X2,     T1,     T2,    TS0,    TS1,
     7     TS2,     TF,     HF,     H1,     H2,    HS1,    HS2,
     8    HVS0,   HVS1,     U1,     U2,     A1,     A2,   WFSS,
     9    W1SS,   W2SS,  WS2SS,    KCF,    KC1,    KC2,    KCS,
     A     TIF,    TI1,    TI2,    TIS,  X1SET,  M1SET,  M2SET,
     B   T2SET
      REAL          M1,     M2,   M1X1,   M2X2,   M1H1,   M2H2,
     1             M1T,    M2T,  M1X1T,  M2X2T,  M1H1T,  M2H2T,
     2             KCF,    KC1,    KC2,    KCS,  M1SET,  M2SET
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING THE SOLUTION
      DIMENSION   TMP(51),  XP(2,51),  TP(2,51),  MP(2,51),
     1            WFP(51),  WP(2,51), WSP(3,51)
      REAL MP
      DATA IP/0/
C...
C...  PRINT A HEADING FOR THE NUMERICAL SOLUTION
      IF((IP/5*5).NE.IP)GO TO 3
      WRITE(NO,1)
1     FORMAT(1H ,/,
     1 5X,1HT,8X,2HM1,8X,2HM2,8X,2HX1,8X,2HX2,8X,2HT1,8X,2HT2,
     2      /,14X,2HWF,8X,2HW1,8X,2HW2,7X,3HWS0,7X,3HWS1,7X,3HWS2,
     3      /,13X,3HM1T,7X,3HM2T,5X,5HM1X1T,5X,5HM2X2T,5X,5HM1H1T,
     4       5X,5HM2H2T)
C...
C...  PRINT THE SOLUTION
      WRITE(NO,2)     T,    M1,    M2,    X1,    X2,    T1,    T2,
     1               WF,    W1,    W2,   WS0,   WS1,   WS2,
     2              M1T,   M2T, M1X1T, M2X2T, M1H1T, M2H2T
2     FORMAT(F6.0,2F10.1,2F10.4,2F10.1,/,
     1                       6X,6F10.3,/,
     2                       6X,4F10.4,2F10.1)
3     IP=IP+1
C...
C...  STORE THE SOLUTION FOR PLOTTING
      TMP(IP)=T
      XP(1,IP)=X1
      XP(2,IP)=X2
      TP(1,IP)=T1
      TP(2,IP)=T2
      MP(1,IP)=M1
      MP(2,IP)=M2
      WFP(IP)=WF
      WP(1,IP)=W1
      WP(2,IP)=W2
      WSP(1,IP)=WS0
      WSP(2,IP)=WS1
      WSP(3,IP)=WS2
      IF(IP.LT.51)RETURN
C...
C...  PLOT THE SOLUTION
      CALL SPLOTS(2,IP,TMP, XP)
      WRITE(NO,10)
10    FORMAT(1H ,//,35H X1, X2 VS T                       )
      CALL SPLOTS(2,IP,TMP, TP)
      WRITE(NO,11)
11    FORMAT(1H ,//,35H T1, T2 VS T                       )
      CALL SPLOTS(2,IP,TMP, MP)
      WRITE(NO,12)
12    FORMAT(1H ,//,35H M1, M2 VS T                       )
      CALL SPLOTS(1,IP,TMP,WFP)
      WRITE(NO,13)
13    FORMAT(1H ,//,35H WF VS T                           )
      CALL SPLOTS(2,IP,TMP, WP)
      WRITE(NO,14)
14    FORMAT(1H ,//,35H W1, W2 VS T                       )
      CALL SPLOTS(3,IP,TMP,WSP)
      WRITE(NO,15)
15    FORMAT(1H ,//,35H WS0, WS1, WS2 VS T                )
      IP=0
      RETURN
      END
