*DECK SEP10
      SUBROUTINE INITAL
C...
C...  PACKED COLUMN ABSORBER DYNAMICS
C...
C...  14.5 KG MOLS/HR OF A GAS MIXTURE CONTAINING 6 PERCENT SO2 AND
C...  94 PERCENT DRY AIR IS TO BE SCRUBBED COUNTERCURRENTLY WITH
C...  985 KG MOLS/HR OF FRESH WATER IN A PACKED COLUMN.  INITIALLY,
C...  THE CONCENTRATION OF SO2 IN THE GAS AND LIQUID PHASES IN THE
C...  COLUMN IS ZERO.  THE COLUMN OPERATES ISOTHERMALLY AT 30.5 C AND
C...  1 ATM.  TWO EQUILBRIUM RELATIONS ARE TO BE INVESTIGATED
C...
C...     YE = 296*X**1.33 (NONLINEAR EQUILIBRIUM)
C...
C...     YE = 25.7*X - 0.00063 (LINEAR EQUILIBRIUM)
C...
C...  THUS, TWO RUNS ARE PROGRAMMED (YE IS THE GAS COMPOSITION IN
C...  EQUILIBRIUM WITH THE LIQUID COMPOSITION X).
C...
C...  THE PACKED COLUMN IS 3 M HIGH AND HAS A CROSS-SECTIONAL AREA
C...  OF 0.45 M**2.  THE MOLAR LIQUID AND GAS HOLDUPS PER UNIT VOLUME
C...  OF COLUMN ARE 8.0 AND 0.16 KG MOLS/M**3, RESPECTIVELY.  THE
C...  OVERALL MASS TRANSFER COEFFICIENT AT THESE CONDITIONS IS
C...
C...     KY*AV = 93.0 KG MOLS/(HR-MOL FRACTION-M**3 COLUMN)
C...
C...  THE EXITING CONCENTRATIONS, X(0,T) AND Y(ZL,T), ARE TO BE
C...  COMPUTED AS A FUNCTION OF TIME.
C...
C...  THE MATHEMATICAL MODEL FOR THE COLUMN CONSISTS OF TWO SIMUL-
C...  TANEOUS MATERIAL BALANCE PDES
C...
C...     H*X  =  (L/S)*X  - KY*AV*(YE - Y)                      (1)
C...        T           Z
C...
C...     G*Y  = -(V/S)*Y  + KY*AV*(YE - Y)                      (2)
C...        T           Z
C...
C...  WHERE
C...
C...     T     TIME, HR
C...
C...     Z     DISTANCE FROM THE GAS INLET, M
C...
C...     Y     MOL FRACTION OF SOLUTE (SO2) IN THE GAS
C...
C...     X     MOL FRACTION OF SOLUTE (SO2) IN THE LIQUID
C...
C...     S     COLUMN CROSS-SECTIONAL AREA, M**2
C...
C...     V     GAS MOLAR FLOW RATE, KG MOL/HR
C...
C...     L     LIQUID MOLAR FLOW RATE, KG MOL/HR
C...
C...     G     GAS MOLAR HOLDUP PER UNIT VOLUME OF COLUMN,
C...           KG MOLS/M**3
C...
C...     H     LIQUID MOLAR HOLDUP PER UNIT VOLUME OF COLUMN,
C...           KG MOLS/M**3
C...
C...  THE INITIAL AND BOUNDARY CONDITIONS FOR EQUATIONS (1) AND (2)
C...  ARE
C...
C...     X(Z,0) = X1(Z), X(ZL,T) = X2(T)                        (3)(4)
C...
C...     Y(Z,0) = Y1(Z), Y( 0,T) = Y2(T)                        (5)(6)
C...
C...  WHERE X1(Z) AND Y1(Z) ARE FUNCTIONS OF Z DEFINED BY THE PROBLEM
C...  STATEMENT, AND X2(T) AND Y2(T) ARE FUNCTIONS OF T DEFINED BY THE
C...  PROBLEM STATEMENT, I.E.,
C...
C...     X1(Z) = 0, X2(T) = 0
C...
C...     Y1(Z) = 0, Y2(T) = 0.06
C...
C...  THE FOLLOWING CODING IS IN DOUBLE PRECISION FOR TWO REASONS:
C...
C...     (1)  THE SHARP FRONTS PROPAGATING ALONG THE COLUMN, AND THE
C...          RANGE OF NUMBERS, ESPECIALLY THE SMALL VALUES OF X AND
C...          Y, MAY REQUIRE DOUBLE PRECISION CALCULATIONS (THIS WILL,
C...          OF COURSE, DEPEND ON THE WORD LENGTH OF THE COMPUTER AND
C...          THEREFORE YOU MAY NOT OBSERVE ANY DIFFERENCE BETWEEN THE 
C...          SINGLE AND DOUBLE PRECISION RESULTS).
C...
C...     (2)  THE ADDITIONAL PROGRAMMING FOR DOUBLE PRECISION FORMAT
C...          IS ILLUSTRATED.
C...
C...          THE FOLLOWING CODING ILLUSTRATES GENERALLY WHAT HAS TO BE
C...          DONE FOR PROGRAMMING IN DOUBLE PRECISION FORTRAN, I.E.,
C...
C...          (A)  TYPING OF VARIABLES EITHER BY LISTING THEM INDIVID-
C...               UALLY IN A DOUBLE PRECISION STATEMENT, OR BY THE
C...               MORE CONVENIENT STATEMENT
C...
C...               IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...          (B)  CONVERSION OF CONSTANTS BY ADDING SOMETHING LIKE D0
C...               (NOTE THAT CONVERSION BACK TO SINGLE PRECISION IS
C...               EASILY DONE WITH AN EDITOR - D0 IS REPLACED WITH E0).
C...
C...          (C)  CONVERSION OF INTRINSIC FUNCTIONS WHICH ARE NOT 
C...               GENERIC, I.E., CANNOT ACCEPT SINGLE AND DOUBLE PRE-
C...               CISION ARGUMENTS AND PRODUCE CORRESPONDING SINGLE
C....              AND DOUBLE PRECISION RESULTS.  AN EXAMPLE IS THE
C...               CONVERSION OF AMAX1 TO DMAX1.
C...
C...          (D)  CONVERSION OF FORMATS, E.G., E10.3 TO D10.3.  NOTE
C...               THAT THE F FORMAT WILL ACCOMMODATE SINGLE AND DOUBLE
C...               PRECISION NUMBERS AND THEREFORE DOES NOT REQUIRE
C...               CONVERSION, EITHER IN INPUT OR OUTPUT.
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     X(NZ),     Y(NZ)
     2      /F/  DXDT(NZ),  DYDT(NZ)
     3      /S/  DXDZ(NZ),  DYDZ(NZ)
     4      /P/         A,         B,         C,         D,        Y1,
     5                 XN,        ZL,        IP
C...
C...  SET THE COLUMN LENGTH
      ZL=3.0D0
C...
C...  SET THE ENTERING LIQUID AND GAS COMPOSITIONS
      XN=0.D0
      Y1=0.06D0
C...
C...  SET THE INITIAL LIQUID AND GAS COMPOSITIONS
      DO 1 J=1,NZ
      X(J)=0.D0
      Y(J)=0.D0
1     CONTINUE
C...
C...  PRECOMPUTE SOME CONSTANTS FOR USE IN THE PDES IN SUBROUTINE DERV.
C...  ALSO, THE MODEL TIME SCALE IS CONVERTED FROM HR TO MIN
      A=985.D0/(60.0D0*0.45D0*8.0D0)
      B=93.0D0/(60.0D0*8.0D0)
      C=14.5D0/(60.0D0*0.45D0*0.16D0)
      D=93.0D0/(60.0D0*0.16D0)
C...
C...  CALL SUBROUTINE DERV TO INITIALIZE THE DERIVATIVE CALCULATIONS
      CALL DERV
C...
C...  INITIALIZE THE COUNTER FOR THE PLOTTED SOLUTION
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     X(NZ),     Y(NZ)
     2      /F/  DXDT(NZ),  DYDT(NZ)
     3      /S/  DXDZ(NZ),  DYDZ(NZ)
     4      /P/         A,         B,         C,         D,        Y1,
     5                 XN,        ZL,        IP
C...
C...  DIMENSION THE ARRAY FOR THE EQUILIBRIUM GAS COMPOSITIONS
      DIMENSION YE(NZ)
C...
C...  CONSTRAIN LIQUID OR GAS COMPOSITIONS TO NONNEGATIVE VALUES
      DO 3 J=1,NZ
      IF(X(J).LT.0.D0)X(J)=0.D0
      IF(Y(J).LT.0.D0)Y(J)=0.D0
3     CONTINUE
C...
C...  SET THE BOUNDARY (ENTERING) COMPOSITIONS
      X(NZ)=XN
      Y(1) =Y1
C...
C...  COMPUTE THE EQUILIBRIUM GAS COMPOSITIONS
      DO 1 J=1,NZ
C...
C...  USE A NONLINEAR (NORUN = 1) OR LINEAR (NORUN = 2) EQUILIBRIUM
C...  RELATIONSHIP
      IF(NORUN.EQ.1)THEN
         YE(J)=296.D0*X(J)**1.33D0
      ELSE IF(NORUN.EQ.2)THEN
         YE(J)=25.7D0*X(J)-0.00063D0
      END IF
      IF(YE(J).LT.0.D0)YE(J)=0.D0
1     CONTINUE
C...
C...  SPATIAL DERIVATIVES IN THE PDES
      CALL DSS020(0.D0,ZL,NZ,X,DXDZ,-1.D0)
      CALL DSS020(0.D0,ZL,NZ,Y,DYDZ, 1.D0)
C...
C...  PDES
      DO 2 J=1,NZ
      DXDT(J)= A*DXDZ(J)-B*(YE(J)-Y(J))
      DYDT(J)=-C*DYDZ(J)+D*(YE(J)-Y(J))
2     CONTINUE
C...
C...  ZERO THE TEMPORAL DERIVATIVES AT THE BOUNDARIES
      DXDT(NZ)=0.D0
      DYDT(1) =0.D0
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     X(NZ),     Y(NZ)
     2      /F/  DXDT(NZ),  DYDT(NZ)
     3      /S/  DXDZ(NZ),  DYDZ(NZ)
     4      /P/         A,         B,         C,         D,        Y1,
     5                 XN,        ZL,        IP
C...
C...  DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION
      DIMENSION   TP(101),  X1P(101),  YNP(101)
C...
C...  PRINT A HEADING AT THE BEGINNING OF THE RUN
      IF(IP.EQ.0)WRITE(NO,4)
4     FORMAT(1H ,/,2H T,/,5H X(J),/,5H Y(J))
C...
C...  PRINT THE NUMERICAL SOLUTION IN EVERY TENTH CALL TO PRINT
      IP=IP+1
      IF(((IP-1)/10*10).EQ.(IP-1))THEN
      NZI=4
      WRITE(NO,1)T,(X(J),J=1,NZ,NZI),(Y(J),J=1,NZ,NZI)
1     FORMAT(1H ,//,F10.1,//,6F10.6,/,10X,5F10.6,
     1                    //,6F10.6,/,10X,5F10.6)
      END IF
C...
C...  STORE THE SOLUTION FOR SUBSEQUENT PLOTTING.  ANY VALUES OF X AND Y
C...  WHICH ARE SLIGHTLY NEGATIVE ARE CONVERTED TO ZERO BY THE CALL TO
C...  DMAX1
      X1P(IP)=DMAX1(0.D0,X(1) )
      YNP(IP)=DMAX1(0.D0,Y(NZ))
       TP(IP)=T
C...
C...  PLOT THE SOLUTION AT THE END OF THE RUN
      IF(IP.LT.101)RETURN
      CALL TPLOTS(1,IP,TP,X1P)
      WRITE(NO,2)
2     FORMAT(1H ,//,28H  X(0,T) (MOL FR) VS T (MIN))
      CALL TPLOTS(1,IP,TP,YNP)
      WRITE(NO,3)
3     FORMAT(1H ,//,28H Y(ZL,T) (MOL FR) VS T (MIN))
C...
C...  CHECK THE STEADY STATE MATERIAL BALANCE
C...
C...  IF EQUATIONS (1) AND (2) ARE ADDED, WE OBTAIN
C...
C...     H*X  - (L/S)*X  = G*Y  + (V/S)*Y                       (3)
C...        T          Z      T          Z
C...
C...  EQUATION (3) COULD BE CHECKED AT EACH POINT ALONG THE SOLUTION 
C...  SINCE ALL FOUR PARTIAL DERIVATIVES ARE READILY AVAILABLE, I.E.,
C...  X  , X  , Y  AND Y .  HOWEVER, HERE WE CHECK EQUATION (3) ONLY
C...   T    Z    T      Z
C...  AT THE END OF THE RUN WHEN THE SYSTEM HAS REACHED STEADY STATE
C...  FOR WHICH X  = Y  = 0.  THUS
C...             T    T
C...
C...     (L/S)*X  - (V/S)*Y = 0                                 (4)
C...            Z          Z
C...
C...  EQUATION (4) IS EVALUATED AT Z = 0 AND Z = ZL IN THE FOLLOWING
C...  CODING
C...
C...  Z = 0
         TERMX=(985.D0/0.45D0)*DXDZ(1)
         TERMY=(14.5D0/0.45D0)*DYDZ(1)   
         SUM=TERMX-TERMY
         WRITE(NO,10)TERMX,TERMY,SUM
10       FORMAT(//,' TERMX = ',D10.3,
     +             ' TERMY = ',D10.3,
     +             '   SUM = ',D10.3,' AT Z = 0')
C...
C...  Z = ZL
         TERMX=(985.D0/0.45D0)*DXDZ(NZ)
         TERMY=(14.5D0/0.45D0)*DYDZ(NZ)   
         SUM=TERMX-TERMY
         WRITE(NO,11)TERMX,TERMY,SUM
11       FORMAT(//,' TERMX = ',D10.3,
     +             ' TERMY = ',D10.3,
     +             '   SUM = ',D10.3,' AT Z = ZL')
      RETURN
      END

