*DECK SEP20
      SUBROUTINE INITAL
C...
C...  COUNTERCURRENT LIQUID-LIQUID EXTRACTOR FOR A PARTIALLY MISCIBLE
C...  SYSTEM
C...
C...  ACETONE IS TO BE RECOVERED FROM AN ACETONE-METHYLISOBUTYLKETONE
C...  (MIK) SOLUTION IN A SIX-STAGE COUNTERCURRENT LIQUID-LIQUID EX-
C...  TRACTOR USING WATER (THAT ALSO CONTAINS SMALL AMOUNTS OF ACETONE
C...  AND MIK) AS A SOLVENT.  THE FEED ENTERS THE FIRST STAGE AT A MASS
C...  FLOW RATE FL0 AND ANALYZES W1L0 AND W2L0 MASS FRACTIONS OF WATER
C...  AND ACETONE, RESPECTIVELY.  THE SOLVENT ENTERS THE SIXTH STAGE AT
C...  A MASS FLOW RATE FH7 AND ANALYZES W1H7 AND W2H7 MASS FRACTIONS OF
C...  WATER AND ACETONE, RESPECTIVELY.
C...
C...  INITIALLY, THE CONCENTRATION OF EACH PHASE IN EACH STAGE CORRES-
C...  PONDS TO THAT OF THE FEED AND THE SOLVENT.  IN THIS SYSTEM, THE
C...  WATER AND MIK ARE PARTIALLY MISCIBLE.  THE MASS FLOW RATES OF THE
C...  TWO PHASES EXITING EACH STAGE ARE GIVEN BY
C...
C...     FHN = K1*MHN, FLN = K2*(MHN + MLN)                    (1)(2)
C...
C...  WHERE MHN AND MLN ARE THE MASSES OF THE HEAVY AND LIGHT PHASES
C...  IN STAGE N.
C...
C...  ASSUME THAT EQUILIBRIUM IS ATTAINED IN EACH STAGE WITH THE
C...  FOLLOWING EQUILBRIUM RELATIONSHIPS
C...
C...     W2H = A + B*W1L + C*W1L**2                               (3)
C...
C...     W1H = M + N*W2H + P*W2H**2                               (4)
C...
C...     W2L = D + E*W1L + F*W1L**2                               (5)
C...
C...  THE DETAILS OF THE MATHEMATICAL MODEL FOLLOW.
C...
C...  THE VARIABLES FOR STAGE N ARE INDICATED IN THE FOLLOWING DIAGRAM
C...  STAGE N ARE INDICATED IN THE FOLLOWING DIAGRAM
C...
C...                      ...............
C...             FLN-1    .  W1LN W2LN  .     FLN
C...          ............+     W3LN    .............+
C...                      .      LN     .
C...                      ...............
C...              FHN     .  W1HN W2HN  .    FHN+1
C...          +............     W3HN    +.............
C...                      .      HN     .
C...                      ...............
C...
C...  A TOTAL MATERIAL BALANCE FOR STAGE N GIVES
C...
C...     D(LN + HN)/DT = FLN-1 + FHN+1 - FLN - FHN                (6)
C...
C...  BALANCES FOR COMPONENTS 1 (WATER) AND 2 (ACETONE) GIVE
C...                                                              (7)
C...     D(LN*W1LN + HN*W1HN)/DT = FLN-1*W1LN-1 + FHN+1*W1HN+1
C...
C...                             - FLN*W1LN - FHN*W1HN
C...
C...     D(LN*W2LN + HN*W2HN)/DT = FLN-1*W2LN-1 + FHN+1*W2HN+1
C...                                                              (8)
C...                             - FLN*W2LN - FHN*W2HN
C...
C...  IF WE SET
C...
C...     MN = LN + HN                                             (9)
C...
C...     MN*W1MN = LN*W1LN + HN*W1HN                             (10)
C...
C...     MN*W2MN = LN*W2LN + HN*W2HN                             (11)
C...
C...  EQUATIONS (9), (10) AND (11) CAN THEN BE SUBSTITUTED IN THE
C...  DERIVATIVES OF EQUATIONS (6), (7) AND (8) TO DEFINE THREE
C...  DEPENDENT VARIABLES OF THE ODES, MN, MN*W1MN AND MN*W2MN.  THE
C...  VARIABLES IN THESE EQUATIONS ARE
C...
C...     W1LN     MASS FRACTION OF COMPONENT 1   (WATER) IN THE LIGHT
C...              (MIK) PHASE IN STAGE N
C...
C...     W2LN     MASS FRACTION OF COMPONENT 2 (ACETONE) IN THE LIGHT
C...              (MIK) PHASE IN STAGE N
C...
C...     W3LN     MASS FRACTION OF COMPONENT 3     (MIK) IN THE LIGHT
C...              (MIK) PHASE IN STAGE N
C...
C...     W1HN     MASS FRACTION OF COMPONENT 1   (WATER) IN THE HEAVY
C...              (WATER) PHASE IN STAGE N
C...
C...     W2HN     MASS FRACTION OF COMPONENT 2 (ACETONE) IN THE HEAVY
C...              (WATER) PHASE IN STAGE N
C...
C...     W3HN     MASS FRACTION OF COMPONENT 3     (MIK) IN THE HEAVY
C...              (WATER) PHASE IN STAGE N
C...
C...     W1MN     MASS FRACTION OF COMPONENT 1   (WATER) IN THE COMBINED
C...              (MIK + WATER) PHASES IN STAGE N
C...
C...     W2MN     MASS FRACTION OF COMPONENT 2 (ACETONE) IN THE COMBINED
C...              (MIK + WATER) PHASES IN STAGE N
C...
C...     W3MN     MASS FRACTION OF COMPONENT 3     (MIK) IN THE COMBINED
C...              (MIK + WATER) PHASES IN STAGE N
C...
C...     FLN      MASS FLOW RATE OF THE LIGHT   (MIK) PHASE FROM STAGE N
C...
C...     FHN      MASS FLOW RATE OF THE HEAVY (WATER) PHASE FROM STAGE N
C...
C...     LN       HOLDUP OF THE LIGHT   (MIK) PHASE IN STAGE N
C...
C...     HN       HOLDUP OF THE HEAVY (WATER) PHASE IN STAGE N
C...
C...     MN       COMBINED HOLDUP OF THE LIGHT AND HEAVY PHASES (MIK +
C...              WATER) PHASES IN STAGE N
C...
C...  WE WILL ALSO NEED AN EQUATION DERIVED BY ELIMINATING MN, HN AND
C...  LN FROM EQUATIONS (9), (10) AND (11)
C...
C...     MN*W2MN = HN*W2HN + (MN - HN)*W2LN
C...
C...     MN*W1MN = HN*W1HN + (MN - HN)*W1LN
C...
C...  OR
C...
C...     (W2MN - W2LN)*MN = (W2HN - W2LN)*HN
C...
C...     (W1MN - W1LN)*MN = (W1HN - W1LN)*HN
C...
C...  DIVISION OF ONE EQUATION BY THE OTHER GIVES
C...
C...     (W2MN - W2LN)/(W1MN - W1LN) = (W2HN - W2LN)/(W1HN - W1LN)
C...
C...  OR SOLVING FOR W2MN
C...
C...     W2MN = W2LN + ((W2LN - W2HN)/(W1LN - W1HN))*(W1MN - W1LN)  (12)
C...
C...  THE POINTS (W1LN,W2LN), (W1HN,W2HN) AND (W1MN,W2MN) LIE ON A
C...  STRAIGHT LINE.  FURTHERMORE, THE FIRST TWO POINTS WILL BE ON THE
C...  EQUIULIBRIUM CURVES FOR EACH PHASE.  THUS, IF WE KNOW ONE POINT,
C...  THE OTHER TWO CAN BE CALCULATED (FROM THE LINEAR EQUATIONS AND
C...  THE EQUILIBRIUM RELATIONSHIPS).
C...
C...  IF WE ASSUME THAT THE POINT (W1MN,W2MN) HAS BEEN DETERMINED BY
C...  INTEGRATING THE ODES NUMERICALLY, IN ORDER TO COMPUTE THE OTHER
C...  TWO POINTS WE NEED THE EQUILIBRIUM RELATIONSHIPS.  AN ALGORITHM
C...  FOR THE CALCULATIONS FOLLOWS
C...
C...     (1)  ASSUME W2LN = 0.5*W2MN
C...
C...     (2)  CALCULATE W1LN (USING THE EQUILIBRIUM RELATIONSHIP FOR
C...                          THE LIGHT PHASE, EQUATION (15))
C...
C...     (3)  CALCULATE W1HN (USING THE EQUILIBRIUM RELATIONSHIP FOR
C...                          THE CONJUGATE LINE, EQUATION (17))
C...
C...     (4)  CALCULATE W2HN (USING THE EQUILIBRIUM RELATIONSHIP FOR
C...                          THE HEAVY PHASE, EQUATION (16))
C...
C...                        (C)
C...     (5)  CALCULATE W2MN   (USING EQUATION (12))
C...
C...                     (C)
C...     (6)  IF ABS(W2MN    - W2MN) LE 0.001, EXIT THE ALGORITHM
C...             (SINCE THE CALCULATION HAS CONVERGED)
C...
C...                     (C)
C...          IF ABS(W2MN    - W2MN) GT 0.001
C...                                        (C)
C...             SET W2LN = W2LN - 0.5*(W2MN    - W2MN) AND GO TO (2)
C...
C...          THE LAST STEP IN (6) FOR THE CALCULATION OF A NEW W2LN IS
C...          BASED ON INTERVAL HALVING.
C...
C...  ONCE THE CALCULATION HAS CONVERGED, WE HAVE THE THREE COMPOSITION
C...  POINTS AND THEREFORE HN AND LN CAN BE COMPUTED FROM
C...
C...     LN = (W1MN - W1HN)/(W1LN - W1HN)*MN                      (13)
C...
C...     HN = MN - LN                                             (14)
C...
C...  THE EQUILIBRIUM RELATIONSHIPS FOR THE WATER-ACETONE-MIK SYSTEM
C...  ARE
C...
C...  MIK PHASE                                                   (15)
C...
C...     W2LN =  0.38658 + 0.66147*W1LN - 1.255*W1LN**2,  W2LN GT 0.4345
C...
C...     W2LN = -0.22626 + 12.852*W1LN - 60.356*W1LN**2,  W2LN LE 0.4345
C...
C...  WATER PHASE                                                 (16)
C...
C...     W2HN = 0.67714 - 0.37644*W1HN - 0.3233*W1HN**2,  W2HN LE 0.35
C...
C...     W2HN = 0.38658 + 0.66147*W1HN - 1.255*W1HN**2,   W2HN GT 0.35
C...
C...  CONJUGATE
C...
C...     W2LN = 0.16279 + 1.4376*W1HN - 1.6344*W1HN**2            (17)
C...
C...  THE MODEL PARAMETERS ARE
C...
C...     FEED COMPOSITION AND FLOW
C...
C...        W1L0 = 0, W2L0 = 0.4, W3L0 = 0.6, FL0 = 5
C...
C...     SOLVENT COMPOSITION AND FLOW
C...
C...        W1H7 = 0.98, W2H7 = 0.02, W3H7 = 0., FH7 = 4
C...
C...     INITIAL CONDITIONS
C...
C...        W1LN(0) = 0, W2LN(0) = 0.4, W1HN(0) = 0.98, W2HN(0) = 0.02
C...
C...        LN(0) = HN(0) = 500, N = 1, 2, ..., 6
C...
C...     FLOWS
C...
C...        FHN = 0.01*HN, FLN = 0.005*(HN + LN), N = 1, 2, ..., 6
C...
C...  SIMULATE THIS SYSTEM TO DETERMINE
C...
C...     (1)  CONCENTRATIONS OF STREAMS EXITING FROM THE EXTRACTION
C...          UNIT
C...
C...     (2)  PERCENT RECOVERY OF ACETONE
C...
C...     (3)  FLOW RATES OUT OF EACH STAGE
C...
C...  COMMON AREA
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      M(6),   MW1M(6),   MW2M(6)
     2      /F/   DMDT(6), DMW1DT(6), DMW2DT(6)
     3      /V/   FL(0:6),   FH(1:7),
     4           W1L(0:6),  W2L(0:6),
     5           W1H(1:7),  W2H(1:7),
     6           W1M(1:6),  W2M(1:6),
     7             L(1:6),    H(1:6),
     8                 NS,        IP
      REAL              M,      MW1M,      MW2M,         L
C...
C...
C...  NUMBER OF STAGES
      NS=6
C...
C...  FEED COMPOSITION AND FLOW
      W1L(0)=0.
      W2L(0)=0.4
      FL(0)=5.
C...
C...  SOLVENT COMPOSITION AND FLOW
      W1H(NS+1)=0.98
      W2H(NS+1)=0.02
      FH(NS+1)=4.
C...
C...  INITIAL CONDITIONS
C...
C...  LIGHT PHASE
      DO 1 N=1,NS
      L(N)=500.
      W1L(N)=0.
      W2L(N)=0.4
C...
C...  HEAVY PHASE
      H(N)=500.
      W1H(N)=0.98
      W2H(N)=0.02
C...
C...  MIXED PHASE
      MW1M(N)=L(N)*W1L(N)+H(N)*W1H(N)
      MW2M(N)=L(N)*W2L(N)+H(N)*W2H(N)
      M(N)=L(N)+H(N)
1     CONTINUE
C...
C...  COMPUTE THE INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      M(6),   MW1M(6),   MW2M(6)
     2      /F/   DMDT(6), DMW1DT(6), DMW2DT(6)
     3      /V/   FL(0:6),   FH(1:7),
     4           W1L(0:6),  W2L(0:6),
     5           W1H(1:7),  W2H(1:7),
     6           W1M(1:6),  W2M(1:6),
     7             L(1:6),    H(1:6),
     8                 NS,        IP
      REAL              M,      MW1M,      MW2M,         L
C...
C...  MODEL ALGEBRA
C...
C...     MIXED PHASE COMPOSITIONS
         DO 2 N=1,NS
         W1M(N)=MW1M(N)/M(N)
         W2M(N)=MW2M(N)/M(N)
2        CONTINUE
C...
C...     EQUILIBRIUM LIGHT AND HEAVY PHASE COMPOSITIONS
         CALL EQUIL
C...
C...     LIGHT PHASE HOLDUP
         CALL LN
C...
C...     HEAVY PHASE HOLDUP
         CALL HN
C...
C...     LIGHT PHASE FLOW RATE
         CALL FLN
C...
C...     HEAVY PHASE FLOW RATE
         CALL FHN
C...
C...  MODEL DIFFERENTIAL EQUATIONS
      DO 1 N=1,NS
      DMDT(N)  =FL(N-1)         +FH(N+1)
     1         -FL(N  )         -FH(N  )
      DMW1DT(N)=FL(N-1)*W1L(N-1)+FH(N+1)*W1H(N+1)
     1         -FL(N  )*W1L(N  )-FH(N  )*W1H(N  )
      DMW2DT(N)=FL(N-1)*W2L(N-1)+FH(N+1)*W2H(N+1)
     1         -FL(N  )*W2L(N  )-FH(N  )*W2H(N  )
1     CONTINUE
      RETURN
      END
      SUBROUTINE EQUIL
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      M(6),   MW1M(6),   MW2M(6)
     2      /F/   DMDT(6), DMW1DT(6), DMW2DT(6)
     3      /V/   FL(0:6),   FH(1:7),
     4           W1L(0:6),  W2L(0:6),
     5           W1H(1:7),  W2H(1:7),
     6           W1M(1:6),  W2M(1:6),
     7             L(1:6),    H(1:6),
     8                 NS,        IP
      REAL              M,      MW1M,      MW2M,         L
      COMMON/IO/       NI,        NO
C...
C...  SUBROUTINE EQUIL COMPUTES THE LIGHT AND HEAVY PHASE EQUILIBRIUM
C...  COMPOSITIONS
C...
C...  SET THE MAXIMUM NUMBER OF ITERATIONS AND TOLERANCE IN THE CAL-
C...  CULATION OF THE LIGHT AND HEAVY PHASE COMPOSITIONS
      NMAX=25
      TOL=0.001
C...
C...  STEP THROUGH THE NS STAGES
      DO 1 N=1,NS
C...
C...  INITIALIZE THE ITERATION COUNTER
      ITER=0
C...
C...  INITIAL ESTIMATE OF W2L(N)
      W2L(N)=0.5*W2M(N)
C...
C...  W1L(N) FROM THE LIGHT PHASE EQUILIBRIUM RELATIONSHIP
2     W1L(N)=ELIGHT(W2L(N))
C...
C...  W1H(N) FROM THE CONJUGATE LINE EQUILIBRIUM RELATIONSHIP
      W1H(N)=ECONJG(W2L(N))
C...
C...  W2H(N) FROM THE HEAVY PHASE EQUILIBRIUM RELATIONSHIP
      W2H(N)=EHEAVY(W1H(N))
C...
C...        (C)
C...  W2M(N)    FROM THE DEFINITIONS OF M(N), W1M(N) AND W2M(N)
      W2MC=W2L(N)+(W2L(N)-W2H(N))/(W1L(N)-W1H(N))*(W1M(N)-W1L(N))
C...
C...  CHECK FOR CONVERGENCE
      IF(ABS(W2MC-W2M(N)).LE.0.001)THEN
C...
C...     CONVERGENCE HAS BEEN ACHIEVED AND THE LIGHT AND HEAVY PHASE
C...     COMPOSITIONS ARE AVAILABLE FOR SUBSEQUENT CALCULATIONS
         GO TO 1
C...
C...     CONVERGENCE HAS NOT BEEN ACHIEVED, SO CONTINUE THE ITERATIVE
C...     CALCULATION OF THE LIGHT AND HEAVY PHASE COMPOSITIONS
      ELSE
         W2L(N)=W2L(N)-0.5*(W2MC-W2M(N))
         ITER=ITER+1
         IF(ITER.LT.NMAX)GO TO 2
         WRITE(NO,3)N
3        FORMAT(1H ,//,30H CONVERENCE FAILURE FOR STAGE ,I2)
         GO TO 1
      END IF
C...
C...  GO ON TO THE NEXT STAGE
1     CONTINUE
      RETURN
      END
      FUNCTION ELIGHT(W2LN)
C...
C...  FUNCTION ELIGHT COMPUTES W1LN FROM W2LN
      IF(W2LN.GT.0.4345)THEN
         A=-1.255
         B=0.66147
         C=0.38658-W2LN
         ELIGHT=(-B+SQRT(B**2-4.0*A*C))/(2.0*A)
         IF(ELIGHT.LT.0.)ELIGHT=0.
         IF(ELIGHT.GT.1.)ELIGHT=1.
         RETURN
      ELSE
         A=-60.356
         B= 12.852
         C=-0.22626-W2LN
         ELIGHT=(-B+SQRT(B**2-4.0*A*C))/(2.0*A)
         IF(ELIGHT.LT.0.)ELIGHT=0.
         IF(ELIGHT.GT.1.)ELIGHT=1.
         RETURN
      END IF
      END
      FUNCTION ECONJG(W2LN)
C...
C...  FUNCTION ECONJG COMPUTES W1HN FROM W2LN
         A=-1.6344
         B= 1.4367
         C=0.16279-W2LN
         ECONJG=(-B-SQRT(B**2-4.0*A*C))/(2.0*A)
         IF(ECONJG.LT.0.)ECONJG=0.
         IF(ECONJG.GT.1.)ECONJG=1.
         RETURN
      END
      FUNCTION EHEAVY(W1HN)
C...
C...  FUNCTION EHEAVY COMPUTES W2HN FROM W1HN
      COMMON/IO/NI,NO
      W2HNA=0.67714-0.37644*W1HN-0.3233*W1HN**2
      W2HNB=0.38658+0.66147*W1HN-1.255 *W1HN**2
      IF(W2HNA.LE.0.35)THEN
         EHEAVY=W2HNA
         GO TO 1
      ELSE IF(W2HNB.GT.0.35)THEN
         EHEAVY=W2HNB
         GO TO 1
      ELSE
         WRITE(NO,2)
2        FORMAT(1H ,//,
     1   51H WH2N NOT CALCULATED CORRECTLY FROM THE HEAVY PHASE,/,
     2   51H EQUILIBRIUM EQUATIONS SO RUN IS TERMINATED        )
         STOP
      END IF
C...
C...  CONSTRAIN THE MASS FRACTIONS TO MEANINGFUL VALUES
1     IF(EHEAVY.LT.0.)EHEAVY=0.
      IF(EHEAVY.GT.1.)EHEAVY=1.
      RETURN
      END
      SUBROUTINE LN
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      M(6),   MW1M(6),   MW2M(6)
     2      /F/   DMDT(6), DMW1DT(6), DMW2DT(6)
     3      /V/   FL(0:6),   FH(1:7),
     4           W1L(0:6),  W2L(0:6),
     5           W1H(1:7),  W2H(1:7),
     6           W1M(1:6),  W2M(1:6),
     7             L(1:6),    H(1:6),
     8                 NS,        IP
      REAL              M,      MW1M,      MW2M,         L
C...
C...  SUBROUTINE LN COMPUTES THE LIGHT PHASE HOLDUP, L(N)
      DO 1 N=1,NS
      L(N)=(W1M(N)-W1H(N))/(W1L(N)-W1H(N))*M(N)
1     CONTINUE
      RETURN
      END
      SUBROUTINE HN
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      M(6),   MW1M(6),   MW2M(6)
     2      /F/   DMDT(6), DMW1DT(6), DMW2DT(6)
     3      /V/   FL(0:6),   FH(1:7),
     4           W1L(0:6),  W2L(0:6),
     5           W1H(1:7),  W2H(1:7),
     6           W1M(1:6),  W2M(1:6),
     7             L(1:6),    H(1:6),
     8                 NS,        IP
      REAL              M,      MW1M,      MW2M,         L
C...
C...  SUBROUTINE HN COMPUTES THE HEAVY PHASE HOLDUP, H(N)
      DO 1 N=1,NS
      H(N)=M(N)-L(N)
1     CONTINUE
      RETURN
      END
      SUBROUTINE FLN
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      M(6),   MW1M(6),   MW2M(6)
     2      /F/   DMDT(6), DMW1DT(6), DMW2DT(6)
     3      /V/   FL(0:6),   FH(1:7),
     4           W1L(0:6),  W2L(0:6),
     5           W1H(1:7),  W2H(1:7),
     6           W1M(1:6),  W2M(1:6),
     7             L(1:6),    H(1:6),
     8                 NS,        IP
      REAL              M,      MW1M,      MW2M,         L
C...
C...  SUBROUTINE FLN COMPUTES THE LIGHT PHASE FLOW RATE, FL(N)
      DO 1 N=1,NS
      FL(N)=0.005*(H(N)+L(N))
1     CONTINUE
      RETURN
      END
      SUBROUTINE FHN
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      M(6),   MW1M(6),   MW2M(6)
     2      /F/   DMDT(6), DMW1DT(6), DMW2DT(6)
     3      /V/   FL(0:6),   FH(1:7),
     4           W1L(0:6),  W2L(0:6),
     5           W1H(1:7),  W2H(1:7),
     6           W1M(1:6),  W2M(1:6),
     7             L(1:6),    H(1:6),
     8                 NS,        IP
      REAL              M,      MW1M,      MW2M,         L
C...
C...  SUBROUTINE FHN COMPUTES THE HEAVY PHASE FLOW RATE, FH(N)
      DO 1 N=1,NS
      FH(N)=0.01*H(N)
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      M(6),   MW1M(6),   MW2M(6)
     2      /F/   DMDT(6), DMW1DT(6), DMW2DT(6)
     3      /V/   FL(0:6),   FH(1:7),
     4           W1L(0:6),  W2L(0:6),
     5           W1H(1:7),  W2H(1:7),
     6           W1M(1:6),  W2M(1:6),
     7             L(1:6),    H(1:6),
     8                 NS,        IP
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION    TP(51),W2HP(2,51),    RP(51)
      REAL              M,      MW1M,      MW2M,         L
C...
C...  PRINT THE NUMERICAL SOLUTION EVERY TENTH CALL TO SUBROUTINE PRINT
      IP=IP+1
      IF(((IP-1)/10*10).NE.(IP-1))GO TO 2
      WRITE(NO,1)T,    W1L(0),   W1L(NS),    W1H(1), W1H(NS+1),
     1                 W2L(0),   W2L(NS),    W2H(1), W2H(NS+1),
     2                  FL(0),    FL(NS),     FH(1),  FH(NS+1),
     3                   L(1),     L(NS),      H(1),     H(NS)
1     FORMAT(1H ,//,9X,1HT,
     1       4X,6HW1L(0),3X,7HW1L(NS),4X,6HW1H(1),1X,9HW1H(NS+1),/,
     2      14X,6HW2L(0),3X,7HW2L(NS),4X,6HW2H(1),1X,9HW2H(NS+1),/,
     3      14X,6H FL(0),3X,7H FL(NS),4X,6H FH(1),1X,9H FH(NS+1),/,
     4      14X,6H  L(1),3X,7H  L(NS),4X,6H  H(1),1X,9H    H(NS),/,
     5 F10.3,4F10.4,/,10X,4F10.4,/,10X,4F10.4,/,10X,4F10.4)
C...
C...  STORE THE SOLUTION FOR PLOTTING
2     TP(IP)=T
      W2HP(1,IP)=W2H(NS+1)
      W2HP(2,IP)=W2H(   1)
      RP(IP)=(W2H(1)*FH(1))/(W2L(0)*FL(0)+W2H(NS+1)*FH(NS+1))*100.
C...
C...  PLOT THE SOLUTION
      IF(IP.LT.51)RETURN
      CALL SPLOTS(2,IP,TP,W2HP)
      WRITE(NO,10)
10    FORMAT(1H ,//,43H 1 - W2H(NS+1), 2 - W2H(1) VS T            )
      CALL SPLOTS(1,IP,TP,  RP)
      WRITE(NO,11)
11    FORMAT(1H ,//,43H PERCENT RECOVERY OF ACETONE VS T          )
      RETURN
      END
