      SUBROUTINE INITAL
C...
C...  MA1 - NH3 ABSORBER
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION LS, KYA
C...
C...  COMMON AREA TO LINK SUBROUTINES
      PARAMETER(NZ=21)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/  XA(NZ),  YA(NZ)
     +      /F/ XAT(NZ), YAT(NZ)
     +      /S/ XAZ(NZ), YAZ(NZ), YAS(NZ)
     +      /C/     XA1,     YA1,     XA2,     YA2,
     +      ZL,      GS,      LS,     KYA,       S
C...
C...  COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE
      COMMON/E/A(2)
C...
C...  INPUT/OUTPUT (I/O) UNIT NUMBERS
      COMMON/IO/ NI, NO
C...
C...  GIVEN END COMPOSITIONS (MOLE FRACTIONS)
C...
C...     ENTERING GAS
         YA1=0.0825D0
C...
C...     ENTERING LIQUID
         XA2=0.0D0
C...
C...     EXITING GAS
         YA2=0.003D0
C...
C...  CONVERT PRECEDING MOL FRACTIONS TO MOL RATIOS (BASED ON AIR AND
C...  WATER AS THE NONTRANSFERRED COMPONENTS)
C...
C...     ENTERING GAS
         YA1=YA1/(1.0D0-YA1)
C...
C...     ENTERING LIQUID
         XA2=XA2/(1.0D0-XA2)
C...
C...     EXITING GAS
         YA2=YA2/(1.0D0-YA2)
C...
C...  FLOW RATES (MOLS/M**2-MIN, SO THE COLUMN CROSS SECTIONAL AREA
C...  IS REQUIRED)
C...
C...     PI
         PI=4.0D0*DATAN(1.0D0)
C...
C...     COLUMN DIAMETER (M)
         D=0.154D0
C...
C...     COLUMN CROSS SECTIONAL AREA (M**2)
         S=PI*(D**2)/4.0D0
C...
C...     GAS VOLUMETRIC FLOW RATE (M**3/MIN)
         V=11.04D0
C...
C...     GAS CONSTANT (M**3-ATM/KG MOL-K) (NOTE - ONE KG MOL OF GAS
C...     OCCUPIES 22.4 M**3 AT 1 ATM AND 273.16 K)
         R=22.4D0*1.0D0/273.16D0
C...
C...     TEMPERATURE (K)
         TK=273.16D0+20.0D0
C...
C...     PRESSURE (ATM)
         P=1.0D0
C...
C...     GAS MOLAR FLOW RATE (KG MOL/MIN-M**2)
         G=V*P/(R*TK)*(1.0D0/S)
C...
C...     GAS MOLAR FLOW RATE ON SOLUTE FREE BASIS (KG MOLS AIR
C...     /MIN-M**2)
         GS=G*(1.0D0/(1.0D0+YA1))
C...
C...     LIQUID MOLAR FLOW RATE ON SOLUTE FREE BASIS (KG MOLS H20
C...     /MIN-M**2)
         LS=11.18D0/18.0D0*(1.0D0/S)
C...
C...  EXITING LIQUID COMPOSITION
      XA1=OPLN(YA1,YA2,XA2,GS,LS)
C...
C...  COLUMN LENGTH (M)
      ZL=3.28D0
C...
C...  WRITE SUMMARY OF OPERATING CONDITIONS
      WRITE(NO,1)YA1,XA2,YA2,XA1,LS,GS
1     FORMAT(' Summary of Operating Conditions',//,
     +       '    Entering gas     = ',F7.4,/,
     +       '    Entering liquid  = ',F7.4,/,
     +       '    Exiting gas      = ',F7.4,/,
     +       '    Exiting liquid   = ',F7.4,/,
     +       '    Liquid flow rate = ',F7.4,' kg mols H20/min-m**2',/,
     +       '    Gas flow rate    = ',F7.4,' kg mols air/min-m**2',/)
C...
C...  COEFFICIENTS IN THE LEAST SQUARES SECOND ORDER POLYNOMIAL
C...  APPROXIMATION OF THE EQUILIBRIUM DATA (THE COEFFICIENTS ARE
C...  WRITTEN INTO COMMON/E/ BY THIS CALL TO COEFF)
      CALL COEFF
C...
C...  MASS TRANSFER COEFFICIENT (KG MOL AIR/M**2-MIN)
      KYA=73.8D0
C...
C...  COLUMN HEIGHT USING LM MOL RATIO CONCENTRATION DIFFERENCE
      WRITE(NO,2)
2     FORMAT(///,' Column height from log mean concentration',
     +           ' differenece')
      Z=HTLM(GS,KYA,XA1,YA1,XA2,YA2)
      WRITE(NO,3)Z
3     FORMAT(/,' Column height = ',F6.2,' m',//)
C...
C...  COLUMN HEIGHT USING INTEGRATION  OF MOL RATIO CONCENTRATION
C...  DIFFERENCE
      WRITE(NO,4)
4     FORMAT(///,' Column height from numerical integration of',
     +           ' concentration difference')
      Z=HTIN(GS,KYA,YA1,YA2)
      WRITE(NO,5)Z
5     FORMAT(/,' Column height = ',F6.2,' m',//)
C...
C...  INITIAL CONDITIONS FOR PDE SOLUTION
      DO 10 I=1,NZ
         XA(I)=0.0D0
         YA(I)=0.0D0
10    CONTINUE
      RETURN
      END
      SUBROUTINE DERV
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION LS, KYA
C...
C...  COMMON AREA TO LINK SUBROUTINES
      PARAMETER(NZ=21)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/  XA(NZ),  YA(NZ)
     +      /F/ XAT(NZ), YAT(NZ)
     +      /S/ XAZ(NZ), YAZ(NZ), YAS(NZ)
     +      /C/     XA1,     YA1,     XA2,     YA2,
     +       ZL,     GS,      LS,     KYA,       S
C...
C...  BOUNDARY CONDITION AT Z = 0
      YA(1)=YA1
C...
C...  BOUNDARY CONDITION AT Z = ZL
      XA(NZ)=XA2
C...
C...  XAZ, YAZ
      CALL DSS020(0.0D0,ZL,NZ,XA,XAZ,-LS)
      CALL DSS020(0.0D0,ZL,NZ,YA,YAZ, GS)
C...
C...  YA* GAS COMPOSITION (FROM EQUILIBRIUM LINE)
      DO 1 I=1,NZ
         YAS(I)=EQUIL(XA(I))
1     CONTINUE
C...
C...  PDES
      DO 2 I=1,NZ
         YAT(I)=-GS*YAZ(I)+KYA*(YAS(I)-YA(I))
         XAT(I)= LS*XAZ(I)-KYA*(YAS(I)-YA(I))
2     CONTINUE
       YAT(1)=0.0D0
      XAT(NZ)=0.0D0
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION LS, KYA
C...
C...  COMMON AREA TO LINK SUBROUTINES
      PARAMETER(NZ=21)
      COMMON/T/       T,   NSTOP,   NORUN
     +      /Y/  XA(NZ),  YA(NZ)
     +      /F/ XAT(NZ), YAT(NZ)
     +      /S/ XAZ(NZ), YAZ(NZ), YAS(NZ)
     +      /C/     XA1,     YA1,     XA2,     YA2,
     +       ZL,     GS,      LS,     KYA,       S
C...
C...  ARRAYS FOR PLOTTING
      PARAMETER(NP=21)
      DIMENSION TP(NP), XAP(NP), YAP(NP)
C...
C...  MONITOR THE OUTPUT
      WRITE(*,*)T,YA(NZ)
C...
C...  WRITE A HEADING FOR THE SOLUTION
      IF(T.LT.0.001D0)THEN
         WRITE(NO,1)
1        FORMAT(/,9X,'t',3X,'Xa(L,t)',3X,'Xa(0,t)',
     +                   3X,'Ya(0,t)',3X,'Ya(L,t)',/)
         IP=0
      END IF
C...
C...  WRITE THE SOLUTION
      WRITE(NO,2)T,XA(NZ),XA(1),YA(1),YA(NZ)
2     FORMAT(F10.2,4F10.5)
C...
C...  STORE THE SOLUTION FOR PLOTTING
      IP=IP+1
      TP(IP)=T
      XAP(IP)=XA(1)
      YAP(IP)=YA(NZ)
C...
C...  PLOT THE SOLUTION
      IF(IP.LT.NP)RETURN
C...
C...     BOTTOM LIQUID COMPOSITION VS T
         CALL SPLOTS(1,NP,TP,XAP)
         WRITE(NO,3)
3        FORMAT(/,' Xa(0,t) vs t')
C...
C...     TOP GAS COMPOSITION VS T
         CALL SPLOTS(1,NP,TP,YAP)
         WRITE(NO,4)
4        FORMAT(/,' Ya(L,t) vs t')
      RETURN
      END
      DOUBLE PRECISION FUNCTION OPLN(YA,YA2,XA2,GS,LS)
C...
C...  Function OPLN computes the bulk liquid composition from
C...  the corresponding gas composition using the equation for
C...  the operating line based on the top compositions (at end
C...  2)
C...
C...     Xa = (Gs/Ls)*(Ya - Ya2) + Xa2
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION LS
C...
C...  EQUATION FOR THE OPERATING LINE BASED ON THE COMPOSITIONS AT
C...  THE TOP OF THE COLUMN (END 2)
      XA=(GS/LS)*(YA-YA2)+XA2
C...
C...  EVALUATE FUNCTION OPLN
      OPLN=XA
      RETURN
      END
      DOUBLE PRECISION FUNCTION EQUIL(XA)
C...
C...  Function EQUIL computes the gas phase composition in equilibroum
C...  with the liquid in the column, using the least squares second
C...  order polynomial evaluated by subroutine COEFF.
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE
      COMMON/E/A(2)
C...
C...  SECOND ORDER LEAST SQUARES POLYNOMIAL
      EQUIL=A(1)*XA+A(2)*(XA**2)
      RETURN
      END
      DOUBLE PRECISION FUNCTION HTLM(GS,KYA,XA1,YA1,XA2,YA2)
C...
C...  Function HTLM compuytes the height of a column using the log-
C...  mean concentration difference, according to the equation
C...
C...                                      *
C...     z = (Gs/Kya)*(Ya1 - Ya2)/(Ya - Ya )
C...                                        lm
C...  where
C...                              *               *
C...             *      (Ya1 - Ya1 )  - (Ya2 - Ya2 )
C...     (Ya - Ya )   = ----------------------------
C...               lm                      *
C...                             (Ya1 - Ya1 )
C...                          ln ------------
C...                                       *
C...                             (Ya2 - Ya2 )
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION KYA
C...
C...  COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE
      COMMON/E/A(2)
C...
C...  INPUT/OUTPUT (I/O) UNIT NUMBERS
      COMMON/IO/ NI, NO
C...
C...  EQUILIBRIUM GAS PHASE COMPOSITION AT THE BOTTOM OF THE COLUMN
      YA1S=EQUIL(XA1)
C...
C...  EQUILIBRIUM GAS PHASE COMPOSITION AT THE TOP OF THE COLUMN
      YA2S=EQUIL(XA2)
C...
C...  LOG MEAN CONCENTRATION DIFFERENCE
         DELTA1=YA1-YA1S
         DELTA2=YA2-YA2S
         DYLM=(DELTA1-DELTA2)/DLOG(DELTA1/DELTA2)
C...
C...  WRITE LOG MEAN CONCENTRATION DIFFERENCE
      WRITE(NO,1)DYLM
1     FORMAT(//,' Log-mean concentration difference = ',F9.6,/)
C...
C...  WRITE NUMBER OF TRANSFER UNITS
      WRITE(NO,2)(YA1-YA2)/DYLM
2     FORMAT(/,' Number of transfer units = ',F6.2,/)
C...
C...  WRITE HEIGHT OF TRANSFER UNIT (FT)
      WRITE(NO,3)GS/KYA
3     FORMAT(/,' Height of transfer unit = ',F6.3,' m',/)
C...
C...  COLUMN HEIGHT
      HTLM=(GS/KYA)*(YA1-YA2)/DYLM
      RETURN
      END
      DOUBLE PRECISION FUNCTION HTIN(GS,KYA,YA1,YA2)
C...
C...  Function HTIN computes the height of a column using the design
C...  integral
C...
C...                  Ya2   dYa
C...     z = (Gs/Kya) INT --------
C...                  Ya1 Ya - Ya*
C...
C...
C...  The integration is done by Simpson's rule, implemented in
C...  function SIMP
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION LS, KYA

C...  COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE
      COMMON/E/A(2)
C...
C...  INPUT/OUTPUT (I/O) UNIT NUMBERS
      COMMON/IO/ NI, NO
C...
C...  EXTERNAL FUNCTION FOR INTEGRAND
      EXTERNAL YAINT
C...
C...  EVALUATE DESIGN INTEGRAL.  NOTE THAT A SIGN REVERSAL IS
C...  REQUIRED SINCE DYA IN THE INTEGRAL IS NEGATIVE WHILE
C...  1/(YA - YA*) IS POSITIVE
      NP=100
      TU=-SIMP(YAINT,YA1,YA2,NP)
C...
C...  WRITE NUMBER OF TRANSFER UNITS
      WRITE(NO,2)TU
2     FORMAT(//,' Number of transfer units = ',F6.2,/)
C...
C...  WRITE HEIGHT OF TRANSFER UNIT (M)
      WRITE(NO,3)GS/KYA
3     FORMAT(/,' Height of transfer unit = ',F6.3,' m',/)
C...
C...  COLUMN HEIGHT
      HTIN=(GS/KYA)*TU
      RETURN
      END
      DOUBLE PRECISION FUNCTION YAINT(YA)
C...
C...  Function YAINT computes the integrand of the integral
C...
C...
C...     Ya2   dYa
C...     INT --------
C...     Ya1 Ya - Ya*
C...
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION LS, KYA
      COMMON/C/     XA1,     YA1,     XA2,     YA2,
     +       ZL,     GS,      LS,     KYA,       S
C...
C...  LIQUID COMPOSITION FROM THE OPERATING LINE
      XA=OPLN(YA,YA2,XA2,GS,LS)
C...
C...  EQUILIBRIUM GAS COMPOSITION
      YAS=EQUIL(XA)
C...
C...  INTEGRAND
      YAINT=1.0D0/(YA-YAS)
      RETURN
      END
      DOUBLE PRECISION FUNCTION SIMP(F,A,B,NP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  EXTERNAL FUNCTION FOR INTEGRAND
      EXTERNAL F
C...
C...  CHECK FOR AN EVEN NUMBER OF PANELS
      IF((NP/2*2).NE.NP)THEN
         WRITE(*,2)NP
2        FORMAT(//,' NP = ',I3,' (NOT EVEN)')
         SIMP=0.0D0
         RETURN
      END IF
C...
C...  INTEGRATION INTERVAL
      DX=(B-A)/DFLOAT(NP)
C...
C...  INITIALIZE SUM
      SIMP=(F(A)+4.0D0*F(A+DX)+F(B))
C...
C...  CONTINUE SUM
      IF(NP.GT.2)THEN
         DO 1 I=2,NP-1,2
            SIMP=SIMP
     +          +2.0D0*F(A+DFLOAT(I  )*DX)
     +          +4.0D0*F(A+DFLOAT(I+1)*DX)

1        CONTINUE
      END IF
C...
C...  INTEGRAL
      SIMP=(DX/3.0D0)*SIMP
      RETURN
      END
      SUBROUTINE COEFF
C...
C...  Subroutine COEFF tests the least squares second order polynomial
C...  fit of the NH3-H20-air equilibrium data in Welty, et al, (1984),
C...  Example 2, pp 680-681.
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE
      COMMON/E/A(2)
C...
C...  INPUT/OUTPUT (I/O) UNIT NUMBERS
      COMMON/IO/ NI, NO
C...
C...  ARRAYS FOR FIVE DATA POINTS
      PARAMETER(N=5)
      DIMENSION X(N), Y(N), YC(N)
C...
C...  EQUILIBRIUM DATA (FROM WELTY, ET AL)
      DATA X/0.0164D0, 0.0252D0, 0.0349D0, 0.0455D0, 0.0722D0/
      DATA Y/0.021D0 , 0.032D0 , 0.042D0 , 0.053D0 , 0.08D0  /
C...
C...  EVALUATE COEFFICIENTS IN THE LEAST SQUARES, SECOND ORDER
C...  POLYNOMIAL
      CALL POLY(X,Y,N)
C...
C...  WRITE THE COEFFICIENTS
      WRITE(NO,1)A(1),A(2)
1     FORMAT(10X,' a1 = ',F8.4,'   a2 = ',F8.4,/)
C...
C...  HEADING FOR THE COMPUTED RESULTS
      WRITE(NO,2)
2     FORMAT('    i',6X,'XNH3',6X,'YNH3',6X,'YNH3',7X,'o/o',/,
     +            5X,7X,'tab', 7X,'tab', 6X,'calc',6X,'diff',/)
C...
C...  COMPUTE VALUES OF THE GAS PHASE COMPOSITION FOR COMPARISON WITH
C...  THE TABULATED VALUES
      DO 5 I=1,5
C...
C...     GAS PHASE COMPOSITION FROM THE SECOND ORDER POLYNOMIAL
         YC(I)=A(1)*X(I)+A(2)*X(I)**2
C...
C...     COMPUTE AND PRINT THE GAS PHASE COMPOSITIONS
         DIFF=(YC(I)-Y(I))/Y(I)*100.0D0
         WRITE(NO,4)I,X(I),Y(I),YC(I),DIFF
4        FORMAT(I5,3F10.4,F10.1)
C...
C...  NEXT POINT
5     CONTINUE
C...
C...  END OF CALCULATION
      RETURN
      END
      SUBROUTINE POLY(X,Y,N)
C...
C...  Subroutine POLY computes the coefficients in the least squares
C...  second order polynomial approximating the equilibrium data.
C...
C...  Consider the least squares analysis of data using a polynomial
C...  model
C...
C...     y = a0 + a1*x + a2*(x**2) + ...                        (1)
C...
C...  For purposes of discussion, we consider a second order poly-
C...  nomial.  The extension to higher order polynomials is straight-
C...  forward.
C...
C...  Application of the least squares principle to N data gives
C...
C...                        N
C...     min E(a0,a1,a2) = sum (y(i) - ym(i))**2
C...                       i=1
C...
C...                N
C...             = sum (y(i) - (a0 + a1*x(i) + a2*(x(i)**2)))**2   (2)
C...               i=1
C...
C...  We can find the minimum of eq. (2), i.e., the values of a0, a1
C...  and a2 which minimize E, by the usual methods of differential
C...  calculus.
C...
C...  Before proceeding, however, we consider the special case for
C...  which the model of eq. (1) is required to pass through the
C...  point y(0) = 0, x(0) = 0.  Thus we see a0 = 0 in eq. (1), and
C...  we therefore have to evaluate only a1 and a2 in eq. (2).
C...
C...  The required partial derivatives (denoted with a "d" because
C...  of the limitations of the character set) are set to zero
C...
C...  dE/da1 =
C...
C...       N
C...    2 sum (y(i) - (a1*x(i) + a2*(x(i)**2)))*(-x(i)) = 0        (3)
C...      i=1
C...
C...
C...  dE/da2 =
C...
C...       N
C...    2 sum (y(i) - (a1*x(i) + a2*(x(i)**2)))*(-x(i)**2) = 0     (4)
C...      i=1
C...
C...  Eqs. (3) and (4) can be rearranged to
C...
C...         N                N             N
C...     a1 sum x(i)**2 + a2 sum x(i)**3 = sum x(i)*y(i)           (5)
C...        i=1              i=1           i=1
C...
C...         N                N             N
C...     a1 sum x(i)**3 + a2 sum x(i)**4 = sum (x(i)**2)*y(i)      (6)
C...        i=1              i=1           i=1
C...
C...  Eqs. (5) and (6) are two linear algebraic equations for a1 and
C...  a2.  The following code computes a1 and a2 from a vector of data
C...  of length N.
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  COEFFICIENTS IN POLYNOMIAL FOR EQUILIBRIUM LINE
      COMMON/E/A(2)
C...
C...  INPUT/OUTPUT (I/O) UNIT NUMBERS
      COMMON/IO/ NI, NO
C...
C...  ARRAYS FOR X, Y
      DIMENSION X(N), Y(N)
C...
C...  COMPUTE THE ELEMENTS OF THE COEFFICIENT MATRIX AND RHS VECTOR OF
C...  EQS. (5) AND (6)
      A11=0.0D0
      A12=0.0D0
      A21=0.0D0
      A22=0.0D0
       B1=0.0D0
       B2=0.0D0
      DO 1 I=1,N
         A11=A11+X(I)**2
         A12=A12+X(I)**3
         A21=A12
         A22=A22+X(I)**4
         B1=B1+X(I)*Y(I)
         B2=B2+(X(I)**2)*Y(I)
1     CONTINUE
C...
C...  COMPUTE A1 AND A2 IN EQ. (1) (BY DETERMINANTS)
      DET=A11*A22-A12*A21
C...
C...  TEST IF THE ALGEBRAIC EQUATIONS ARE NUMERICALLY SINGULAR
      IF(DABS(DET).LT.1.0D-20)THEN
C...
C...     DETERMINANT IS SMALL SO TERMINATE THE CALCULATION.  THE
C...     PRECEDING THRESHOLD OF 1.0D-20 WAS SELECTED ARBITRARILY
C...     AND IS DETERMINED PRIMARILY BY THE WORDLENGTH OF THE
C...     COMPUTER
         WRITE(NO,2)DET
2        FORMAT(' Determinant of coefficient matrix = ',D11.3)
         STOP
      END IF
C...
C...  SYSTEM IS NOT NUMERICALLY SINGULAR (ILL-CONDITIONED), SO CONTINUE
C...  THE SOLUTION OF THE ALGEBRAIC SYSTEM
      DET1=B1*A22-B2*A12
      DET2=B2*A11-B1*A21
C...
C...  COMPUTE COEFFICIENTS A1 AND A2 FROM THE DETERMINANTS (BY KRAMER'S
C...  RULE, WHICH IS GENERALLY NOT RECOMMENDED FOR THE SOLUTION OF
C...  SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS, BUT IS ADEQUATE FOR THIS
C...  SMALL 2 X 2 PROBLEM)
      A(1)=DET1/DET
      A(2)=DET2/DET
      RETURN
      END
