

*    *******************************************************************
*    *                                                                 *
*    *                           PROGRAM GASWORKS                      *
*    *                              VERSION 1.2                        *
*    *                                                                 *
*    *******************************************************************

*     This program evolved from program chiller, originally written by
*      M. Reed, U.C. Berkeley, 1976.
*     Chiller was further developed by N. Spycher and M. Reed
*     (1980-1988) then adapted to GASWORKS in late 1988 by M. Reed,
*     M. Hirschmann, R. Symonds, and N. Spycher

*    *******************************************************************
*     This program needs data file GASTHERM to run
*    *******************************************************************
*       This version designed to run with GASTHERM which contains
*       regression coefficients as well as log K'S.
*       individual log K'S are not used in GASWORKS.
*     ******************************************************************
*     GASWORKS is used to calculate properties of overall equilibrium in
*       solid-liquid-gas systems upon:
*        1.  Heating or cooling a gas (or fluid-solid-gas mixture) at
*             constant composition.
*        2.  Heating or cooling, as in (1.), but with solid fractionatio
*             (bulk composition change).
*        3.  Mixing of two gases of different temperature and/or
*             composition.
*        4.  Reaction  of a gas with a solid (e.g. rock),
*             causing a chemical reaction.
*     ******************************************************************
*     ******************************************************************
*     Revised:
*       Dec 88 M.R.  reorganized real and integer statements to
*        eliminate a tropical cockroach ghost that munched on FO2.  Also
*        corrected OXY and made numerous changes in main for constant
*        FO2 calculations
*       Dec 88 M.R.  fixed output formats; tried new oxy form (failed,
*        so far); squished underflow bugs involving output of ALA.
*      22 Feb 89. This version sent to Symonds
*      28 Feb 89. LOOC = 3 to stop execution.  MR
*      25 April 89. Gave the FO2 bees their honey (changed PF to PFLUID)
*        and changed OXY back to the old format
*      Jan-Feb 90: (RS) Converted K'S to log K'S for minerals. Changed
*        arrays to handle 60 components. IPUNCH = 2 to plot log(moles)
*        for gases, otherwise log(fugacities) are plotted.
*      Mar-April 90: (RS) updated GASWORKS to handle the new GASTHERM.
*        Changes size of the mineral arrays and added a new POWERK
*        subroutine to calculate log K as F(T).
*      Oct 90  (RS) reset mintry of kicked out mineral to zero (just
*        after 777 statement in main and 195 statement in MINSAT)
*        modified 407 statement so LOOC=3 print plotting file
*      Nov 90 (RS) fixed the "DO 319" loop so that MTRY's of derived
*        species are set correctly before first iteration
*      Dec 90 (RS) converted GASWORKS to Fortran 77.  Changed
*        convergence criteria.
*      March 91 (RS) converted GASWORKS to compile with Lahey Fortran
*         compiler.
*      May 91 (JD) adjustable dimensions for some some arrays
*      May 91 (MR) adjusted output to 80-column width from chiller
*      August 91 (RS) added GASTOT to matrix and to the WORKRUN file,
*         shortened the 313 to 358 zone
*      Oct 91 (RS) modified convergence criteria for FO2 equation
*      Oct 91 (MR,JD) input formats for SINC and TOTMIX changed to allow
*         exponential notation.  added new switches and separate
*         routines for Symonds' and Reed's (MINPLOT) plotting programs
*         (see IPUNCH comment).  Bufout now eliminates arrays (ALGG,
*         ALGM, etc.).  BG3 parameter installed to adjust MAXSAQ.
*      Jan 92  (RS) added BG4 parameter and some documentation on BG
*         parameters.  Corrected error in mole fraction & partial
*         pressure output
*      Mar 92 (RS) corrected error in the calculation of MTOT(H2) if FO2
*         is input.  Also corrected a precision-related error in
*         calculating AQM(NS) values (in the DO 350 loop).
*      May 02 (MR) change I/O unit numbers for IOUT1 and IOUT2
*      2022 May 25: Change format statements to allow lines longer
*         than 72 characters in Gastherm.dat. J.Palandri

*     ******************************************************************

      Program GASWORKS
      Implicit  double precision  (A-H,O-Z)

*---------------------------------------------------------------------|
*     BG1 =  parameter to scale MAXNM  (minerals)                     |
*     BG2 =  parameter to scale MAXNS  (total species)                |
*     BG3 =  parameter to scale MAXSAQ (component species)            |
*            MAXIN, AND MAXOX                                         |
*     BG4 =  parameter to scale MAXN = BG3 + BG3 + 1                  |
*     BG5 =  parameter to scale the size of the SATMIN array          |
*           see box below for more definitions                        |
*---------------------------------------------------------------------|

      Integer BG1,BG2,BG3,BG4,BG5
      Parameter(BG1=450,BG2=650,BG3=42,BG4=85)
      Parameter(BG5=BG3+1)
      Double precision  LGK,MTRY,MTOT,LGKM,MINTRY,LOGACT
      Integer SAQ, QAS, SAQTOT, SPECOX,OFLAG
      Integer ERLOOP, SPEC, BSTOT, SEQ, BSTP1, SPECM,SOLS,QES,BSTP2
      Character*1 BLNK,SUP
      Character*8 AMOLE,BLANK,DUMMY,ETOIL,MIN,NAME,NOMOX,SATMIN,SUPNAM
      Character*8 TEXT,TYP,TITLE1,TITLE2
      Character*80 BUFOUT(BG3)
      Dimension name(BG2),TITLE1(10),TITLE2(10),IT(BG2),DENS(BG1),
     1 W(BG3),X(BG4),SUPNAM(20),FOLD(BG4),
     2 FACT(BG4),CONV(BG4), SATMIN(BG5),COX(8),ISPOX(8)
**     Dimension ACTM(20)
      Common/MAXSIZ/ MAXSAQ,MAXNS,MAXNM,MAXSTO,MAXSSL,MAXN,MAXMEM,MAXIN
      Common/PRECIS/ RMIN, RLOGMN, RLOGMX
      Common /George/ AZERO(BG2),XI, NS, NST
      Common/Morsel/ GAMIN(BG1),B(BG1),SUMSOL(15),SOLTRY(15),
     1 DGDN(BG3,20),SOLS(15,20),ISOLS(BG1),LIMSOL,JSOLS(BG1),IST(15)
      Common/Morris/  MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW
     1 (BG1), tempK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1)
      Common/Mort/MTRY(BG2),GAMMA(BG2),WATM,FO2,FO2log,TERMLN,KK
     1,OFLAG
      Common/Trix/ TMIN(BG3)
      Common/Oxen/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8)
      Common/Bull/ BSTOT,BSTP1,JGSTOT
      Common/Turnip/ SAQ(BG3), QAS(BG3), SAQTOT
      Common/Slough/ A(BG4,BG4),F(BG4), T(BG4,BG4),IV(BG4), N,ITREF
      Common/Froid/ SINC,BIT,SLIM,AKLOG,APFCL(BG2,5),
     1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4),
     2 XDH(4),ALK(5),LGK(BG2),temp,tempC, HEATC,COMTOT(BG3), HEATM,
     3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3),
     4 CVCCF(7,3), IFRAC
      Common/Minnie/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,Plog
     1,GASTOT,NMTT,LOOC
**    Common/Fuga/GASlog(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3,
**   1 NGAS,IMIX,IDEAL,ITGAS
      Common/Io/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
      Common/PP/ TS(200),PS(200),NPP,INCP
      Common/Bess/ WTPC(BG3), COEFOX(BG3,8), XMW(BG3),
     1 NOMOX(BG3),ITOX(BG3),NOOX, SPECOX(BG3,8),MINSOL
      Common/Val/BVAL(20), CVAL(20)
      data ETOIL/'********'/, BLNK/' '/, BLANK/'        '/
      data SUPNAM/20*'        '/, SATMIN/BG5*'        '/,
     1 TYP/'grams   '/,AMOLE/'moles   '/,BUFOUT/BG3*' '/
      data fact/BG4*1./, fold/BG4*1.D+65/

*=====================================================================|
*    ASSIGNMENTS OF DEVICE NUMBERS FOR I/O OPERATIONS                 |
*     INP1 : WORKRUN (data specific to each run, including new        |
*            data after each loop                                     |
*     INP2 : GASTHERM (thermochemical data for gases, solids, and     |
*            liquids)                                                 |
*     INP3 : INPUTS PRESSURE CURVE DATA. *AT PRESENT THIS FILE DOES   |
*            NOT EXIST* THE PROGRAM DOES NOT LOOK FOR THIS FILE AND   |
*            IT SHOULD NOT BE ASSIGNED IN THE WORKEX EXEC(IBM) OR     |
*            .COM (VAX) FILE.                                         |
*     INP4 : oxide or solid data for MINSOL option                    |
*     IPUN : plot file                                                |
*     IOUT1: main output                                              |
*     IOUT2: terminal remarks                                         |
*     IOUT3: pickup file written after a mineral is added or          |
*            subtracted from the matrix                               |
*-------------------------------------------------------------------- |

      INP1 = 4
      INP2 = 1
**    INP3 = 7
      INP4 = 3
      IOUT1 = 9
      IOUT2 = 6
      IOUT3 = 8
      IPUN = 2

*---------------------------------------------------------------------|
*     OPEN FILES INTERNALLY                                           |
*---------------------------------------------------------------------|
*
      OPEN(UNIT=1,FILE='gastherm.dat',STATUS='OLD')
      OPEN(UNIT=2,FILE='workplot.dat',STATUS='UNKNOWN')
      OPEN(UNIT=3,FILE='workox.dat',STATUS='OLD')
      OPEN(UNIT=4,FILE='workrun.dat',STATUS='OLD')
      OPEN(UNIT=9,FILE='workout.dat',STATUS='UNKNOWN')
**    OPEN(UNIT=7,FILE='curve.dat',STATUS='OLD')
      OPEN(UNIT=8,FILE='workrun.sav',STATUS='UNKNOWN')

*---------------------------------------------------------------------|
*     ARRAY MAXIMUM SIZES:                                            |
*        MAXSAQ = component species                                   |
*        MAXNS  = derived gas species                                 |
*        MAXMN  = solids + liquids                                    |
*        MAXTO  = stoichimetric components                            |
*        MAXN   = matrix size                                         |
*        MAXSSL = solid solutions                                     |
*        MAXMEM = solid solution end memembers                        |
*        MAXIN  = maximum number of phases or phase end members       |
*                 (solids + liquids) in matrix                        |
*        MAXOX  = oxides                                              |
*---------------------------------------------------------------------|

      MAXSAQ = BG3
      MAXNS = BG2
      MAXNM = BG1
      MAXSTO = 8
      MAXN = BG4
      MAXSSL = 15
      MAXMEM = 20
      MAXIN = BG3
      MAXOX = BG3

*---------------------------------------------------------------------|
*     COMPUTER PRECISION LIMITS TO AVOID OVERFLOWS, UNDERFLOWS ETC.   |
*     For 386 PCS, smallest double precision number is about          |
*     10**-300 AMD ABOUT 10**300.  For IBM mainframes, these are      |
*     lower.                                                          |
*        RlogMN = minimum log10                                       |
*        RlogMX = maximum log10                                       |
*        RMIN   = minimum real                                        |
*---------------------------------------------------------------------|

      RMIN = 1.D-300
      RlogMN = -300.
      RlogMX = 300.

      ERLOOP= 0
      MINCT= 0
      NNMT = 0
      NST = 0
      NMT = 0
      NDUMP= 0
      JSTOP = 0
      GRAMSO = 0.
      IDEAL = 0
      OFLAG = 0
      FO2 = 0.
      BIT = 0.
      BASE = 10.
      TERMLN = 1./Dlog(BASE)
      DO 10 I = 1,MAXOX
      NOMOX(I) = BLANK
      WTPC(I) = 0.
   10 XMW(I) = 0.
      DO 15 I = 1,MAXSAQ
      TMIN(I) = 0.
      W(I) = 0.
      COMTOT(I) = 0.
      MTOT(I) = 0.
      SAQ(I) = 0
      AQM(I) = 0.
   15 QAS(I) = 0
      DO 16 NS = 1,MAXNS
      DO 6 N = 1,5
    6 APFCL(NS,N) = 0.
      DO 8 I = 1,MAXSTO
      SPEC(NS,I) = 0
    8 COEF(NS,I) = 0.
      NAME(NS) = BLANK
      AZERO(NS) = 0.
      IT(NS) = 0
      MTRY(NS) = 0.
      LGK(NS) = 0.
      ACT(NS) = 0.
   16 GAMMA(NS) = 0.
      DO 17 I = 1,MAXMEM
      CVAL(I) = 0.
   17 BVAL(I) = 0.
      DO 27 NM = 1, MAXNM
      DO 23 N = 1,5
   23 BPFCL(NM,N) = 0.
      DENS(NM) = 0.
      SCALE(NM) = 0.
      GAMIN(NM) = 1.0
      QES(NM) = 0
      LGKM(NM) = 0.
      SEQ(NM) = 0
      B(NM) = 0.
      ISOLS(NM) = 0
      ITM(NM) = 0
      XMOLW(NM) = 0.
      MIN(NM) = BLANK
      JSOLS(NM) = 0
   27 MINTRY(NM) = 0.

      READ(INP1,870) TITLE1
      READ(INP1,870) TITLE2

*---------------------------------------------------------------------|
*     TITLE IN TWO, 80-CHARACTER LINES                                |
*---------------------------------------------------------------------|

      READ(INP1,816)
      READ(INP1,834) ERPC, FO2log,PFLUID, temp, SINC,SLIM,GASTOT
      READ(INP1,816)
      READ(INP1,836) tempC,TOTMIX,SOLMIN,RM,GASGRM, SUPRNT

      IF (FO2log.NE.0.) FO2 = 10.**FO2log

*---------------------------------------------------------------------|
*     ERPC IS THE FRACTIONAL CONVERGENCE LIMIT ON IONIC STRENGTH.
*       IF IT IS NOT SPECIFICALLY INPUT A VALUE OF 1D-12 IS USED
*     IF FO2log EQUALS NOT ZERO, FO2 IS FIXED AT SPECIFIED VALUE. USUALL
*       FO2log IS SET TO ZERO, IN WHICH CASE THE SPECIFIED MTOT(H2)
*       WILL BE USED IN H2 MASS BALANCE EQUATION AND FO2 IS CALCULATED.
*       THE INITIAL MTOT(H2) IS USUALLY DRAWN FROM AN EARLIER SOLVGAS
*       RUN ON THE DESIRED WATER COMPOSITION.
*     PFLUID , IN ATM, IS THE INITIAL PRESSURE OF THE SYSTEM.
*       IF USING WITH NON IDEAL GASES, CHECK FOR PRESSURE LIMITS
*       FOR COMPUTATION OF FUGACITY COEFFICIENTS (SUBROUTINE PHI).
*     SLIM IS THE LAST STEP AFTER WHICH
*       WE WANT TO STOP THE CALCULATION. DEFAULT IS 0 C
*       IF COOLING, 1227 C IF HEATING.
*     Temp IS THE INITIAL tempERATURE FOR THE UNMIXED/UNCOOLED
*       SOLUTION.  INPUT IN DEGREES C.
*     SINC CAN BE (+) OR (-), AND IS THE STEP INCREMENT (T, P, GRAMS, OR
*       MIX FRACTION).  TYPICALLY SINC = -16 DEG.C IS A WORKABLE
*       MAXIMUM STEP.
*       IF USING GASMIX, SINC IS THE FRACTIONAL
*       TITRATION FACTOR, E.G. .01.  THE CHANGE OF MTOTS IS
*       SINC * COMTOT, AS CALCULATED IN MIXER.
*       SINC CANNOT BE ZERO! IF ZERO, IT IS DEFAULTED TO -5 FOR
*       T CHANGE, OR 0.1 IF MIXING
*     tempC IS THE tempERATURE OF THE PHASE TO BE TITRATED IN
*       (USUALLY THE COOLER ONE, BUT NOT NECESSARILY) WHEN USING
*       GASMIX. IF NOT USING GASMIX, SET tempC = 0 !!, THIS IS THE
*       TEST FOR ACTIVATING THE GASMIX OPTION.
*     TOTMIX IS THE TOTAL FRACTION OF MIXING GAS OR OTHER SUBSTANCE
*       (E.G. GAS, ROCK) ADDED WHEN USING MIXING (tempC OR MINSOL
*       NON ZERO).  FOR AN INITIAL RUN,TOTMIX=0; THEREAFTER, TOTMIX
*       IS SIMPLY THE SUM OF ALL PREVIOUS
*       MIXING INCREMENTS (SINCS).  NOTICE THAT TOTMIX (AND SINC, WHEN
*       USING THE MIXING OPTION, BECAUSE TOTMIX IS A SUMMATION OF
*       SINCS) IS SIMPLY A DIMENSIONLESS NUMBER THAT DETERMINES THE
*       TITRATION STEP SIZE FOR GAS-GAS MIXING OR ROCK-GAS MIXING.
*       THE ACTUAL QUANTITY OF TITRATED MATERIAL (IN MOLES OR GRAMS) IS
*       FIXED BY THE VALUES OF COMTOT(I) (SEE BELOW) USED.
*     SOLMIN  MINIMUM MINTRY VALUE ALLOWED FOR GASES OR MINERALS
*      DURING ITERATION.  IF COMPUTED MINTRY GETS SMALLER DURING
*      ITERATIONS, THE MINERAL OR GAS WILL BE KICKED OUT OF THE
*      MATRIX DURING THE ITERATION PROCESS (SAVES TIME, EASES CONVER-
*      GENCE). MUST BE 0 OR NEGATIVE.  DURING ITERATION, A NEGATIVE
*      MINTRY GENERALLY MEANS THE PHASE DOES NOT BELONG ANYMORE.  THIS I
*      NOT ALWAYS THE CASE HOWEVER, AND A TOO LARGE SOLMIN MIGHT RESULT
*      IN SWAPPING A PHASE IN, AND OUT OF THE MATRIX INDEFINITELY. IF
*      SET TO A LARGE NEGATIVE VALUE, EQUILIBRATION WITH NEGATIVE MINTRY
*      IS ALLOWED, REQUIRING ONE MORE EQUILIBRATION STEP TO REMOVE
*      THE PHASE WITH NEGATIVE MINTRY (OLD METHOD). SOLID SOLUTIONS ARE
*      KICKED OUT IF MINTRY IS LESS THAN ZERO.
*      SET SOLMIN > 0 TO SKIP THIS OPTION (EQUILIBRATION W/ NEGATIVE
*      AMOUNTS ALLOWED); SOLMIN > 0 WILL KEEP SOLID SOLUTIONS IN
*      MATRIX, HOWEVER. IN THE EVENT ONE SOL-SOL SHOULD UNDERSATURATE,
*      CONVERGENCE WILL FAIL W/ SOLMIN > 0.
*     RM IS NUMBER BETWEEN 0 AND 1 TO MULTIPLY F'S TO HELP CONVERGENCE.
*       SET TO 0 UNLESS THERE ARE CONVERGENCE PROBLEMS.
*     GASGRM  GRAMS OF STARTING GAS PHASE, TO COMPUTE WATER/GAS
*       RATIO IF MINSOL OPTION IS SELECTED
*     SUPRNT  DERIVED SPECIES MINIMUM MOLE NUMBER FOR OUTPUT.  ONLY
*       SPECIES WITH MOLE NUMBER GREATER OR EQUAL TO SUPRNT WILL PRINT
*       DEFAULT IS 1.D-10
*---------------------------------------------------------------------|

      IF(SUPRNT.EQ.0.) SUPRNT = 1.D-10

*---------------------------------------------------------------------|
*     PMAX IS USED IN CHILLER TO SET UPPER PRESSURE LIMITS FOR TREATMENT
*     OF NON-IDEAL GASES.  NOT IN USE HERE SO SET TO 10 KB.
*---------------------------------------------------------------------|

      PMAX = 10000.

      IF(PFLUID.EQ.0.) PFLUID = 1.
      IF(ERPC.EQ.0.) ERPC = 1.0D-12
      READ(INP1,816)
*     COLUMN           5     10      15     20     25
      READ(INP1,800) BSTOT, IFRAC, IPUNCH, NLOOP, ISTEP

     1, LIMSOL, LOOC, ITREF, IDEAL, INCREM, INCP, MINSOL, OFLAG
*     COL 30     35     40     45     50     55     60      65

      IF(IPUNCH.EQ.1.OR.IPUNCH.EQ.2) THEN
      WRITE(IPUN,870) TITLE1
      WRITE(IPUN,870) TITLE2
      END IF

*---------------------------------------------------------------------|
*     BSTOT IS THE NUMBER OF COMPONENT SPECIES IN GASTHERM (E.G. 35).
*     IFRAC = 0 RESULTS IN NO FRACTIONATION DURING COOLING
*     IFRAC = 1 RESULTS IN FRACTIONATION OF MINERALS.
*     IPUNCH NOT ZERO RESULTS IN CREATING A PLOT FILE, BELOW.
*       There are two basic plotting options controlled by the value of
*       IPUNCH.  Each option has two internal variations, thus four
*       non-zero values of IPUNCH are used to control plotting,
*       as follows:
*        IPUNCH       Routine     Effect
*        ------       -------     ------
*         1           Symonds   PLOT FILE WILL HAVE log(FUGACITIES)
*                               FOR GAS SPECIES
*         2           Symonds   PLOT FILE WILL CONTAIN log (MOLES)
*         3           Minplot   Use 3 only for very first equilibration
*                               of very first run.  It puts the header
*                               on the plot file, then IPUNCH is reset
*                               to 4.
*         4           Minplot   Use for all "pickup" runs, after the
*                               first plot file segment has been
*                               written.
*     NLOOP = 0 NO CALCULATIONS ARE PERFORMED. PROGRAM STOPS AFTER
*       ECHOING ALL INPUT DATA, ALL SELECTED INPUT PARAMETERS, AND AN
*       INDEXED LIST OF ALL GAS SPECIES, AND MINERALS. USEFUL TO
*       GET MINERAL INDEXES (SEQ), AND CHECK ALL INPUT DATA.
*       NLOOP > 0 LIMIT ON THE NUMBER OF NEWTON ITERATIONS ALLOWED
*       BEFORE ARBITRARILY STOPPING.  70 IS A GOOD LIMIT.  USUALLY 8-20
*       LOOPS ARE USED.  SOMETIMES MORE THAN 70 ARE NEEDED, EVEN FOR
*       A SOLVABLE CASE, E.G., 120 MIGHT SOMETIMES BE NEEDED AND MAY BE
*       TRIED ON OCCASSION TO WORK THROUGH A COMPUTATIONAL KNOT.
*     ISTEP NOT ZERO RESULTS IN PRINTING OF EXTRA OUTPUT FOR CLOSER
*       FOLLOWING OF CALCULATION PROGRESS (USE ONLY FOR DEBUGGING).
*     LIMSOL NOT ZERO ALLOWS INCLUSION OF SOLID SOLUTIONS (OR MIXED
*       GASES). IF LIMSOL = 0, THE AQUEOUS PHASE ONLY IS OF VARIABLE
*       COMPOSITION (THEREFORE ENDMEMBERS ARE TREATED AS REGULAR PURE
*       MINERALS).
*       LIMSOL = 2 RESULTS IN KICKING OUT OF THE MATRIX ENDMEMBERS
*       WITH CALCULATED AMOUNTS LESS THAN 0.  THIS IS THE ONLY WAY TO
*       THROW OUT A SOLSOL THAT BECOMES UNDERSATURATED, BUT MIGHT
*       RESULT IN KICKING OUT SOLSOLS THAT ARE STILL SATURATED.
*       LIMSOL = 1 WILL KEEP SOLSOL IN, AND CONVERGENCE WILL FAIL IF
*       A SOLSOL BECOME UNDERSAT.
*     LOOC  DETERMINES MINERAL SELECTION METHOD IF MORE THAN ONE
*       MINERAL SUPERSATURATE.  LOOC = 0 RESULTS IN SELECTION OF MOST
*       SUPERSATURATED MINERAL.  LOOC NOT 0 AND NOT 3
*       RESULTS IN BACK STEPPING
*       (REVERSE T INCREMENT OR "BACK-MIX") TO THE POINT WHERE ONLY ONE
*       MINERAL SUPERSATURATES.
*       LOOC = 3 CAUSES EXECUTION TO STOP AFTER THE FIRST EQUILIBRATION.
*     ITREF = 0  NO ITERATIVE REFINEMENT IS USED WHEN SOLVING
*                THE EQUATIONS. THIS SHOULD ALMOST ALWAYS WORK.
*           = 1  NO ITERATIVE REFINEMENT IS USED, BUT THE
*                RESIDUAL ERRORS ARE CALCULATED. ISTEP SHOULD
*                BE SET TO 1 TO PRINT THESE VALUES.
*           > 1  THE RESIDUAL ERRORS ARE CALCULATED AND MINIMISED
*                BY ONE OR MORE ITERATIVE REFINEMENTS. THIS
*                OPTION SHOULD BE USED ONLY IN CASES WHEN CONVER-
*                GENCE DOES NOT OCCUR DUE TO MACHINE PRECISION.
*      IDEAL  = 0  NON IDEAL GASES AND NON IDEAL MIXING  ASSUMED
*             = 1  NON IDEAL GASES AND IDEAL MIXING ASSUMED
*             = 2  IDEAL GASES AND IDEAL MIXING ASSUMED
*     INCREM NOT ZERO PRODUCES IMMEDIATE TITRATION OF SINC-SIZED CHUNK
*       OF REACTANT, BEFORE ANY EQUILIBRATIONS. RELEVANT ONLY
*       IF tempC IS NON ZERO.  USE TO JUMP AHEAD IN TITRATION.
*       CANNOT BE USED WITH FRACTIONATION, AND WILL RESET IFRAC TO 0
*     INCP = 1: P DECREASE AT CONSTANT T. USE SINC AND SLIM WITH
*       P VALUES INSTEAD OF T.  THIS OPTION IS OVERWRITTEN BY OTHER P
*       VARIATIONS OPTIONS (SEE COOLER).
*     MINSOL:  USED TO ENABLE SOLID TITRATION FOR WATER/ROCK REACTIONS
*        tempC MUST BE SET TO ZERO FOR THIS OTION
*        = 0  DISABLED
*        = 1  TITRATE SINC AMOUNT OF SOLIDS IN GRAMS
*        = 2  SAME AS 1, BUT STOPS TITRATING A MINERAL IF EQUILIBRATED
*        = 3  TITRATE SINC AMOUNT OF SOLIDS IN MOLES
*        = 4  SAME AS 3, BUT STOPS TITRATING A MINERAL IF EQUILIBRATED
*
*     OFLAG:  DETERMINES WHETHER FO2 IS ALLOWED TO CHANGE AFTER FIRST
*        EQUILIBRATION.
*        >0  FO2 DOES NOT CHANGE
*        =0  FO2 IS CALCULATED AFTER EACH tempERATURE CHANGE OR
*            TITRATION
*       THERE ARE 4 POSSIBILITIES:
*        1.  OFLAG> 0 AND FO2 IS GIVEN IN WORKRUN.
*        2.  OFLAG> 0 AND FO2 IS NOT GIVEN IN WORKRUN.
*        3.  OFLAG= 0 AND FO2 IS GIVEN IN WORKRUN.
*        4.  OFLAG= 0 AND FO2 IS NOT GIVEN IN WORKRUN.
*       THE RESULTS OF WHICH WILL BE:
*        1.  FO2 IS HELD AT GIVEN VALUE THROUGHOUT RUN.
*        2.  FO2 IS CALCULATED DURING FIRST EQUILIBRATON AND HELD
*               AT THAT VALUE FOR THE REST OF THE RUN.
*        3.  FO2 IS AT GIVEN VALUE DURING FIRST EQUILIBRATION, BUT
*               ALLOWED TO VARY DURING REST OF RUN.
*        4.  FO2 IS ALWAYS RECALCULATED.
*----------------------------------------------------------------------|

      IF(SINC.NE.0.) GO TO 22
      SINC = -5.
      IF(tempC.NE.0..OR.MINSOL.NE.0) SINC = 0.1
   22 IF(SLIM.NE.0.) GO TO 25
      IF(MINSOL.EQ.0.AND.tempC.EQ.0..AND.SINC.LT.0.)  SLIM = 25.
      IF(MINSOL.EQ.0.AND.tempC.EQ.0..AND.SINC.GT.0.)  SLIM = 1227.
   25 CONTINUE
      IF(MINSOL.GT.2) TYP = AMOLE
      IF(INCREM.GT.0) IFRAC = 0
      IMINT=BSTOT
      LIMSS=LIMSOL
*      tempC NON ZERO  GASMIX OPTION ON, CANNOT HAVE MINSOL ON TOO
      IF(tempC.NE.0.) MINSOL = 0
   24 CONTINUE

*     NOW WE READ DATA RELATED TO THE CHEMICAL COMPOSITION OF THE SYSTEM

      READ(INP1,816)
      I = 1
  100 READ(INP1,803)  SAQ(I), MTOT(I), MTRY(I), GAMMA(I), COMTOT(I)

*     SAQ SPECIFIES THE SEQUENCE POSITIONS OF THE COMPONENT SPECIES
*       TO BE DRAWN FROM GASTHERM FOR THIS CALCULATION.  E.G., IF YOU
*       WANT H2,H2O,CO2,H2S,HCL,HF,HBR,NACL,KCL,FECL3,CD,ALF3,CUCL2,ETC.
*       PUT IN 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, ETC.
*       THESE NUMBERS ARE PICKED FROM THE ORDER OF COMPONENT
*       SPECIES LISTED IN SOLTHERM. H2 AND H2O (1, AND 2) MUST BE
*       ENTERED FIRST, IN THAT ORDER FOR PRESENT VERSION OF GASTHERM.
*       OTHER SPECIES MAY BE ENTERED IN ANY ORDER.
*     MTOT IS THE TOTAL NUMBER OF MOLES OF EACH  BASIS SPECIES IN THE
*       STARTING GAS OR IN STARTING (GAS PLUS SOLIDS).
*     MTRY (I) CONTAINS TRIAL MOLE NUMBERS FOR COMPONENT SPECIES. USUALL
*       THESE ARE KNOWN FROM AN EARLIER SOLVGAS RUN.
*     GAMMA ARE FUGACITY COEFFICIENTS. IF 0 THEY ARE RESET TO 1.
*       IN SOME CASES WHERE NUMERICAL STABILITY IS A PROBLEM, GAMMA(I)
*       NOT EQUAL TO 1 MAY BE NEEDED FOR INITIAL VALUES.
*     COMTOT CONTAINS TOTAL MOLES OF COMPONENT SPECIES IN MIXER
*       PHASE (THE ONE TO BE TITRATED IN), FOR USE IN GASMIX.
**      COMTOT MAY ALTERNATIVELY BE USED TO SPECIFY COMPOSITION OF ANY
**      REACTANT TO BE TITRATED IN.  E.G. TO REACT SOLUTION WITH ONE
**      KILOGRAM OF DOLOMITE, SPECIFY
**        COMTOT (1) = -10.8454 (NEGATIVE 10.8454 MOLES OF H+ IN ONE
**                               KILOGRAM OF DOLOMITE)
**          "    (2) = 0 (NO H20 IN DOLOMITE)
**          "    (3) = 0 (NO CL- IN DOLOMITE)
**          "    (4) = 0 (ETC. FOR ZERO VALULES OF COMTOT(I))
**          "    (5) = 10.8454 (POSITIVE 10.8454 MOLES OF BICARBONATE IN
**                              ONE KILOGRAM OF DOLOMITE)
**          "    (6) = 0
**          "    (7) = 0
**          "    (8) = 0
**          "    (9) = 5.4227 (POSITIVE 5.4227 MOLES OF CALCIUM IN ONE
**                             KILOGRAM OF DOLOMITE)
**          "   (10) = 5.4227 (POSITIVE 5.4227 MOLES OF MAGNESIUM IN ONE
**                             KILOGRAM OF DOLOMITE)
**          " (11 - 29) = 0
**      I.E. 1 KG. OF DOLOMITE CONTAINS -10.8454 MOLES OF H+, +10.454
**      MOLES OF HCO3-, +5.4227 MOLES OF CA++ AND MG++.  THE LIMSOL
**      OPTION, HOWEVER, IS A LOT EASIER TO USE FOR THIS TYPE OF
**      REACTION.

*      READS UNTIL BLANK LINE
      IF(SAQ(I).EQ.0) GO TO 144
      IF(SAQ(I).GT.MAXSAQ) GO TO 143
      J = SAQ(I)
      QAS(J) = I
*     SET DEFAULT VALUES
      IF(GAMMA(I).LE.0.) GAMMA(I) = 1.
      IF(MTRY(I).LE.0.) MTRY(I) = 1.D-10
      I = I + 1
      IF(I.LE.MAXSAQ) GO TO 100
      READ(INP1,805) IDUM
*      (IN CASE WE HAVE EXACTLY MAXSAQ COMPONENTS)
       IF(IDUM.EQ.0) GO TO 144
         WRITE(IOUT1,630) MAXSAQ
         WRITE(IOUT2,630) MAXSAQ
  630 FORMAT(//5X,'WARNING! THE NUMBER OF COMPONENTS GIVEN IN WORKRUN
     1 EXCEEDS THE'/5X,'CURRENT ARRAY DIMENSION (',I2,').  CHECK',
     2 ' WORKRUN DATA FILE.'/)
         STOP
  143    WRITE(IOUT1,631) I, SAQ(I), MAXSAQ
         WRITE(IOUT2,631) I, SAQ(I), MAXSAQ
  631 FORMAT(//5X,'WARNING! SAQ(',I2,') = ',I4,';  MAXIMUM VALUE'
     1 ,' ALLOWED IS ', I3/)
         STOP
**************************************
*     THE DO 123 LOOP THAT FOLLOWS IS TO ESTIMATE THE INITIAL NUMBER
*     OF MOLES OF GAS.  GASTOT IS THE TOTAL NUMBER OF MOLES OF GAS.
  144 SAQTOT = I - 1
      IF(GASTOT.NE.0.) GO TO 124
      GASTOT=0.0
      DO 123 I=1, SAQTOT
  123 GASTOT = GASTOT +MTOT(I)
  124 GSTOLD = GASTOT
**************************************
      tempK = temp + 273.15

*     READS EQUILIBRIUM ASSEMBLAGE (MINERALS AND OR GASES) AND AMOUNTS
*     ACT IS USED HERE AS A DUMMY INPUT VARIALBLE FOR MINTRY, AFTER
*      STATEMENT 165 MINTRY IS SET = ACT
*     MINTRY CONTAINS THE NUMBER OF MOLES OF MINERALS CURRENTLY IN
*       EQUILIBRIUM AT P, T, AND X. IF THIS IS THE INITIAL RUN, NORMALLY
*       NO MINERALS ARE KNOWN TO BE IN EQUILIBRIUM.  IN ANY CASE, EVEN
*       IF CERTAIN MINERALS ARE WRONGLY INCLUDED, GASWORKS WILL THROW
*       THEM OUT.

      READ(INP1,817)
      I = 1
  153 READ(INP1,823) SATMIN(I), ACT(I)

      IF(ACT(I).EQ.0.) ACT(I) = 1.D-9
      IF(SATMIN(I).EQ.BLANK) GO TO 154
      I = I + 1
      IF(I.LE.MAXIN) GO TO 153
      READ(INP1,823) DUMMY
        IF(DUMMY.EQ.BLANK) GO TO 154
        WRITE(IOUT1,652) MAXIN
        WRITE(IOUT2,652) MAXIN
  652   FORMAT(/5X,'TOO MANY SATURATED MINERALS SPECIFIED IN WORKRUN.'
     $  ,' LIMIT =' ,I3//)
        STOP
  154 NMT= I-1

*     READS NAMES AND AMOUNTS OF SOLID REACTANT FOR MINSOL OPTION.
*     MAXOX SOLIDS MAXIMUM. IF NONE, INPUT 3 BLANK LINES.
*     WTPC IN WEIGHT % IF MINSOL IS 1 OR 2, IN MOLES % IF MINSOL IS 3,4
      I = 1
      READ(INP1,817)
  133 READ(INP1,823) NOMOX(I), WTPC(I)

      IF(NOMOX(I).EQ.BLANK) GO TO 132
      I = I + 1
      IF(I.LE.MAXOX) GO TO 133
      READ(INP1,823) DUMMY
        IF(DUMMY.EQ.BLANK) GO TO 132
        WRITE(IOUT1,654) MAXOX
        WRITE(IOUT2,654) MAXOX
  654   FORMAT(/5X,'TOO MANY SOLID REACTANTS SPECIFIED IN WORKRUN.'
     $  ,' LIMIT =' ,I3//)
        STOP
  132 NOOX = I-1

*     NOW WE READ NAMES OF SPECIES, GASES OR MINERALS TO BE SUPPRESSED
*     READS UNTIL BLANK OR EOF
      READ(INP1,817)
      DO 70 NSUP = 1,20
      READ(INP1,822,END=72) SUPNAM(NSUP)

   70 IF(SUPNAM(NSUP).EQ.BLANK) GO TO 72
   72 IF(NSUP.NE.20) NSUP = NSUP-1

*     IF THE MINSOL OPTION IS SELECTED, THE OXIDE OR SOLID DATA
*     FILE IS READ BELOW, AND THE NAMES OF SOLIDS GIVEN ABOVE ARE
*     MATCHED. THE COMPOSITIONS OF NOMOX (REACTANTS) ARE ARBITRARY. USE
*     MINERALS, OXIDES, ETC.  WILL READ UNTIL ALL NAMES MATCH OR EOF.
*     *** WARNING: A REACTANT CANNOT BE ENTERED TWICE IN THE OXIDE DATA
*     FILE.  IF IT IS, THIS MIGHT RESULT IN SKIPPING READING OF FURTHER
*     OXIDE DATA (MESSES UP ICOUNT !).
      IF(MINSOL.EQ.0) GO TO 161
      ICOUNT = 0
  152 READ(INP4,854,END=156) DUMMY,SUP,WTX,ITOT,
     1  (COX(J),ISPOX(J), J = 1, ITOT)

      IF(ITOT.LE.MAXSTO) GO TO 150
          WRITE(IOUT1,659) DUMMY, ITOT, MAXSTO
  659     FORMAT(/5X,'SOLID ',A8,' IN MINOX HAS TOO MANY COMPONENTS'
     $    /5X,'ITOT FOR THAT SOLID IS ',I3,' ; MAXIMUM IS ',I3//)
          STOP
  150 CONTINUE
      DO 155 I = 1,NOOX
  155 IF(NOMOX(I).EQ.DUMMY) GO TO 158
      GO TO 152
*     IF EOF REACHED, WE CHECK WHICH SOLIDS WERE NOT FOUND IN MINOX
  156 CONTINUE
      DO 157 I = 1,NOOX
      IF(XMW(I).NE.0.) GO TO 157
      WRITE(IOUT1,658) NOMOX(I)
      WRITE(IOUT2,658) NOMOX(I)
  658 FORMAT(/5X,'SOLID: ',A8,' NOT FOUND IN MINOX DATA FILE'/)
  157 CONTINUE
      STOP
*
  158 ICOUNT = ICOUNT + 1
      XMW(I) = WTX
      ITOX(I) = ITOT
      DO 159 J = 1, ITOT
      COEFOX(I,J) = COX(J)
      JS = ISPOX(J)
      IF(QAS(JS).GT.0) GO TO 159
       WRITE(IOUT1,661) NOMOX(I), JS
       WRITE(IOUT2,661) NOMOX(I), JS
  661  FORMAT(/5X,'REACTANT ',A8, ' CONTAINS SPECIES NUMBER ',I2,
     1  ' WHICH IS NOT',/5X,'INCLUDED IN THE LIST OF COMPONENT SPECIES',
     2  ' FOR THIS RUN.'//)
          STOP
  159 SPECOX(I,J) = QAS(JS)
      IF(ICOUNT.LT.NOOX) GO TO 152

  161 NPP=0
**    IF THERE IS A FILE CONTAINING A P-T CURVE, TO BE USED FOR A
**      BOILING PATH, PFLUID WILL BE CALCULATED FROM THE
**      PRESSURE CURVE GIVEN IN INP3.
**      THE CORRESPONDING PRESSURE WILL BE INTERPOLATED BETWEEN THE
**      TWO CLOSEST GIVEN T-P POINTS.
*     ** AT PRESENT THERE IS NO SUCH INPUT FILE FOR GASWORKS **
*     THERE IS ANOTHER WAY TO GET A P-T CURVE AFTER STATEMENT 160
**    DO 60 I=1,200
**    READ(INP3,835,ERR=65) TS(I),PS(I)
*     *********************
**    NPP=NPP+1
** 60 CONTINUE

   65 CONTINUE
      BSTP1 = BSTOT + 1
* *********************************************************************
      WRITE(IOUT1,999)
*
*     STARTS READING GASTHERM DATA FILE. INDIVIDUAL log K'S ARE
*       SKIPPED (SECOND ENTRY FOR GASES AND MINERALS,
*       FIRST ENTRY FOR OTHERS.) THESE ARE READ AS DUM OR SKIPPED
*       BY A SLASH IN READ FORMATS.
*     THE FOLLOWING SKIPS THE COMMENTS AT THE TOP OF THE GASTHERM
*       DATA FILE UNTIL A LINE OF STARS IS REACHED (8 STARS MINIMUM,
*       FROM 1ST COL.)

   62 READ(INP2,822) TEXT
      IF(TEXT.NE.ETOIL) GO TO 62

**    DO 160 J = 1, 9
**    IF(J.LT.9) READ(INP2,822) DUM
**    READ(INP2,817) (XPFC(J,I),I=1,5)
**    DO 162 I = 1, 5
**162 COE(I) = XPFC(J,I)
**    CALL POWER(COE,temp,AKLOG)
*     **********
**    XF(J) = AKLOG
**160 CONTINUE
*     THE LINES BELOW FORTRAN LINE 62 ARE PRESENTLY DISABLED, BUT CAN
*      BE USED IN THE FUTURE TO FILL ARRAY XF FROM GASTHERM.
*      THE ARRAY XF CONTAINS 11 ENTRIES WHICH ARE USED IN *CHILLER*
*      PRIMARILY FOR DEBYE-HUCKLE PARAMETERS AND TO SPECIFY GAMMAS FOR
*      DISSOLVED CO2.  THESE ARE UNUSED IN THE PRESENT VERSION OF
*      *GASWORKS*, BUT MAY IN THE FUTURE BE USED FOR NON-IDEAL TREATMENT
*      OF GASES.  XF(9) IS DESIGNATED FOR A POWER FUNCTION GIVING
*      P(FLUID) AS A FUNCTION OF T, THOUGH NONE IS PRESENTLY INSTALLED.
*      ALL DATA DESTINED FOR XF SHOULD BE ADDED TO GASTHERM BELOW THE
*      LINE OF STARS BUT BELOW THE LIST OF COMPONENT SPECIES.


**    IF(XF(9).NE.0.) PFLUID=XF(9)
*     IF(NPP.NE.0) CALL PRESS(PFLUID,temp)

      Plog = Dlog10(PFLUID)
      DO 110 NJ = 1, BSTOT
*
*     DO 110 LOOP READS COMPONENT SPECIES FROM GASTHERM, SKIPPING THOSE
*       NOT SPECIFIED BY SAQ, ABOVE.

      IF(QAS(NJ).NE.0)GO TO 103
      READ(INP2,820)

      GO TO 110
  103 CONTINUE
      NS = QAS(NJ)
      READ(INP2,822) NAME(NS), W(NS)

*     W CONTAINS THE MOLECULAR WEIGHT OF THE BASIS SPECIES

  110 CONTINUE
      NS = SAQTOT

*     READS PROPERTIES OF DERIVED SPECIES FROM GASTHERM.  NAME IS
*     AS ABOVE.
  113 CONTINUE
      NS = NS + 1
      IF(NS.LE.MAXNS) GO TO 114
        NS = NS-1
        WRITE(IOUT1,660) MAXNS, NAME(NS), NS
  660   FORMAT(/5X,'NUMBER OF SPECIES READ IN GASTHERM EXCEEDS MAXIMUM'
     $  ,' OF ',I3, /5X,'LAST NAME READ IS:  ',A8,' (NS=',I3,')'//)
        STOP

  114 READ(INP2,820) NAME(NS),SUP,ITOT,(COEF(NS,I),
     1 SPEC(NS,I), I = 1,ITOT)
*     SUP- SEE COMMENTS BELOW, SUP = NOT BLANK SUPPRESSES SPECIES.
*     ITOT CONTAINS NUMBER OF COMPONENT SPECIES IN THE DERIVED SPECIES,
*       PLUS ONE.
*     COEF(NS,I) CONTAINS STOICHIOMETRIC COEFFICIENT FOR SPECIES "SPEC
*       (NS,I)" IN THE DERIVED SPECIES.
*     SPEC(NS,I) IDENTIFIES COMPONENT SPECIES IN THE DERIVED SPECIES.
*
*     THE I = 1 POSITION IS LEFT BLANK.
*     IF H2 IS IN THE DERIVED GAS, IT SHOULD BE IN THE I=2 POSITION.

*       TEST FOR END OF DERIVED SPECIES LIST BY BLANK LINE.
      IF(NAME(NS).NE.BLANK) GO TO 115
      NST = NS - 1
      GO TO 120
  115 CONTINUE
      IT(NS) = ITOT
      IF(ITOT.LE.MAXSTO) GO TO 112
        WRITE(IOUT1,662) NAME(NS), ITOT, MAXSTO
  662   FORMAT(/5X,'SPECIES ',A8,' IN GASTHERM HAS TOO MANY COMPONENTS'
     $  /5X,'ITOT FOR THAT SPECIES IS ',I3,' ; MAXIMUM IS ',I3//)
        STOP

*     THE I = 1 POSITION IS LEFT BLANK
*     IF H2 IS IN THE COMPLEX, IT SHOULD BE IN THE I=2 POSITION

  112 READ(INP2,818) (COE(I),I = 1, 5)

*     COE CONTAINS POWER FUNCTION COEFFICIENTS TO GENERATE log K AS F(T)

      DO 121 I = 1,5
  121 APFCL(NS,I) = COE(I)
      CALL POWERK(COE,tempK,AKLOG)
*
*       GASTHERM CONTAINS ONE SET OF COEFFICIENTS TO CALCULATE log K AS
*       F (T) WHERE T IS IN DEGREES KELVIN.  THE ABOVE STATEMENTS
*       STORE THE POWER FUNCTION COEFFICIENTS THAT WILL BE NEEDED LATER
*       IN ARRAYS APFCL.
*
  119 LGK(NS) = AKLOG
*
*     ALL GAS SPECIES LISTED WITH A NON BLANK CHARACTER IN COL 16 (SUP)
*       WILL BE SUPPRESSED FROM THE CALCULATION, DISPLAYED ON THE SCREEN
*       AND PRINTED ON THE OUTPUT.
*
      IF(SUP.EQ.BLNK) GO TO 117
      WRITE(IOUT2,885) NAME(NS)
      WRITE(IOUT1,885) NAME(NS)
      GO TO 114
  117 CONTINUE
*
*     WE ALSO SUPPRESS THE GASES MENTIONED IN THE INPUT FILE
      DO 130 I = 1,NSUP
      IF(SUPNAM(I).NE.NAME(NS)) GO TO 130
      WRITE(IOUT2,886) NAME(NS)
      WRITE(IOUT1,886) NAME(NS)
      GO TO 114
  130 CONTINUE
*
*     SORT OUT GASES CONTAINING COMPONENTS BEING USED
*       AND REDEFINE GAS COMPONENT INDICIES,
*     ALSO SUMS UP ALL THE GAS SPECIES COEFFICIENTS (COSUM)
*
      COSUM = 0.
      DO 118 I = 2, ITOT
      JN = SPEC(NS,I)
      IF(QAS(JN).EQ.0) GO TO 114
      SPEC(NS,I) = QAS(JN)
      COSUM = COSUM + COEF(NS,I)
  118 CONTINUE
      COEF(NS,1) = -1.
      SPEC(NS,1) = NS
      GAMMA(NS) = 1.0
      AZERO(NS) = COSUM - 1.
      GO TO 113
*
  120 CONTINUE
      read (inp2,*)
      DO 126 I = 1, SAQTOT
      AQM(I)=MTOT(I)
  126 TMIN(I) = 0.
*
**    CALL HELHUC
*     ***********
      NM = 0
      DO 26 I = 1,MAXSSL
   26 SUMSOL(I) = 0.0
**
**    HERE WE READ COEFFICIENTS TO COMPUTE CROSS VIRIAL COEFFICIENTS FOR
**     H2O-CO2-CH4 GAS MIXTURE. B12 IS FOR B(H2O-CO2), B13 FOR B(H2O-CH4
**     B23 FOR B(CO2-CH4), C112 FOR C(H2O-H2O-CO2) ETC. SEE SUBROUTINE
**     PHI. THESE COEFFICIENTS ARE ASSIGNED AS FOLLOW:
**
**       CVB(1) : B12
**       CVB(2) : B13
**       CVB(3) : B23
**       CVC(1) : C112
**       CVC(2) : C122
**       CVC(3) : C113
**       CVC(4) : C133
**       CVC(5) : C223
**       CVC(6) : C233
**       CVC(7) : C123
**
**    AS OF FEB. 86 WE DO NOT HAVE DATA FOR C113, C133 AND C123. THESE
**      VALUES ARE ASSUMED TO BE ZERO. THIS SEEMS TO BE A BETTER
**      ASSUMPTION THAN TO CALCULATE THEM WITH A MIXING RULE (LIKE
**      A GEOMETRIC AVERAGE).
**
**    DO 28 I=1,3
**     READ(INP2,817) (CVBCF(I,J), J=1,3)
**     ***************
** 28  CVB(I) = CVBCF(I,1)/tempK/tempK + CVBCF(I,2)/tempK + CVBCF(I,3)
**    DO 29 I=1,7
**     READ(INP2,817) (CVCCF(I,J), J=1,3)
**     ***************
** 29  CVC(I) = CVCCF(I,1)/tempK/tempK + CVCCF(I,2)/tempK + CVCCF(I,3)
**
**
   33 CONTINUE
      NM = NM + 1
      IF(NM.LE.MAXNM) GO TO 30
        NM = NM-1
        WRITE(IOUT1,664) MAXNM, MIN(NM)
  664   FORMAT(/5X,'NUMBER OF MINERALS READ IN SOLTHERM EXCEEDS MAXIMUM'
     $  ,' OF ',I3,/5X,'LAST NAME READ IS:  ',A8//)
        STOP
*
   30 READ(INP2,854) MIN(NM),SUP,XMOLW(NM),ITOT,(COEFM(NM,I),SPECM(NM,I)
     1 ,I = 1, ITOT), B(NM), JSOLS(NM)
*     ***************
*     READ MINERAL CHARACTERISTICS FROM GASTHERM.
*     IF SUP NOT BLANK OR NOT 'G' = SUPPRESS MINERAL.  SUP = 'G' FOR
*      GASES (WILL CAUSE CHILLER TO COMPUTE AND OUTPUT FUGACITY).
*     XMOLW = MOLECULAR WEIGHT.
*     COEFM & SPECM ARE AS FOR DERIVED SPECIES, ABOVE.
*     B(NM) CONTAINS EXPONENT FOR MINERAL END MEMBER TO BE USED TO
*       CALCULATE ACTIVITY USING IDEAL MULTI-SITE MIXING.
*     JSOLS(NM) IDENTIFIES THE PARTICULAR SOLID SOLUTION, IF ANY, THAT
*       THIS MINERAL BELONGS TO.
      ITM(NM) = ITOT
*     BLANK LINE IN GASTHERM ENDS READING OF MINERALS.
      IF(MIN(NM).EQ.BLANK) GO TO 165
*
      IF(ITOT.LE.MAXSTO) GO TO 34
        WRITE(IOUT1,666) MIN(NM), ITOT, MAXSTO
  666   FORMAT(/5X,'MINERAL ',A8,' IN GASTHERM HAS TOO MANY COMPONENTS'
     $  /5X,'ITOT FOR THAT MINERAL IS ',I3,' ; MAXIMUM IS ',I3//)
        STOP
*
*    SPECIES  LISTED WITH CHARACTER OTHER THAN BLANK IN COL 9
*     ARE SUPPRESSED FROM THE CALCULATION. THESE MINERALS ARE DISPLAYED
*     ON THE SCREEN AND PRINTED ON OUTPUT.
*
   34 IF(SUP.EQ.BLNK) GO TO 35
      WRITE(IOUT2,885) MIN(NM)
      WRITE(IOUT1,885) MIN(NM)
      READ(INP2,815)
      GO TO 30
   35 CONTINUE
*     WE ALSO SUPPRESS THE SPECIES MENTIONED IN THE INPUT FILE
      DO 135 I = 1,NSUP
      IF(SUPNAM(I).NE.MIN(NM)) GO TO 135
      WRITE(IOUT2,886) MIN(NM)
      WRITE(IOUT1,886) MIN(NM)
      READ(INP2,815)
      GO TO 30
  135 CONTINUE
*
*     DETERMINES IF MINERAL WAS GIVEN IN SATURATED MINERAL LIST
*     AND SET SEQ NUMBERS
*     SEQ CONTAINS THE SEQUENCE POSITIONS OF MINERALS CURRENTLY
*       INCLUDED IN EQUILIBRIUM ASSEMBLAGE.  THESE INTEGERS ARE NOT
*       SIMPLY SEQUENCE POSITIONS IN SOLTHERM, BUT SEQUENCE POSITIONS
*       IN THE RELEVANT MINERAL LIST (AS SORTED BY GASWORKS FROM
*       GASTHERM ACCORDING TO THE COMPOSITION OF THIS PARTICULAR RUN).
      IF(NMT.EQ.0) GO TO 38
      DO 37 I = 1, NMT
      IF(MIN(NM).NE.SATMIN(I)) GO TO 37
      SEQ(I) = NM
      GO TO 38
   37 CONTINUE
*
   38 SCALE(NM)=0.
      DO 40 I = 1, ITOT
      JS = SPECM(NM,I)
*     SKIP IF STOICHIOMETRY INCLUDES SPECIES THAT WERE NOT GIVEN
      IF(QAS(JS).NE.0) GO TO 42
      READ(INP2,815)
      GO TO 30
   42 SPECM(NM,I) = QAS(JS)
      SCALE(NM)=SCALE(NM)+DABS(COEFM(NM,I))
*
*    SCALE(NM) IS USED IN MINSAT TO SCALE THE log Q/K VALUES FOR MINERAL
*      ACCORDING TO THE NUMBER OF COMPONENT SPECIES IN THE MINERAL.
*
   40 CONTINUE
*
      READ(INP2,818) (COE(I),I = 1,5)
*     *********************
*     COE CONTAINS POWER FUNCTION COEFFICIENTS FOR GENERATING log K=F(T)
*
      DO 41 I = 1,5
   41 BPFCL(NM,I)=COE(I)
      CALL POWERK(COE,tempK,AKLOG)
*     **********
*     BPFCL STORES POWER FUNCTION COEFFICIENTS FOR LATER RE-USE
*     THE ABOVE SEQUENCE OF GO TO STATEMENTS DO FOR MINERALS WHAT LINES
*     121 AND 122 DO FOR GASES ABOVE.
*
   45 CONTINUE
      IF(LIMSOL.EQ.0) JSOLS(NM) = 0
*    TURNS OFF SOLID SOLUTIONS IF ALL ARE TURNED OFF BY LIMSOL = 0.
*
      LGKM(NM) = AKLOG
      J = 0
      GO TO 33
  165 CONTINUE
      IF(IDEAL.GT.2) IDEAL = 2
*
      NMTT = NM - 1
      IF(NMT.EQ.0) GO TO 175
      DO 48 I = 1, NMT
*
*    DO 48 SORTS OUT MINERALS CURRENTLY INCLUDED, AS SPECIFIED IN
*     WORKRUN WHERE MINERAL NAMES ARE INPUT
*
      NM = SEQ(I)
      IF(NM.NE.0) GO TO 49
       WRITE(IOUT1,656) SATMIN(I)
       WRITE(IOUT2,656) SATMIN(I)
  656  FORMAT(/3X,A8,', SPECIFIED IN WORKRUN,' ,
     $ ' WAS NOT FOUND IN GASTHERM. CHECK SPELLING.'//)
       STOP
   49 QES(NM)=I
      MINTRY(NM)=ACT(I)
*
*    SETS TRIAL MOLES OF MINERALS AS INPUT IN ACT ARRAY IN CHILLRUN.
*
   48 CONTINUE
*
  175 CONTINUE
*
      BSTOT = SAQTOT
      BSTP1 = SAQTOT + 1
      JSTOP = 0
*
*
*     ******************************************************************
*     END INPUT    BEGIN DATA ECHO
*     ******************************************************************
*     ******************************************************************
*
*
      WRITE(IOUT1,970)
      WRITE(IOUT1,943)
      WRITE(IOUT1,960) TITLE1
      WRITE(IOUT1,960)  TITLE2
      WRITE(IOUT1,943)
      WRITE(IOUT1,970)
      WRITE(IOUT1,923)
      WRITE(IOUT1,940)
      WRITE(IOUT1,945) (NAME(I),I,MTOT(I), MTRY(I), GAMMA(I),
     1 COMTOT(I), I=1,BSTOT)
*
   50 WRITE(IOUT1,611) ERPC, FO2log, PFLUID, SLIM, temp,SINC,tempC,
     1 TOTMIX, SOLMIN, RM, GASGRM,SUPRNT,GASTOT
      WRITE(IOUT2,611) ERPC, FO2log, PFLUID, SLIM, temp,SINC,tempC,
     1  TOTMIX, SOLMIN, RM,GASGRM,SUPRNT,GASTOT
  611 FORMAT(//4X,'THE FOLLOWING PARAMETERS WERE GIVEN:'
     1 /6X,'ERPC  = ',D10.4,'  FO2log = ',F10.5,'  PFLUID = ',F10.5,
     2 /6X,'SLIM  = ',G10.4,'  temp   = ',F10.5,'  SINC   = ',G10.4,
     3 /6X,'tempC = ',F10.5,'  TOTMIX = ',G10.4,'  SOLMIN = ',D10.4,
     4 /6X,'RM    = ',F10.5,'  GASGRM = ',F10.5,'  SUPRNT = ',D10.4,
     5 /6X,'GASTOT= ',F10.5)
      WRITE(IOUT1,612)  IFRAC, IPUNCH, NLOOP, ISTEP
     1, LIMSS, LOOC, ITREF, IDEAL,INCREM,INCP,MINSOL,OFLAG
      WRITE(IOUT2,612)  IFRAC, IPUNCH, NLOOP, ISTEP
     1, LIMSS, LOOC, ITREF, IDEAL, INCREM,INCP,MINSOL,OFLAG
  612 FORMAT(//4X,'THE FOLLOWING OPTIONS WERE SELECTED:'
     1 /6X,'IFRAC =',I5,'  IPUNCH =',I5,'  NLOOP  =',I5,'  ISTEP =',I5,
     2 /6X,'LIMSOL=',I5,'  LOOC   =',I5,'  ITREF  =',I5,'  IDEAL =',I5,
     3 /6X,'INCREM=',I5,'  INCP   =',I5,'  MINSOL =',I5,'  OFLAG =',I5)
      IF(NPP.NE.0) WRITE(IOUT1,614)  NPP
  614 FORMAT(6X,'GIVEN PRESSURE GRID. ',I3,' T-P POINTS READ'/)
      WRITE(IOUT1,923)
      WRITE(IOUT1,617) SUPRNT
  617 FORMAT(2X,'OF THE DERIVED SPECIES LISTED BELOW, ONLY THOSE WITH'
     1 ,/2X,'MOLE NUMBERS ABOVE ', D10.4,' WILL APPEAR ON THE OUTPUT'/)
      WRITE(IOUT1,950)
*     DOUBLE COLUMN ECHO OF DERIVED SPECIES AND MINERALS FROM CHILLER
*     IMPORTED HERE
      NSTMBS = NST - BSTOT
      K = NSTMBS/2
      IF(MOD(NSTMBS,2).NE.0) K = K + 1
      WRITE(IOUT1,955) (NAME(I), I, LGK(I), NAME(I+K),
     1 I+K, LGK(I+K), I = BSTP1,BSTOT+K)
      IF(NMTT.EQ.0) GO TO 176
      WRITE(IOUT1,900)
      L = NMTT/2
      IF(MOD(NMTT,2).NE.0) L = (NMTT+1)/2
      WRITE(IOUT1,910) (MIN(I), I, LGKM(I), MINTRY(I), MIN(I+L),
     1 I+L, LGKM(I+L), MINTRY(I+L), I = 1,L)
*
      IF(MINSOL.EQ.0) GO TO 176
      IF(MINSOL.LE.2) WRITE(IOUT1,974)
      IF(MINSOL.GT.2) WRITE(IOUT1,976)
      WRITE(IOUT1,975) (NOMOX(I), WTPC(I), I = 1, NOOX)
      GRMLG = -100.
      IF(TOTMIX .GT.0.) GRMLG = Dlog10(TOTMIX)
      WRITE(IOUT1,971) TOTMIX, TYP, GRMLG
      WR = 0.
      IF(TOTMIX.NE.0.) WR = GASGRM/TOTMIX
      WRITE(IOUT1,973) GASGRM, WR
*     ******************************************************************
*
  176 CONTINUE
      IF(JSTOP.EQ.1) GO TO 388
*
      CALL SQUOMP
*     ***********
*
      IF(NLOOP.NE.0) GO TO 170
      WRITE(IOUT2, 616)
      WRITE(IOUT1, 616)
  616 FORMAT(/6X,'END DATA ECHO. SET NLOOP > 0 TO START CALCULATIONS'/)
      STOP
  170 CONTINUE
*
**    TO EASE CONVERGENCE W/ PICKUP RUN, WE CALCULATE GAMMAS SOLID SOL.
**    NM = SOLS(10,1)
**    JS = ISOLS(NM)
*     ********
  189 IF(INCREM.NE.0) GO TO 560
*
*     STEP STARTS HERE
*     ===============================================================
*
  190 CONTINUE
*     N IS THE TOTAL NUMBER OF UNKNOWNS
  191 CONTINUE
      ERLOOP = 0
      NDUMP= 0
      MINCT= 0
      JSTOP= 0
*     GASTOT EQUATION IS IN THE BSTOT + 1 POSITION IN THE MATRIX
      JGSTOT = BSTOT + 1
      N = BSTOT + NMT  + 1
        IF(N.LE.MAXN) GO TO 192
        WRITE(IOUT1,668) N, MAXN
  668   FORMAT(/5X,'SIZE OF MATRIX REACHED LIMIT.  N = ',I3,'  MAXIMUM'
     $  ,' SIZE IS ',I3//)
        STOP
  192   IF(NMT.LE.MAXIN) GO TO 193
        WRITE(IOUT1,670) NMT, MAXIN
  670   FORMAT(/5X,'TOO MANY PHASES IN MATRIX (',I3,'); MAXIMUM'
     $  ,' NUMBER ALLOWED IS ',I3//)
        STOP
*
*      THE PRECEDING "IF" SETS AN APPROXIMATE FH2, ASSUMING THAT H20 IS
*        THE OVERWHELMINGLY DOMINANT GAS
*     HERE, WE CALCULATE CONVERGENCE LIMITS FOR F'S BELOW
*     ERPC SHOULD BE D-12 OR LARGER (32 BIT MACHINE)
  193 DO 194 I = 1,N
      IF (I.EQ.JGSTOT) THEN
         CONV(I)=ERPC*GSTOLD
      ELSE IF (I.GT.JGSTOT) THEN
         II = I-JGSTOT
         NM = SEQ(II)
         CONV(I) = ERPC*DABS(LGKM(NM))
         IF(DABS(LGKM(NM)).EQ.0.) CONV(I) = ERPC*10.
      ELSE IF (I.EQ.1.AND.FO2.NE.0.) THEN
      CONV(1) = ERPC*DABS(LGK(BSTP1))
      ELSE IF (DABS(MTRY(I)).GT.DABS(MTOT(I))) THEN
         CONV(I) = ERPC*DABS(MTRY(I))
      ELSE IF (MTOT(I).EQ.0.) THEN
         CONV(I) = ERPC
      ELSE
         CONV(I) = ERPC*DABS(MTOT(I))
      ENDIF
  194 CONTINUE
*
*     ******************************************************************
*
*      HERE WE GO TO 313 IN FIRST LOOP TO COMPUTE MTRYS OF DERIVED
*      SPECIES  BEFORE SETTING THE MATRIX
       GO TO 313
*
*      >>>>>   ITERATIONS START HERE  <<<<<
*      ==================================================
*
  200 CONTINUE
*
  202 ERLOOP = ERLOOP + 1
      NSTOP = 0
      IF(FO2log.EQ.0.) GO TO 210
*     ABOVE IF SKIPS SUMMATION OF H2-BEARING SPECIES TO GET MTOT(H2) IF
*     FO2 IS INPUT
      IF(ERLOOP.EQ.1) GO TO 225
      H2TOT = MTRY(1)
*
*     SUMS H2 SPECIES FOR MTOT OF H2 IF FO2 IS INPUT
      DO 204 NS = BSTP1, NST
      IF(SPEC(NS,2).NE.1) GO TO 204
      H2TOT = H2TOT + MTRY(NS)*COEF(NS,2)
  204 CONTINUE
      H2TOT = H2TOT + TMIN(1)
*
      MTOT(1) = H2TOT
  210 CONTINUE
  225 CONTINUE
*
*     ******************************************************************
*
      DO 250 I = 1,N
      DO 245 J = 1,N
      A(I,J) = 0.0
  245 CONTINUE
      F(I) = 0.
      A(I,I) = 1.
  250 CONTINUE
*     INITIALIZE THE MATRIX VALUES FOR THE GASTOT EQUATION
      A(JGSTOT,JGSTOT)=-1.
      F(JGSTOT)=-GASTOT
      DO 300 KK = 1, BSTOT
      F(JGSTOT) = F(JGSTOT) + MTRY(KK)
      A(JGSTOT,KK) = 1.
*
*
*     THE FOLLOWING IF'S CONTROL PRODUCTION OF MATRIX ENTRIES
*     FOR O2
      IF(KK.EQ.1.AND.FO2log.NE.0.) GO TO 300
*
*
*
      F(KK) = -MTOT(KK) + MTRY(KK)
      DO 290 NS = BSTP1, NST
*     THE DO 290 LOOP THROUGH ALL DERIVED SPECIES FOR INCLUSION IN
*     APPROPRIATE MASS BALANCE EQUATIONS.
      ITOT = IT(NS)
      DO 265 IC = 2,ITOT
      IF(SPEC(NS,IC).EQ.KK) GO TO 267
  265 CONTINUE
      GO TO 290
  267 CONTINUE
      CFF = COEF(NS,IC)
      A(KK,JGSTOT) = A(KK,JGSTOT) - CFF*MTRY(NS)*AZERO(NS)/GASTOT
      A(JGSTOT,KK) = A(JGSTOT,KK) + CFF*MTRY(NS)/MTRY(KK)
      DO 288 II = 2, ITOT
      J = SPEC(NS,II)
*     DTERM IS THE DERIVATIVE OF MOLALITY EXPRESSION FOR DERIVED SPECIES
*     WITH RESPECT TO COMPONENT SPECIES J (= SPEC(NS,II)
      DTERM = COEF(NS,II) * MTRY(NS)/MTRY(J)
      A(KK,J) = A(KK,J) + CFF*DTERM
  288 CONTINUE
      F(KK) = F(KK) + CFF*MTRY(NS)
  290 CONTINUE
  300 CONTINUE
      DO 301 NS = BSTP1,NST
      F(JGSTOT) = F(JGSTOT) + MTRY(NS)
      A(JGSTOT,JGSTOT) = A(JGSTOT,JGSTOT) - AZERO(NS)*
     1 MTRY(NS)/GASTOT
  301 CONTINUE
*     MINERA WILL SET UP MASS ACTION EQS. FOR MINERALS
      IF(NMT.NE.0) CALL MINERA
*                  ***********
      IF(FO2log.NE.0.) CALL OXY
*                  ********
*     PROCEDURE TO HELP CONVERGENCE (SUPPOSEDLY !!). SKIP IF RM=0 OR 1
*     TRY ONLY IF THERE ARE CONVERGENCE PROBLEMS (LAST RESOURCE).
*     ASSURES THAT ALL F'S DECREASE AFTER EACH ITERATION BY MULTIPLYING
*     F'S BY POWERS OF RM (GENERALLY RM IS BETWEEN .1 AND .9).
      IF(RM.EQ.0..OR.RM.GE.1.) GO TO 297
      DO 292 I = 1,N
      FFF = F(I)
        IF(DABS(F(I)).LT.DABS(FOLD(I))) GO TO 296
         FACT(I) = FACT(I) * RM
         F(I) = FFF * FACT(I)
         IF(FACT(I).LT.1.D-4) FACT(I) = 1.D-4
        GO TO 294
  296  FACT(I) = 1.0
  294  FOLD(I) = FFF
  292 CONTINUE
*
*
  297 IF(ERLOOP.LE.5) GO TO 365
*     CONVERGENCE CRITERION.  SEES THAT ALL F(I) <  CONV(I) DEFINED
*     AFTER STATEMENT 191
      DIFTOT = GASTOT
      DO 360 NS = 1, NST
      DIFTOT = DIFTOT - MTRY(NS)
  360 CONTINUE
      IF(DABS(DIFTOT).GT.ERPC) GO TO 365
      DO 363 I = 1,N
      IF(DABS(F(I)).GT.CONV(I)) GO TO 365
  363 CONTINUE
      GO TO 410
*
  365 IF(ISTEP.EQ.0) GO TO 308
      IF(ERLOOP.NE.2) GO TO 308
  306 NSTOP = 1
      GO TO 388
  307 NSTOP = 0
      NDUMP = 1
      CALL DUMP
*     *********
  308 CONTINUE
      IF(ERLOOP.EQ.1.AND.ISTEP.EQ.1) CALL DUMP
      CALL SLUD
*     *********
      CALL SIR(X)
*     ********
      IF(ERLOOP.EQ.1.AND.ISTEP.EQ.1) THEN
      WRITE(IOUT1,977) (X(I), I = 1, N)
      ENDIF
      IF(NDUMP.EQ.1.AND.ISTEP.EQ.1) CALL DUMP2
*                                   **********
      NDUMP = 0
      FMAX =  PFLUID * GSTOLD
      PlogP1 = Dlog10(FMAX)
      DO 310 KK = 1, BSTOT
      COLD = MTRY(KK)
      MTRY(KK) = MTRY(KK) + X(KK)
*     LIMIT SIZE OF MTRY VALUES TO AID CONVERGENCE
      IF(MTRY(KK).GT.FMAX)  MTRY(KK) = COLD
      XX = X(KK)
  309 CONTINUE
      IF(MTRY(KK).GT.0.) GO TO 320
      XX = XX/2.
      MTRY(KK) = COLD + XX
      GO TO 309
  320 CONTINUE
      IF(MTRY(KK).GT.RMIN)  GO TO 310
      WRITE(IOUT2,676) NAME(KK), RMIN
  676 FORMAT(5X,'** WARNING !  MOLALITY OF ',A8,' LESS THAN ', D10.3)
      MTRY(KK) = RMIN
  310 CONTINUE
*     BRACKETS GASTOT TO EASE CONVERGENCE
      GASTOT = GASTOT + X(JGSTOT)
      GASMAX = GSTOLD*10.
      GASMIN = GSTOLD/10.
      IF (GASTOT.GT.GASMAX) GASTOT = GASMAX
      IF (GASTOT.LT.GASMIN) GASTOT = GASMIN
      NNMT = NMT
      IF(NMT.EQ.0) GO TO 313
      IF(PFLUID.GT.PMAX) PFLUID = PMAX
*
      Plog = Dlog10(PFLUID)
      DO 414 I = 1,MAXSSL
  414 SOLTRY(I) = 0.0
      IQ = 0
      DO 314 I = 1, NMT
      IQ = IQ + 1
      NM = SEQ(IQ)
      M = BSTOT + I + 1
      MINTRY(NM) = MINTRY(NM) + X(M)
*
      IF(ISOLS(NM).EQ.0) GO TO 417
      JU = ISOLS(NM)
      IF(MINTRY(NM).LE.0.) GO TO 311
       SOLTRY(JU) = SOLTRY(JU) + MINTRY(NM)
       GO TO 314
*     WE KICK OUT SOLSOL IF LE 0 AND LIMSS = 2
  311 IF(ERLOOP.GT.6.AND.LIMSS.GT.1) GO TO 420
*      WE RESET NEGATIVE SOL SOL TO POSITIVE IF NO KICK OUT
       MINTRY(NM) = 1.D-20
       SOLTRY(JU) = SOLTRY(JU) + MINTRY(NM)
       GO TO 314
*
*     PURE MINERAL: WE KICK OUT MINERAL IF MINTRY LT SOLMIN
  417 IF(ERLOOP.LE.4.OR.SOLMIN.GE.0.) GO TO 314
      IF(MINTRY(NM).GE.SOLMIN) GO TO 314
*
*     MINERAL KICK OUT SCHEME BELOW
  420 WRITE(IOUT2,777) MIN(NM), MINTRY(NM)
      WRITE(IOUT1,777) MIN(NM), MINTRY(NM)
  777 FORMAT(5X,A8,' MINTRY= ',D10.3,' MINERAL REMOVED FROM MATRIX')
      SEQ(IQ) = 0
      QES(NM) = 0
      MINTRY(NM) = 1.D-300
*     SORT QES AND SEQ AND CONV WITHOUT CURRENT NM MINERAL
      NNMT = NNMT - 1
*       IF NO MINERALS ARE LEFT IN MATRIX, WE SKIP RESORTING
      IF(NNMT.EQ.0) GO TO 314
       DO 418 L = IQ,NNMT
          LL = L+1
          J = BSTOT + L + 1
          JJ = J + 1
          SEQ(L) = SEQ(LL)
          CONV(J) = CONV(JJ)
          F(J) = F(JJ)
  418  QES(SEQ(L)) = L
      IQ = IQ - 1
      IF(ISOLS(NM).EQ.0) GO TO 314
*
*     IF WE KICKED OUT A SOLID SOLUTION ENDMEMBER
      ISOLS(NM) = 0
      JS = JSOLS(NM)
      IS = IST(JS)
      ISTOT = 0
*     WE COUNT HOW MANY ENDMEMBERS LEFT IN SOLID SOLUTION (ISTOT)
       DO 336 IL = 1,IS
          LS = SOLS(JS,IL)
          IF(QES(LS).EQ.0) GO TO 336
          ISTOT = ISTOT + 1
          LSBR = LS
  336  CONTINUE
*
*      IF ONLY ONE END MEMBER LEFT, WE BREAK SOLID SOLUTION
*      BY SETTING ISOLS TO 0
       IF(ISTOT.NE.1) GO TO 314
       WRITE(IOUT1,778) MIN(LSBR)
       WRITE(IOUT2,778) MIN(LSBR)
  778  FORMAT(5X,A8,' NOT IN SOLID SOLUTION ANYMORE')
       ISOLS(LSBR) = 0
*
  314 CONTINUE
*
*
*     NMT AND N HAVE TO BE RESET IN CASE WE KICKED OUT MINERALS
      NMT = NNMT
      N = NNMT + BSTOT + 1
  313 CONTINUE
*
**    HERE WE COMPUTE ACTIVITY COEFFICIENTS FOR
**    NON IDEAL SOLID SOLUTIONS THAT ARE SATURATED (CURRENTLY IN
**    MATRIX).  WE HAVE ONLY GOLD SO FAR  SOL.SOL ID NO. 10
**
**    NM = SOLS(10,1)
**    JS = ISOLS(NM)
**    IF(JS.NE.0) CALL ELECTR
**    ***********
      BSTP2 = BSTP1
      JK = 1
      IF(FO2.EQ.0.) GO TO 315
      MTRY(BSTP1) = FO2 * GASTOT/(GAMMA(BSTP1) * PFLUID)
      ACT(BSTP1) = GAMMA(BSTP1) * MTRY(BSTP1) * PFLUID/GASTOT
      ACT(1) = MTRY(1)*GAMMA(1)*PFLUID/GASTOT
      JK = 2
      BSTP2 = BSTP1 + 1
  315 DO 319 NS=BSTP2,NST
      CSAVE= -Dlog10(GAMMA(NS)) - LGK(NS)
      ITOT = IT(NS)
      DO 321 I = 2, ITOT
      J = SPEC(NS,I)
      CSAVE2 = (Dlog10(GAMMA(J) * MTRY(J))) * COEF(NS,I)
      CSAVE= CSAVE + CSAVE2
  321 CONTINUE
      PTERM=AZERO(NS) * Dlog10(PFLUID/GASTOT)
      RSL=CSAVE + PTERM
      ATERM = RSL + Dlog10(GAMMA(NS) * PFLUID/GASTOT)
      IF(RSL.LT.PlogP1) GO TO 316
      RSL = PlogP1
      ATERM = PlogP1
  316 IF(RSL.GT.RlogMN) GO TO 317
      RSL = RlogMN
      ATERM = RlogMN
  317 MTRY(NS)=10.**RSL
      ACT(NS) = 10.**ATERM
  319 CONTINUE
      IF(ERLOOP.EQ.0) GO TO 200
*
*     ******************************************************************
      WRITE(IOUT2,650) ERLOOP, GASTOT
  650 FORMAT(5X,'LOOP NO.',I4,'  GASTOT: ',D15.7)
      IF(ERLOOP.EQ.NLOOP) GO TO 410
*     COMPUTE THE FUGACITIES (ACT) OF THE COMPONENT SPECIES
  352 DO 358 NS = JK, BSTOT
      FTERM = DLOG10(MTRY(NS))
      ATERM = FTERM + DLOG10(GAMMA(NS) * PFLUID/GASTOT)
      ACT(NS) = 10.**ATERM
  358 CONTINUE
      IF(ISTEP.EQ.1) WRITE(IOUT1,905) ERLOOP,GASTOT
      IF(ISTEP.EQ.0) GO TO 366
      DO 362 NS = 1, BSTOT
  362 WRITE(IOUT1,920) NS, NAME(NS), MTRY(NS)
      WRITE(IOUT1,977) (X(I), I = 1, N)
  366 CONTINUE
*     WE COMPUTE AQM, THE QUANTITY OF GASEOUS COMPONENTS
      DO 350 I=1,SAQTOT
      AQM(I) = MTRY(I)
      DO 370 NS = BSTP1, NST
      ITOT = IT(NS)
      DO 364 IC = 2,ITOT
      IF(SPEC(NS,IC).EQ.I) GO TO 367
  364 CONTINUE
      GO TO 370
  367 CONTINUE
      CFF = COEF(NS,IC)
      AQM(I) = AQM(I) + CFF*MTRY(NS)
  370 CONTINUE
  350 CONTINUE
      GO TO 200
*
*
  410 CONTINUE
      GSTOLD = GASTOT
*     SUM H2-BEARING SPECIES IF FO2 IS INPUT
      IF (FO2log.EQ.0.) GO TO 412
      H2TOT = MTRY(1)
      ACT(1) = MTRY(1)*GAMMA(1)*PFLUID/GASTOT
      DO 411 NS=BSTP1,NST
      IF(SPEC(NS,2).NE.1) GO TO 411
      H2TOT = H2TOT + COEF(NS,2)*MTRY(NS)
  411 CONTINUE
      H2TOT = H2TOT + TMIN(1)
      MTOT(1) = H2TOT
  412 JSTOP = 1
  388 CONTINUE
      WRITE(IOUT1,999)
      WRITE(IOUT1,960) TITLE1
      WRITE(IOUT1,960) TITLE2
      WRITE(IOUT1,943)
      WRITE(IOUT1,970)
      WRITE(IOUT1,946) ERLOOP,NLOOP
      WRITE(IOUT1,943)
      WRITE(IOUT1,933) temp, PFLUID, TOTMIX
      WRITE(IOUT1,927) GASTOT
      WRITE(IOUT1,944)
      WRITE(IOUT1,914)
      WRITE(IOUT1,915)
      GRT=0.
      DO 400 NS = 1, NST
      IF(NS.LE.BSTOT) GRT=GRT+ AQM(NS)*W(NS)
      IF(NS.GT.BSTOT) GO TO 396
      ALA = Dlog10(ACT(NS))
**    ALGAM = Dlog10(GAMMA(NS))
*     ALM IS THE PARTIAL PRESSURE
      ALM = ACT(NS)/GAMMA(NS)
*     ALMF IS THE MOLE FRACTION
      ALMF = ALM/PFLUID
      WRITE(IOUT1,920) NS,NAME(NS),MTRY(NS),ALMF,ALM,ACT(NS),ALA,
     1 GAMMA(NS)
      GO TO 400
  396 IF(MTRY(NS).LT.SUPRNT.AND.NS.NE.BSTP1) GO TO 400
      ALA = Dlog10(ACT(NS))
**    ALGAM = Dlog10(GAMMA(NS))
      ALM = ACT(NS)/GAMMA(NS)
      ALMF = ALM/PFLUID
      WRITE(IOUT1,920) NS,NAME(NS),MTRY(NS),ALMF,ALM,ACT(NS),ALA,
     1 GAMMA(NS)
  400 CONTINUE
      GRAMSO = GRT
*
  402 WRITE(IOUT1,923)
      GRT = 0.
      IF(NMT.EQ.0) GO TO 398

      WRITE(IOUT1,957)
      DO 404 NS = 1, BSTOT
      TMOLAL = AQM(NS)/GASTOT
      PPM=1.D06 * AQM(NS) * W(NS)/GRAMSO
*     GRAMSO CONTAINS GRT, WHICH IS TOTAL MASS OF GAS PHASE,
*      THUS PPM IS PPM IN GAS PHASE OF GIVEN COMPONENT GAS SPECIES
*     TMOLAL CONTAINS A BULK MOLE FRACTION (OF SORTS)
      WRITE(IOUT1,962) NS, NAME(NS), MTRY(NS), TMOLAL, PPM
  404 CONTINUE
      WRITE(IOUT1,965) GRAMSO
      IF(GASGRM.EQ.0.) GASGRM = GRAMSO
       WRITE(IOUT1,961)
      DO 405 NS = 1, BSTOT
*     TMIN WAS COMPUTED IN MINERAL. IT IS THE TOTAL NUMBER
*      OF MOLES OF COMPONENTS PRESENT IN SOLIDS
      WRITE(IOUT1,962) NS, NAME(NS), MTOT(NS), AQM(NS), TMIN(NS)
  405 CONTINUE
      WRITE(IOUT1,923)
      WRITE(IOUT1,930)
      WRITE(IOUT2,923)
      WRITE(IOUT2,930)
      GRT = 0.
      DO 427 I = 1, NMT
      NM = SEQ(I)
      GRT = GRT + MINTRY(NM) * XMOLW(NM)
  427 CONTINUE
      VOLMIN = 0.
      DO 395 I = 1, NMT
      NM = SEQ(I)
      GRAMS = MINTRY(NM) * XMOLW(NM)
      ALGG=1.
      ALGM=1.
      IF(MINTRY(NM).LT.0.) GO TO 429
      ALGM = Dlog10(MINTRY(NM))
      ALGG = Dlog10(GRAMS)
  429 CONTINUE
      XMIN = 1.0
      JS = ISOLS(NM)
      IF(JS.NE.0) XMIN = MINTRY(NM)/SOLTRY(JS)
      ACTMIN = GAMIN(NM)*XMIN
      IF(JS.NE.0) ACTMIN = ACTMIN**B(NM)
      LOGACT = Dlog10(ACTMIN)
      IF(GRT.NE.0.) PCMIN = GRAMS/GRT * 100.
**    IF(NM.GT.ITGAS) GO TO 435
**    PI = XMIN*PFLUID
**    FI = Dlog10(PI*GAMIN(NM))
**     CALCULATES VOLUME FROM PV/RT = Z, Z CALCULATED FROM BVAL AND
**     AND CVAL IN SUBROUTINE PHI
**    ZGAS = 1 + BVAL(NM)*PFLUID + CVAL(NM)*PFLUID*PFLUID
**    CM3 = tempK * 83.1432 * ZGAS * MINTRY(NM) / PFLUID
*     WRITES SATURATED MINERALS DATA
  435 IF(DENS(NM).EQ.0.) DENS(NM) = 2.5
      CM3 = GRAMS/DENS(NM)
      VOLMIN = VOLMIN + CM3
      WRITE(IOUT1,935) MIN(NM), MINTRY(NM), ALGM, GRAMS,
     1 ALGG,PCMIN,CM3
      WRITE(BUFOUT(I),937) MIN(NM),XMIN,ACTMIN,GAMIN(NM),LOGACT
*      TERMINAL OVERLOOK
  395 WRITE(IOUT2,935) MIN(NM), MINTRY(NM), ALGM, GRAMS,
     1 ALGG,PCMIN,CM3
      WRITE(IOUT1,931)
      WRITE(IOUT1,'(A80)') (BUFOUT(I),I=1,NMT)
  408 WRITE(IOUT1,943)
*
  398 IF(tempC.NE.0.) WRITE(IOUT1,936) TOTMIX
      WRITE(IOUT1,972) GRT, VOLMIN
      WRITE(IOUT1,968)
      IF(tempC.EQ.0.) GO TO 382
*     COMPUTE WEIGHT OF MIXING SOLUTION
      TOTWT = 0.
      DO 383 I = 1, BSTOT
  383 TOTWT = TOTWT + TOTMIX * COMTOT(I) * W(I)
      IF(TOTMIX.NE.0.) WR = GASGRM/TOTWT
      GO TO 380
*
  382 IF(MINSOL.EQ.0) GO TO 380
      GRMLG = -100.
      IF(TOTMIX.GT.0.) GRMLG=Dlog10(TOTMIX)
      WRITE(IOUT1,971) TOTMIX, TYP, GRMLG
      IF(TOTMIX.NE.0.) WR = GASGRM/TOTMIX
      WRITE(IOUT1,973) GASGRM, WR
  380 IF(NSTOP.EQ.1) GO TO 307
**
**    NON IDEAL SOLID SOLUTIONS. WE NEED TO COMPUTE ACTIVITIES BEFORE
**    CALL TO MINSAT, EVEN IF MINSAT WILL RECOMPUTE THESE ACTIVITIES.
**    THIS STEP IS ONLY FOR PHASES THAT ARE NOT YET SATURATED.
**    WE NEED ALL END MEMBERS ACTIVITIES TO CHECK IF THE SOLID
**    SOLUTION SHOULD BE SUPERSATURATED.
**    SO FAR, WE HAVE ONLY ELECTRUM    SOL-SOL NO. = 10
**    SILVER SHOULD BE ENTERED AFTER GOLD IN SOLTHERM  IS=1 AU, =2 AG
**    NM = SOLS(10,1)
**    LS = IST(10)
**    GOLD NOT IN THE CALCULATION, OR ALREADY SATURATED => SKIP
**    IF(NM.EQ.0.OR.QES(NM).NE.0) GO TO 515
**    DO 500 IS = 1,LS
**    NM = SOLS(10,IS)
**    CSAVE = 0.
**    ITOT=ITM(NM)
**    DO 505 I=1,ITOT
**      CFM = COEFM(NM,I)
**      SPM = SPECM(NM,I)
**      CSAVE1 = CFM*Dlog10(ACT(SPM))
**505 CSAVE = CSAVE + CSAVE1
**    ACTM(IS) = CSAVE - LGKM(NM)
**500 CONTINUE
**    CALL ELSAT(ACTM(1),ACTM(2),tempK)
**    **********
**515 CONTINUE
**    NOW WE NEED TO COMPUTE GAS FUGACITIES, SATURATION PRESSURE AND
**    FUGACTITY COEFFICIENTS, IF GASES ARE NOT SUPERSATURATED,
**    BEFORE WE CALL MINSAT. WE TEST IF GASES ARE SUPERSATURATED
**    (THAT MEANS IN THE MATRIX) ON H2O GAS. IF H2O GAS IS NOT IN
**    THE MATRIX (QES(1)=0) THEN GASES ARE NOT SUPERSATURATED.
**    WE REPEAT SOME OF THE ABOVE W/ SOLID SOL, BUT IT IS CLEANER
**    IF (QES(1).NE.0.OR.IPSAT.EQ.0) GO TO 399
**    DO 438 NM=1,NGAS
**    CSAVE = 0.
**    ITOT=ITM(NM)
**    DO 440 I=1,ITOT
**      CFM = COEFM(NM,I)
**      SPM = SPECM(NM,I)
**      CSAVE1 = CFM*Dlog10(ACT(SPM))
**440 CSAVE = CSAVE + CSAVE1
**    GASlog(NM) = CSAVE - AKGASL(NM)
**438 CONTINUE
**    CALL PSAT
*     *********
*
      IF(LOOC.EQ.3) GO TO 406
      CALL MINSAT
*     ***********
      WRITE(IOUT1,944)
      IF(ERLOOP.LT.NLOOP) GO TO 406
      WRITE(IOUT1,620) NLOOP
      WRITE(IOUT2,620) NLOOP
  620 FORMAT(//5X,'NUMBER OF ITERATIONS REACHED THE GIVEN LIMIT (',I4,
     1    ')',/5X,'NON-CONVERGENCE ASSUMED. EXECUTION ABORTED.'//)
      STOP
*
  406 IF(IPUNCH.EQ.0) GO TO 407
      IF(BIT.NE.SINC.AND.LOOC.NE.3) GO TO 407
*
*     WRITE PLOTTING FILE ON OPTION
*
      IF(IPUNCH.EQ.1.OR.IPUNCH.EQ.2) THEN
*     SYMONDS PLOT FILE
      WRITE(IPUN,832)
      WRlog = 0.
      IF(TOTMIX.NE.0.) WRlog = Dlog10(WR)
      ALA1=Dlog10(ACT(BSTP1))
      DO 409 NS = 1,NST
      ALA = Dlog10(ACT(NS))
*     PLOT log MOLES IF INPUNCH EQUALS 2
      IF(IPUNCH.EQ.2) ALA = Dlog10(MTRY(NS))
      WRITE(IPUN,832) NAME(NS),ALA,temp,PFLUID,ALA1,WRlog
  409 CONTINUE
      DO 413 I = 1, NMTT
*+    XMIN=1.
*+    JS=ISOLS(SEQ(I))
*+    IF(JS.NE.0) XMIN=MINTRY(SEQ(I))/SOLTRY(JS)
*+    GRAMS = MINTRY(SEQ(I)) * XMOLW(SEQ(I))
*+    IF(GRT.NE.0.) PCMIN=GRAMS/GRT*100.
      ALA = MINTRY(I)
      IF(ALA.EQ.0.) ALA = 1.D-300
      ALA = Dlog10(ALA)
      WRITE(IPUN,832) MIN(I),ALA,temp,PFLUID,ALA1,WRlog
  413 CONTINUE
*
      ELSE IF(IPUNCH.EQ.3.OR.IPUNCH.EQ.4) THEN
*
*     New plotting code inserted here, modified from CHILLER.  This will
*      allow GASWORKS to run with MINPLOT plotting program.
*
*     CONTROL WRITING OF TITLE AND NST, BSTOT and space holder (1) AT
*       TOP OF FIRST PLOT SEGMENT
      IF(IPUNCH.EQ.4) GO TO 415
      IPUNCH = 4
      WRITE(IPUN,800) BSTOT,NST,1
      WRITE(IPUN,870) TITLE1
      WRITE(IPUN,870) TITLE2
  415 TOT = TOTMIX
      WRITE(IPUN,833) temp,PFLUID,ACT(BSTP1),TOT
      WRITE(IPUN,830) (MTRY(I), I=1,NST)
      WRITE(IPUN,830) (MTOT(I), I=1,BSTOT)
      WRITE(IPUN,830) (AQM(I), I=1,BSTOT)
      WRITE(IPUN,838) 1,'bogus', 0.
      DO 416 I = 1, NMT
      XMIN=1.
      JS=ISOLS(SEQ(I))
      IF(JS.NE.0) XMIN=MINTRY(SEQ(I))/SOLTRY(JS)
      GRAMS = MINTRY(SEQ(I)) * XMOLW(SEQ(I))
      IF(GRT.NE.0.) PCMIN=GRAMS/GRT*100.
      WRITE(IPUN,837) SEQ(I), MIN(SEQ(I)), MINTRY(SEQ(I)),
     1 GRAMS, PCMIN, XMIN
  416 CONTINUE
      WRITE (IPUN, 943)
      END IF
*
  407 IF(LOOC.EQ.3) STOP
      INP=IOUT3
      IF(BIT.EQ.SINC) INP=INP1
*
*     A PICKUP FILE IS CREATED AFTER EACH EQUILIBRATION AND IS WRITTEN
*      IN THE FILE IOUT3 (E.G. EACH TIME WE PUT A NEW MINERAL IN, OR
*      KICK A MINERAL OUT). HOWEVER, WHEN WE FINALLY TRUELY EQUILIBRATE,
*      (NO MORE MINERALS ARE SUPER- OR UNDERSATURATED, AND WE ARE
*      READY FOR NEXT COOLING OR MIXING STEP), THE PICKUP FILE IS
*      WRITEN OVER THE LAST PICKUP FILE FOR THE LAST COMPLETED STEP.
*      THEREFORE, IOUT3 IS CREATED ONLY WHEN MINERALS ARE ADDED OR
*      OR KICKED OUT IN A GIVEN STEP.  THE FOLLOWING CREATES PICKUP
*      FILES.
      REWIND INP
      WRITE(INP,870) TITLE1
      WRITE(INP,870) TITLE2
      WRITE(INP,801)
      WRITE(INP,834) ERPC,FO2log,PFLUID,temp,SINC,SLIM,GASTOT
      WRITE(INP,802)
      WRITE(INP,836) tempC,TOTMIX,SOLMIN, RM,GASGRM,SUPRNT
      WRITE(INP,810)
      WRITE(INP,800) IMINT,IFRAC,IPUNCH,NLOOP,ISTEP,
     1 LIMSS,LOOC,ITREF,IDEAL,INCREM,INCP,MINSOL,OFLAG
      WRITE(INP,804)
      DO 520 I = 1,BSTOT
  520 WRITE(INP,805) SAQ(I),NAME(I),MTOT(I),MTRY(I),GAMMA(I),COMTOT(I)
      WRITE(INP,806)
      IF(NMT.EQ.0) GO TO 525
      WRITE(INP,823) (MIN(SEQ(I)), MINTRY(SEQ(I)), I=1,NMT)
  525 WRITE(INP,807)
      IF(NOOX.EQ.0) GO TO 535
      DO 530 I = 1, NOOX
  530 WRITE(INP,823)  NOMOX(I), WTPC(I)
  535 WRITE(INP,808)
      DO 540 I = 1,NSUP
  540 WRITE(INP,822) SUPNAM(I)
      WRITE(IOUT2,952)
*
      IF(BIT.EQ.0.) GO TO 190
*     TERMINAL OVERLOOK OF CALCULATIONS
*%    FO2=ACT(BSTP1)
      FO2LGP=Dlog10(ACT(BSTP1))
      WRITE(IOUT2,841) FO2LGP
  560 WRITE(IOUT2,840) temp,PFLUID
      IF(INCREM.NE.0) BIT = SINC
*
  598 CALL COOLER
*     ***********
*
*     STOPS ARE BUILT IN COOLER WHENEVER SLIM IS REACHED
*     MIXER, BULK CALLED FROM COOLER
*
      GO TO 190
*     STOP IN COOLER
*
  800 FORMAT (16I5)
  801 FORMAT(/'<  ERPC  >< logFO2 >< PFLUID ><  temp  ><  SINC  >'
     1 ,'<  SLIM  >< GASTOT >')
  802 FORMAT(/'<  tempC >< TOTMIX >< SOLMIN ><   RM   >< GASGRM >'
     1 ,'< SUPRNT >')
  803 FORMAT(I5,1X,8X,4X,2D16.9,F8.4,2X,D10.4)
  804 FORMAT(/,' saq> < name >    <     mtot     ><     mtry     >'
     1 ,'<gamma >  < comtot >')
  805 FORMAT(I5,1X,A8,4X,2D16.9,F8.4,2X,D10.4)
  806 FORMAT(/'< min  >  < mintry >')
  807 FORMAT(/'<nomox >  <  wtpc  >')
  808 FORMAT(/'<SUPNAM>')
  810 FORMAT(/' bsto ifra ipun nloo iste lims looc itre idea'
     1 ,' incr incp mins ofla')
  815 FORMAT (//)
  816 FORMAT (/,A80)
  817 FORMAT (8X, 6D12.5)
  818 FORMAT (//8X, 5D14.7)
  820 FORMAT (A8,7X,A1,I1, 7(F6.3,I2),F5.3,I2)
  822 FORMAT (A8, 5X,F9.4)
  823 FORMAT (A8, 2X,D10.4)
  830 FORMAT (8D10.4)
  832 FORMAT(A8,5F10.5)
  833 FORMAT(2F10.5,2D10.4,2F10.5)
  834 FORMAT(D10.4,3F10.5,2G10.4,2F10.5)
  836 FORMAT(F10.5,G10.4,D10.4,2F10.5,D10.4)
  835 FORMAT(5D16.9)
  837 FORMAT (I3,2X,A8,2X,7D10.5)
  838 FORMAT(I3,2X,A8,2X,F10.5)
  840 FORMAT(/1X,'CALCULATIONS DONE FOR ',F8.3,' DEG.C.;'
     1 ,1X,F8.3,' ATM')
  841 FORMAT(/5X,'THE log(FO2) IS: ',F7.2)
  854 FORMAT (A8,A1,1X,F7.2,I2,1X,7(F6.3,I2),F4.3,I2,F6.3,I2,F6.3)
  870 FORMAT (10A8)
  885 FORMAT(5X,A8,' IS SUPPRESSED IN GASTHERM DATA FILE')
  886 FORMAT(5X,A8,' IS SUPPRESSED IN THIS RUN')
  900 FORMAT (/1X, 'MINERAL   NO.    log K    TRIAL MOLES',
     1'    MINERAL   NO.    log K    TRIAL MOLES'/)
  905 FORMAT (2H  ,'NO. OF LOOPS = ',I3,5X,'MOLES OF GAS=',D12.7)
  906 FORMAT (2H  , 'MOLES OF GAS=',D12.7)
  910 FORMAT (2(1X, A8, 1X, I3, 3X, F8.3, 3X, D10.4, 4X))
  914 FORMAT (2H  , '                           MOLE      '
     1,'PARTIAL')
  915 FORMAT (2H  , 'SPECIES     MOLE NUMBER  FRACTION    '
     1,'PRESSURE    FUGACITY     log(F)    GAMMA'/)
  920 FORMAT (1X,I3,1X,A8,1X,4(D11.5,1X),F9.4,1X,F8.5)
  923 FORMAT (/)
  924 FORMAT (1X,I3, 1X, A8, 24X, 2(D11.5,1X,F10.4,1X), D11.5)
  927 FORMAT (2H  , 7X, 'TOTAL MOLES OF GAS = ',D12.7)
  930 FORMAT (2H  ,'MINERAL      MOLES    log MOLES    GRAMS    ',
     1'log GRAMS  WT PERCENT   VOL (CM3)',/)
  931 FORMAT (//2H  ,'MINERAL      MOLE FRACTION',6X,'ACTIVITY  '
     1 ,3X,'ACTIVITY COEF',2X,'log ACTIVITY',/)
  933 FORMAT (1H  , 'tempERATURE = ',F8.3,' C.; P(FLUID) =',F8.3,
     1' ATM; MIXER FRACTION = ',G10.4,/)
  935 FORMAT (2X  , A8 ,3X, 2(D10.4,2X,F7.3, 3X), 2(D10.4,3X))
  936 FORMAT(/2X,'TOTAL FRACTION OF MIXING SOLUTION ADDED: ',F12.6)
  937 FORMAT (2X  , A8 ,7X,2(D10.4,6X),F7.4,8X,F7.3)
  940 FORMAT (2H  ,'SPECIES    NO.    TOT.MOLES   MTRY    ',
     1'       GAMMA   MIXING SOLUTION',/)
  943 FORMAT (/)
  944 FORMAT (/)
  945 FORMAT (2H  , A8, I5, 4X, 2D12.5, F10.4,3X,D10.4)
  946 FORMAT('  NUMBER OF LOOPS USED = ',I4,'   LOOP LIMIT = ',I4)
  950 FORMAT (2H  ,'  SPECIES     NO.      log K',
     1 '             SPECIES     NO.      log K',/)
  952 FORMAT(10X,'PICKUP FILE WRITTEN ON DISK')
  955 FORMAT (4X, A8, 3X, I4, 3X, F10.4, 11X, A8, 3X, I4, 3X, F10.4)

  957 FORMAT (1X, ' COMPONENT     MOLE NUMBER     MOLE FRACTION',
     1'         PPM',/, 35X,'(TOTAL)'/)
  960 FORMAT (2H  , 10A8)
  961 FORMAT (1X, ' COMPONENT     TOTAL MOLES    GAS PHASE MOLES',
     1'   SOLID MOLES'/)
  962 FORMAT (1X, I2, 1X, A8, 2X, 3(D15.9, 2X))
  965 FORMAT(/16X,'WEIGHT OF GASEOUS SOLUTION: ',F12.6,' GRAMS'//)
  968 FORMAT(2X,'IF UNKNOWN, MINERAL DENSITIES ARE ASSUMED TO BE 2.5'/)
  970 FORMAT (2H  , 78('+'))
  971 FORMAT (2H  ,'AMOUNT OF SOLID REACTANT CONSUMED:', F12.6,1X,A8
     1  ,' log=',F9.4)
  972 FORMAT (1H  ,'AMOUNT OF SOLID PRODUCTS PRODUCED:', F12.6,
     1' GRAMS;',1X,'VOLUME:',F14.8,' CM3')
  973 FORMAT (2H  ,'AMOUNT OF ORIGINAL SOLUTION:      ',F12.6,' GRAMS'
     1 /,2X,'ORIGINAL SOLUTION / ROCK ADDED = ',F14.4,' (= GAS/ROCK)'/)
  974 FORMAT (/1H  ,  'INITIAL COMPOSITION OF SOLID PHASE (WT%):')
  976 FORMAT (/1H  ,  'INITIAL COMPOSITION OF SOLID PHASE (MOLE%):')
  975 FORMAT (/20(3(2X,A8,' = ',D10.4, ',')/))
  977 FORMAT (2H  , 7D11.4)
  999 FORMAT (1H1)
      END
      SUBROUTINE OXY
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DOUBLE PRECISION MTRY,MTOT,LGK,LGKM
      INTEGER BSTP1, BSTOT,OFLAG,SPEC
      COMMON/SLOUGH/ A(BG4,BG4),F(BG4),T(BG4,BG4),IV(BG4), N,ITREF
      COMMON/FROID/ SINC,BIT,SLIM,AKLOG,APFCL(BG2,5),
     1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4),
     2 XDH(4),ALK(5),LGK(BG2),temp,tempC, HEATC,COMTOT(BG3), HEATM,
     3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3),
     4 CVCCF(7,3), IFRAC
      COMMON /GEORGE/ AZERO(BG2),XI, NS, NST
      COMMON/MORT/MTRY(BG2),GAMMA(BG2),WATM,FO2,FO2log,TERMLN,KK
     1,OFLAG
      COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,Plog
     1,GASTOT,NMTT,LOOC
      COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8)
      COMMON/BULL/BSTOT,BSTP1,JGSTOT
      COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
      IF(ISTEP.EQ.1) WRITE(IOUT1,911)
  911 FORMAT (2H  , 40X, 'SUBROUTINE OXY WAS USED')
*
*
*
*     PSEUDO SUBROUTINE FH2CALC.
*      NO PF/GASTOT TERM IN "F" EQUATION BECAUSE OF FO2log CANCELS IT.
      F(1) = -.5*LGK(BSTP1)+Dlog10(GAMMA(2)*MTRY(2))
     1 -Dlog10(GAMMA(1)*MTRY(1))-0.5*FO2log
*     DERIVATIVES W/R H2 AND H2O
      A(1,1) = -1./MTRY(1)*TERMLN
      A(1,2) =  1./MTRY(2)*TERMLN
      DO 57 J = 3,N
   57 A(1,J) = 0.0
      DO 100 NS = BSTP1,NST
      CFF=0.
      IF(SPEC(NS,2).NE.1) GO TO 100
      CFF=COEF(NS,2)
      A(JGSTOT,1) = A(JGSTOT,1) + CFF*MTRY(NS)/MTRY(1)
  100 CONTINUE
      RETURN
      END
*
      SUBROUTINE MINERA
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DOUBLE PRECISION  MSL ,MTRY, MINTRY,MTOT,LGKM
      CHARACTER*8 MIN
      INTEGER SOLS, QES,OFLAG
      INTEGER SPM, SPECM, BSTOT,BSTP1
      INTEGER SPEC, SEQ
*     THIS SUBROUTINE SETS UP ALL PARTIAL DERIVATIVE MATRIX ELEMENTS FOR
*       MINERALS AND SOLID SOLUTIONS
*     EQUATIONS BSTP1 THROUGH BSTP1 + NMT ARE HYDROLYSIS Q EXPRESSIONS
*       IN log FORM E.G. F = -log K + log Q
      COMMON/MAXSIZ/ MAXSAQ,MAXNS,MAXNM,MAXSTO,MAXSSL,MAXN,MAXMEM,MAXIN
      COMMON/MORRIS/  MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW
     1 (BG1), tempK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1)
      COMMON/TRIX/ TMIN(BG3)
      COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,Plog
     1,GASTOT,NMTT,LOOC
      COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
      COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8)
      COMMON/BULL/BSTOT,BSTP1,JGSTOT
      COMMON/SLOUGH/ A(BG4,BG4),F(BG4), T(BG4,BG4),IV(BG4), N,ITREF
      COMMON/MORT/MTRY(BG2),GAMMA(BG2),WATM,FO2,FO2log,TERMLN,KK
     1,OFLAG
      COMMON/MORSEL/ GAMIN(BG1),B(BG1),SUMSOL(15),SOLTRY(15),
     1 DGDN(BG3,20),SOLS(15,20),ISOLS(BG1),LIMSOL,JSOLS(BG1),IST(15)
**    COMMON/FUGA/GASlog(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3,
**   1 NGAS,IMIX,IDEAL,ITGAS
      COMMON/PRECIS/ RMIN, RlogMN, RlogMX
*
      N = BSTOT + NMT + 1
      MINCT = MINCT + 1
      DO 76 I = 1, MAXSAQ
   76 TMIN(I) = 0.0
      DO 101 M = 1, NMT
      NM = SEQ(M)
      IF(NM.EQ.0) GO TO 105
      NMA = M + BSTOT + 1
      ITOT = ITM(NM)
      MSL = 0.
      DO 90 I = 1, ITOT
      SPM = SPECM(NM,I)
      CFM = COEFM(NM,I)
*     ADD MOLES OF MINERAL TO MASS BALANCE EQUATION
      F(SPM) = F(SPM) + CFM * MINTRY(NM)
      A(SPM,NMA) = CFM
*     COMPUTE Q FOR MINERAL CHARACTERISTIC REACTION
*     ACT CONTAINS FUGACITY OF GASES
      CHK=CFM * Dlog10(MTRY(SPM)*GAMMA(SPM)*PFLUID/GASTOT)
      MSL = MSL + CHK
*     FOR MASS ACTION EQS IN log FORM, DERIVATIVES ARE SAME FOR GASEOUS
*     AND AQUEOUS FORMS
      A(NMA,SPM) = CFM * TERMLN/MTRY(SPM)
*     DERIVATIVE OF MINERAL MASS ACTION EQUATION W/R TO GASTOT
      A(NMA,JGSTOT) = A(NMA,JGSTOT) - CFM*TERMLN/GASTOT
*     TMIN CONTAINS MOLES OF COMPONENTS IN SOLIDS
   92 TMIN(SPM) = TMIN(SPM) + CFM * MINTRY(NM)
   90 CONTINUE
      IF (MSL.LT.RlogMN.OR.MSL.GT.RlogMX) WRITE(IOUT1,201) MIN(NM),MSL
      IF(ISOLS(NM).EQ.0) GO TO 95
      JU = ISOLS(NM)
      JS = JSOLS(NM)
*     COMPUTE Q FOR MINERAL SOLID SOLUTIONS. NEGATIVE MINTRY NOT ALLOWED
*     SEE MAIN. IF B IS USED, GAMIN MUST BE 1, IF GAMIN IS NOT 1, B MUST
*     BE 1, THUS WE CAN PUT ALL THE EXPRESSION INTO ONE Dlog10.
      ACM =  Dlog10 (MINTRY(NM)/SOLTRY(JU) *GAMIN(NM))*B(NM)
      IF (ACM.LT.RlogMN.OR.ACM.GT.RlogMX)  WRITE(IOUT1,202) MIN(NM),ACM,
     1MINTRY(NM),SOLTRY(JU)
      MSL = MSL - ACM
   95 CONTINUE
*     F(NMA) IS MASS ACTION EQUATION FOR MINERAL (IN log FORM)
      F(NMA) = -LGKM(NM) + MSL
      A(NMA,NMA) = 0.0
      IF(ISOLS(NM).EQ.0) GO TO 101
*     BELOW FOR SOLID SOLUTIONS ONLY
*      NOW WE COMPUTE PARTIAL DERIVATIVES OF MASS ACTION EQ. W/ REPECT
*      TO MOLES OF MINERAL FOR SOLID SOLUTIONS.
      A(NMA,NMA) = -B(NM)*TERMLN*(1./MINTRY(NM) - 1./SOLTRY(JU))
      IS = IST(JS)
      DO 110 I = 1, IS
      LS = SOLS(JS,I)
*      HERE WE INCLUDE DERIVATIVE OF GAMIN(NM) W/ RESPECT TO MOLES OF
*      NM. THIS IS TO EASE CONVERGENCE IF NEEDED. DGDN ARE COMPUTED IN
*      THE ROUTINES THAT COMPUTE GAMIN. THESE CAN BE ZERO IN MOST CASES.
      IF(LS.NE.NM) GO TO 99
      A(NMA,NMA) = A(NMA,NMA) - TERMLN/GAMIN(NM)*DGDN(QES(NM),I)
      GO TO 110
*
   99 CONTINUE
      DO 111 KI = 1, NMT
*     WE SKIP IF END MEMBER NOT IN MATRIX
      IF(SEQ(KI).EQ.0) GO TO 112
      IF(LS.EQ.SEQ(KI)) NMD = KI + BSTOT + 1
      IF(LS.EQ.SEQ(KI)) A(NMA,NMD) = B(NM)*TERMLN/SOLTRY(JU)
     .                    - TERMLN/GAMIN(NM)*DGDN(QES(NM),I)
*      ABOVE LINE CONTAINS DERIVATIVE OF GAMIN(NM) W/ REPECT TO
*      MOLES OF ENDMEMBER DIFFERENT THAN NM.  SEE COMMENTS ABOVE.
  111 CONTINUE
  112 CONTINUE
  110 CONTINUE
  101 CONTINUE
  105 CONTINUE
*
  201 FORMAT (/1X,'log Q FOR ',A8,' IS ',F7.3)
  202 FORMAT (1X,'log A FOR ',A8,'IS',F7.3/1X,'MINTRY=',D12.5,'SOLTRY=
     1',D12.5)
      RETURN
      END
      SUBROUTINE MINSAT
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DOUBLE PRECISION MSL,MTRY,MINTRY, MTOT, LGKM ,LGK
      CHARACTER*8 MIN
      INTEGER SOLS, QES,OFLAG
      INTEGER SPM, SPECM, BSTOT,BSTP1,SPEC, SEQ
*     SUBROUTINE MINSAT DETERMINES WHETHER THERE IS SUPERSATURATION OF
*       THE PRESENT SOLUTION WITH ANY MINERALS OR SOLID SOLUTIONS.
*       IF MORE THAN ONE MINERAL IS SUPERSATURATED, IT ADDS THE ONE WITH
*       THE LOWEST AFFINITY TO THE MATRIX.  IT ALSO REMOVES ANY MINERALS
*       PREVIOUSLY IN MATRIX THAT ARE NOW UNDERSATURATED
      DIMENSION NMSAV(BG3),XLGQM(BG1), MTOSS(5),XM(15,20)
      COMMON/MAXSIZ/ MAXSAQ,MAXNS,MAXNM,MAXSTO,MAXSSL,MAXN,MAXMEM,MAXIN
      COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8)
      COMMON/BULL/BSTOT,BSTP1,JGSTOT
      COMMON/SLOUGH/ A(BG4,BG4),F(BG4), T(BG4,BG4),IV(BG4), N,ITREF
      COMMON/MORT/MTRY(BG2),GAMMA(BG2),WATM,FO2,FO2log,TERMLN,KK
     1,OFLAG
      COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,Plog
     1,GASTOT,NMTT,LOOC
      COMMON/FROID/ SINC,BIT,SLIM,AKlog,APFCL(BG2,5),
     1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4),
     2 XDH(4),ALK(5),LGK(BG2),temp,tempC, HEATC,COMTOT(BG3), HEATM,
     3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3),
     4 CVCCF(7,3), IFRAC
      COMMON/MORRIS/  MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW
     1 (BG1), tempK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1)
      COMMON/MORSEL/ GAMIN(BG1),B(BG1),SUMSOL(15),SOLTRY(15),
     1 DGDN(BG3,20),SOLS(15,20),ISOLS(BG1),LIMSOL,JSOLS(BG1),IST(15)
**    COMMON/FUGA/GASlog(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3,
**   1 NGAS,IMIX,IDEAL,ITGAS
      COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
      COMMON/PRECIS/ RMIN, RlogMN, RlogMX
      RT = 1.98719 * tempK
      AFOLD = 0.
      LREAP = 0
      MCNT = 0
      NEWCT = 0
      NMT = 0
      JACK=0
      DO 57 I = 1, MAXIN
   57 SEQ(I) = 0
      DO 201 NM = 1, NMTT
      MSL = 0.0
      ITOT = ITM(NM)
*     WE COMPUTE log Q
      DO 93 I = 1,ITOT
      CFM = COEFM(NM,I)
      SPM  = SPECM(NM,I)
      CSAVE=CFM*Dlog10(ACT(SPM))
      MSL = MSL + CSAVE
   93 CONTINUE
      JS = JSOLS(NM)
      IF(JS.EQ.0) GO TO 96
      GAMlog = Dlog10(GAMIN(NM))
      IF(ISOLS(NM).NE.0) GO TO 94
*
*     GASES OR MINERALS  NOT YET SUPERSATURATED
*     COMPUTES SUM OF MOLE FRACTIONS FOR SOLID SOLUTIONS THAT ARE NOT
*     YET SUPERSATURATED. IF SUMM ADDS TO 1 OR GREATER, WE WILL SUPERSAT
      QTERM =  (MSL - LGKM(NM))/B(NM) - GAMlog
      IF (QTERM.LE.RlogMN) QTERM=RlogMN
      IF (QTERM.GE.RlogMX) QTERM=RlogMX
      QTERM =10.**QTERM
      SUMSOL(JS) = SUMSOL(JS) + QTERM
      IS = IST(JS)
      JACK=1
      DO 73 J=1,IS
      IF(SOLS(JS,J).EQ.NM) XM(JS,J) = QTERM
   73 CONTINUE
      GO TO 96
*
*     COMPUTE Q FOR SATURATED SOLID SOLUTIONS
   94 JU = ISOLS(NM)
      ACM = Dlog10(DABS(MINTRY(NM)/SOLTRY(JU)))*B(NM) + GAMlog
      MSL = MSL - ACM
*
   96 XLGQM(NM) = MSL
      DIFM = XLGQM(NM) - LGKM(NM)
*     REMOVE UNDERSATURATED MINERALS FROM MATRIX
*     ================================================
*     IF SINGLE MINERAL WAS PREVIOUSLY IN MATRIX BUT GOT UNDERSATURATED
*     (NEGATIVE AMOUNT OF MOLES WERE COMPUTED) WE REMOVE IT
      IF(QES(NM).EQ.0) GO TO 97
      IF(MINTRY(NM).LT.0.) GO TO 195
*
*     IN CASE MASS BALANCE CONVERGES, BUT MASS ACTION FAILS (SHOULD NOT
*     HAPPEN WITH ERPC SUFFICIENTLY LOW), THE AFFINITY (-RTLN Q/K )
*     WILL BE TOO LARGE => MINERAL MIGHT NOT BELONG, OR PROBLEM. WE
*     KICK MINERAL OUT (USEFUL TO DUMP TOO SMALL SOL-SOL ENDMEMBERS)
      AF = -RT*DIFM*2.303
      IF(DABS(AF).LT.1.D-6) GO TO 101
      WRITE(IOUT1,800) MIN(NM), AF
  800 FORMAT(/5X,A8,' DOES NOT SATISFY ANYMORE THE MASS ACTION EQUATION'
     .        ,/13X,' WE REMOVE IT FROM MATRIX.  AFFINITY =', D10.4)

      GO TO 195
*
*     CHECK FOR SUPERSATURATION OF SINGLE MINERALS
*     ============================================
   97 CONTINUE
      DIFM=DIFM/SCALE(NM)
*     IF UNDERSATURATED, GO TO 201
      IF(DIFM.LT.1.D-10) GO TO 201
*     NOW, WE SELECT ONLY ONE MINERAL WITH SMALLEST AFFINITY IF
*     SEVERAL MINERALS ARE SUPERSATURATED
      IF(DIFM.LT.AFOLD) GO TO 99
      AFOLD = DIFM
      MSTORE = NM
   99 NEWCT = NEWCT + 1
      GO TO 201
*
*     101-195, RE-COUNT SATURATED MINERALS
  101 MCNT = MCNT + 1
      NMSAV(MCNT) = NM
      SEQ(MCNT) = NM
      QES(NM) = MCNT
      NMT = MCNT
      GO TO 201
*
*     HERE WE KICK MINERALS OUT OF THE MATRIX
  195 LREAP = LREAP + 1
      MASH = NM
      MTOSS(LREAP)= NM
      QES(NM) = 0
      MINTRY(NM) = 1.D-300
*
  201 CONTINUE
*     IF MINERALS WERE REMOVED FROM MATRIX
      IF(LREAP.NE.0) GO TO 107
*     IF NO NEW MINERALS SUPERSATURATE, AFOLD = 0
      IF(AFOLD.EQ.0.) GO TO 207
*     MORE THAN ONE NEW SUPERSAT MINERAL: TAKE THE MOST SUPERASTURATED
      IF(LOOC.EQ.0) GO TO 105
*     MORE THAN ONE NEW SUPERSAT MINERAL: BACK TITRATE
      IF(NEWCT.GE.2) GO TO 107
*
*     WE PUT IN THE MOST SUPERSAT. MINERAL (INDEX = MSTORE)
  105 MCNT = MCNT + 1
      NM = MSTORE
      SEQ(MCNT) = MSTORE
      QES(MSTORE) = MCNT
      NMSAV(MCNT) = MSTORE
      NMT = MCNT
      IF(JSOLS(NM).EQ.0) GO TO 107
*     IF WE JUST PUT IN A MINERAL WHICH BELONG TO A SOLID SOLUTION, THEN
*     WE NEED TO PUT IN THE REMAINING END MEMBERS.
      JS = JSOLS(NM)
      IS = IST(JS)
      MINTRY(NM) = 1. D-7
      DO 106 J = 1, IS
      M = SOLS(JS,J)
*     IF OTHER ENDMEMBERS WERE ALREADY IN MATRIX, WE SKIP
      IF(QES(M).NE.0) GO TO 106
      IF(M.NE.NM) MINTRY(M) = MINTRY(NM)*XM(JS,J)/SUMSOL(JS)
      IF(M.EQ.NM) GO TO 106
      MCNT = MCNT + 1
      SEQ(MCNT) = M
      QES(M) = MCNT
      NMSAV(MCNT) = M
      NMT = MCNT
  106 CONTINUE
      GO TO 107
*
*     NOW IN CASE A SOLID SOLUTION, WITH ALL END MEMBERS
*     INDEPENDENTLY UNDERSATURATED, SUPERSATURATES, WE PUT IT IN
  207 JSAV = 0
      IF(JACK.EQ.0) GO TO 107
      SOLQ = 1.0
      DO 213 NM = 1, NMTT
      JS = JSOLS(NM)
      IF(JS.EQ.0) GO TO 213
      IS = IST(JS)
      DO 210 J = 1, IS
      JK = SOLS(JS,J)
      IF(QES(JK).NE.0) GO TO 213
  210 CONTINUE
      IF(SUMSOL(JS).LT.SOLQ) GO TO 213
      SOLQ = SUMSOL(JS)
      JSAV = JS
  213 CONTINUE
      IF(JSAV.EQ.0) GO TO 107
      NEWCT = NEWCT + 1
      IS = IST(JSAV)
      DO 225 J = 1, IS
      M = SOLS(JSAV,J)
      MCNT = MCNT + 1
      SEQ(MCNT) = M
      QES(M) = MCNT
      NMSAV(MCNT) = M
      MINTRY(M) = 1.D-8*XM(JSAV,J)/SUMSOL(JSAV)
      NMT = MCNT
  225 CONTINUE
*
  107 CONTINUE
      LBOTH = LREAP + NEWCT
      IF(LBOTH.GE.2) GO TO 233
      IF(LREAP.NE.1) GO TO 311
      MINTRY(MASH) = 1.D-300
      GO TO 311
  233 IF(LOOC.EQ.0) GO TO 311
      IF(BIT.EQ.0.) BIT = SINC
      BIT =-(BIT/2.)
*     IF MORE THAN TWO MINERALS ARE UNDER OR SUPER SATURATED (LBOTH.GE.
*       2) THEN MUST REPLACE UNDERSAT ONES AND BACK UP
      IF(LREAP.EQ.0) GO TO 108
      DO 237 I = 1, LREAP
      NM = MTOSS(I)
      MCNT = MCNT + 1
      SEQ(MCNT) = NM
      QES(NM) = MCNT
      NMSAV(MCNT) = NM
  237 NMT = MCNT
      GO TO 108
*
  311 IF(LBOTH.NE.0) BIT = 0.
      IF(LBOTH.EQ.0) BIT = SINC
  108 CONTINUE
      CALL SQUIMP
      WRITE(IOUT1,923)
      WRITE(IOUT1,970)
      WRITE(IOUT1,923)
      IF(MCNT.EQ.0) GO TO 110
      WRITE(IOUT1,975)
      WRITE(IOUT1,980)
      DO 109 I = 1, MCNT
      NM = NMSAV(I)
      XLK = LGKM(NM)
      AF = (XLK - XLGQM(NM)) * RT * 2.303
      QKlog =  XLGQM(NM) - XLK
      QKSlog=QKlog/SCALE(NM)
      WRITE(IOUT1,977) NM,MIN(NM),XLK,XLGQM(NM),QKlog,QKSlog,AF
      XLGQM(NM) = 666.
  109 CONTINUE
      WRITE(IOUT1,923)
  110 WRITE(IOUT1,979)
      WRITE(IOUT1,980)
      DO 112 NM = 1, NMTT
      IF(XLGQM(NM).EQ.666.) GO TO 112
      XLK =LGKM(NM)
      AF = (XLK - XLGQM(NM)) * RT * 2.303
      QKlog = XLGQM(NM) - XLK
      QKSlog=QKlog/SCALE(NM)
**    IF(NM.GT.ITGAS) GO TO 130
**    GO TO 112
  130 IF(QKlog.GE.-5.)
     1   WRITE(IOUT1,977) NM,MIN(NM),XLK,XLGQM(NM),QKlog,QKSlog,AF
  112 CONTINUE
*
  923 FORMAT (/)
  970 FORMAT (1H  , 78('+'))
  975 FORMAT (2H  , 'THE FOLLOWING MINERALS ARE PRESENTLY IN MATRIX OR W
     1ERE JUST ADDED',//)
  977 FORMAT (2X,I3,2X,A8,2X,3(3X,F8.2),3X,F8.3,4X,D13.5,4X,F8.4)
  979 FORMAT (2H  , ' THE FOLLOWING MINERALS ARE PRESENTLY EXCLUDED FROM
     1 MATRIX',/3X,'MINERALS WITH log(Q/K) LESS THAN -5 ARE NOT LISTED'
     2 ,' BELOW'//)
  980 FORMAT    (5X,'   MINERAL        log K      log Q     log(Q/K)',
     1'  log(Q/K)/S    AFFINITY',/)
      RETURN
      END
      SUBROUTINE SQUOMP
*     TO SORT AND INITIALIZE SOLID SOLUTION ARRAYS.
*     NM  MINERAL INDEX: CURRENT SEQUENCE POSITION AS READ FROM
*         SOLTHERM FOR EACH PARTICULAR RUN.
*     JSOLS(NM)  SOLID SOLUTION INDEX: FIXED ID NUMBER (E.G.GASES = 6)
*     B(NM)  COEFFICIENT FOR MULTISITE MIXING
*     SOLS(JS,L) RETURNS (CONTAINS) NM VALUE, WITH JS = JSOLS(NM) AND
*         L = POSITION OF END ENDMEMBER FOR THAT JS SOLID SOLUTION
*         (1ST, 2ND ETC. END MEMBER)
*     IST(JS) TOTAL NUMBER OF END MEMBERS IN SOL SOL JS = JSOLS(NM)
*     ISOLS(NM) RETURNS (CONTAINS) THE POSITION, IN ARRAY OF CURRENTLY
*         SATURATED SOL SOL, OF THE SOLID SOLUTION CONTAINING NM
*     SOLTRY(J)  SUM OF MOLES OF ALL END MEMBERS IN SATURATED SOL SOL
*         J = ISOLS(NM)
*     SUMSOL(J)  SUM OF Q/K'S OF ALL END MEMBERS IN SATURATED SOL SOL
*         J = ISOLS(NM)
*     DGDN(J,I) DERIVATIVE OF GAMIN(NM) W/ RESPECT TO MOLES OF
*         ENDMEMBER I, WITH J = QES(NM) AND I = L ABOVE
*
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DOUBLE PRECISION LGKM,MINTRY
      CHARACTER*8 MIN
      INTEGER SPECM, SEQ, QES, SOLS
      COMMON/MAXSIZ/ MAXSAQ,MAXNS,MAXNM,MAXSTO,MAXSSL,MAXN,MAXMEM,MAXIN
      COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,Plog
     1,GASTOT,NMTT,LOOC
      COMMON/MORRIS/  MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW
     1 (BG1), tempK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1)
      COMMON/MORSEL/ GAMIN(BG1),B(BG1),SUMSOL(15),SOLTRY(15),
     1 DGDN(BG3,20),SOLS(15,20),ISOLS(BG1),LIMSOL,JSOLS(BG1),IST(15)
      COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
      WRITE(IOUT1,500)
      DO 180 J = 1,MAXIN
      DO 180 I = 1,MAXMEM
  180 DGDN(J,I) = 0.
      DO 183 J = 1,MAXSSL
      DO 182 I = 1,MAXMEM
  182 SOLS(J,I) = 0
      IST(J) = 0
  183 CONTINUE
      DO 186 NM = 1, NMTT
      JS = JSOLS(NM)
      IF(JS.EQ.0) GO TO 186
      L = 0
       IF(JS.LE.MAXSSL) GO TO 185
       WRITE(IOUT1,400) MIN(NM), JS, MAXSSL
  400  FORMAT(/5X,A8,' HAS A TOO LARGE SOLID SOLUTION INDEX (',I3,')',
     $ /5X,'MAXIMUM VALUE ALLOWED IS ',I3//)
       STOP
*
  185 L = L + 1
       IF(L.LE.MAXMEM) GO TO 187
       WRITE(IOUT1,410) MIN(NM), L, MAXMEM
  410  FORMAT(/5X,A8,' IS PART OF A SOLID SOLUTION WITH TOO MANY'
     $ ,' END MEMBERS (',I3,')'/5X,'MAXIMUM NUMBER ALLOWED IS ',I3//)
       STOP
*
  187 IF(SOLS(JS,L).NE.0) GO TO 185
      SOLS(JS,L) = NM
      IST(JS) = L
      WRITE(IOUT1,510) MIN(NM), JS, B(NM)
  186 CONTINUE
  500 FORMAT(/6X,'THE FOLLOWING MINERALS AND/OR GASES WITH IDENTICAL',
     1/6X,'INDEX ARE IN SOLID SOLUTION:'//,8X,
     2 'MINERAL    INDEX    B(NM)')
  510 FORMAT(8X,A8,3X,I5,3X,F6.3)
*
      ENTRY SQUIMP
      DO 163 NM = 1, NMTT
  163 ISOLS(NM) = 0
      L = 0
      DO 189 I = 1,MAXSSL
*     RESET SUMSOL TO ZERO FOR MINSAT
      SUMSOL(I) = 0.
      M1 = SOLS(I,1)
      IF(IST(I).GE.2) GO TO 188
      IF(IST(I).EQ.0) GO TO 189
      SOLS(I,1) = 0
      JSOLS(M1) = 0
      IST(I) = 0
      GO TO 189
  188 CONTINUE
      IS = IST(I)
      L = L + 1
      SOLTRY(L) = 0.0
      DO 196 J = 1, IS
      LS = SOLS(I,J)
*     SKIP IF ONE END MEMBER IS NOT IN MATRIX
      IF(QES(LS).EQ.0) GO TO 196
      SOLTRY(L) = SOLTRY(L) + MINTRY(LS)
      ISOLS(LS) = L
  196 CONTINUE
  189 CONTINUE
      LIMSOL = L
      RETURN
      END
      SUBROUTINE COOLER
*     COOLER IS USED TO GO ON TO THE NEXT STEP, WHICH CAN BE:
*     T CHANGE, P CHANGE, P-T CHANGE, OR COMPOSITION CHANGE
*     NOTE, CURRENTLY P CHANGE ONLY AFFECTS GASES
*
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DOUBLE PRECISION MINTRY,MTRY,MTOT,LGKM,LGK
      CHARACTER*8 MIN,NOMOX
      INTEGER QES,SEQ,QAS,SAQ,SPECM,BSTOT,SPEC,SAQTOT,SPECOX,BSTP1
     1,OFLAG
      COMMON /GEORGE/AZERO(BG2),XI, NS, NST
      COMMON/TRIX/ TMIN(BG3)
      COMMON/MORT/MTRY(BG2),GAMMA(BG2),WATM,FO2,FO2log,TERMLN,KK
     1,OFLAG
      COMMON/MORRIS/  MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW
     1 (BG1), tempK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1)
      COMMON/TURNIP/ SAQ(BG3), QAS(BG3), SAQTOT
      COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8)
      COMMON/BULL/BSTOT,BSTP1,JGSTOT
      COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,Plog
     1,GASTOT,NMTT,LOOC
      COMMON/FROID/ SINC,BIT,SLIM,AKlog,APFCL(BG2,5),
     1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4),
     2 XDH(4),ALK(5),LGK(BG2),temp,tempC, HEATC,COMTOT(BG3), HEATM,
     3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3),
     4 CVCCF(7,3), IFRAC
**    COMMON/FUGA/GASlog(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3,
**   1 NGAS,IMIX,IDEAL,ITGAS
      COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
      COMMON/PP/TS(200),PS(200),NPP,INCP
      COMMON/BESS/ WTPC(BG3), COEFOX(BG3,8), XMW(BG3),
     1 NOMOX(BG3),ITOX(BG3),NOOX, SPECOX(BG3,8),MINSOL
*
*
*     DECREASE P AT CONSTANT T
   65 IF(INCP.EQ.0) GO TO 80
      PFLUID = PFLUID + BIT
      WRITE(IOUT1,904) BIT , PFLUID
      WRITE(IOUT2,846) BIT
      Plog = Dlog10(PFLUID)
      IF(SINC.GT.0..AND.PFLUID.LE.SLIM) GO TO 160
      IF(SINC.LT.0..AND.PFLUID.GE.SLIM) GO TO 160
      WRITE(IOUT1,628) SLIM, PFLUID
      WRITE(IOUT2,628) SLIM, PFLUID
      STOP
*
*     TITRATE SOLIDS AT CONSTANT T
   80 IF(MINSOL.EQ.0) GO TO 85
      CALL BULK
      IF(SINC.GT.0..AND.TOTMIX.LE.SLIM) GO TO 160
      IF(SINC.LT.0..AND.TOTMIX.GE.SLIM) GO TO 160
      WRITE(IOUT1,627) SLIM, TOTMIX
      WRITE(IOUT2,627) SLIM, TOTMIX
      STOP
*
*     MIX THE SOLUTION, WITH OR WITHOUT T CHANGE
   85 IF(tempC.EQ.0.) GO TO 90
      CALL MIXER
      IF(SINC.GT.0..AND.TOTMIX.GT.SLIM) GO TO 88
      IF(SINC.LT.0..AND.TOTMIX.LT.SLIM) GO TO 88
      IF(tempC.EQ.temp) GO TO 160
      GO TO 100
   88 WRITE(IOUT1,627) SLIM, TOTMIX
      WRITE(IOUT2,627) SLIM, TOTMIX
      STOP
*
*     CHANGE tempERATURE
   90 temp = temp + BIT
      WRITE(IOUT1,903) BIT , temp
      WRITE(IOUT2,842) BIT
      IF(temp.LE.1227.) GO TO 92
      WRITE(IOUT1,625) temp
      WRITE(IOUT2,625) temp
      STOP
   92 IF (OFLAG.NE.0) GO TO 91
      FO2 = 0.
      FO2log = 0.
   91 IF(SINC.GT.0..AND.temp.LE.SLIM) GO TO 100
      IF(SINC.LT.0..AND.temp.GE.SLIM) GO TO 100
      WRITE(IOUT1,626) SLIM,temp
      WRITE(IOUT2,626) SLIM, temp
      STOP
*
*     BELOW WE RECALCULATE ALL THERMO DATA IF T OR P WAS CHANGED
  100 CONTINUE
      tempK=temp + 273.15
**    DO 50 I=1,3
** 50  CVB(I) = CVBCF(I,1)/tempK/tempK + CVBCF(I,2)/tempK + CVBCF(I,3)
**    DO 52 I=1,7
** 52  CVC(I) = CVCCF(I,1)/tempK/tempK + CVCCF(I,2)/tempK + CVCCF(I,3)
**    DO 54 I=1,NGAS
**     VB(I)= VBCOE(I,1)/tempK/tempK + VBCOE(I,2)/tempK + VBCOE(I,3)
** 54  VC(I)= VCCOE(I,1)/tempK/tempK + VCCOE(I,2)/tempK + VCCOE(I,3)
**
**    DO 127 J = 1, 9
**    DO 123 I = 1,5
**123 COE(I) = XPFC(J,I)
**    CALL POWER(COE,temp,AKlog)
**    XF(J) = AKlog
  127 CONTINUE
*     RECOMPUTE DERIVED SPECIES DATA WITH NEW T
      NBSP1 = SAQTOT + 1
      DO 107 NS = NBSP1, NST
      DO 109 I = 1, 5
  109 COE(I) = APFCL(NS,I)
  110 CALL POWERK(COE,tempK,AKlog)
      LGK(NS) = AKlog
  107 CONTINUE
*
*     THERE ARE TWO WAYS CHANGE P AS F(T): POWER FUNCTION P AS
*     F(T), WHOSE COEFFICIENTS ARE READ IN GASTHERM; A SET OF GIVEN
*     P-T POINTS READ FROM AN ACCESSORY FILE (INP3).
*     IF NEITHER OF THESE OPTIONS ARE USED, P IS CONSTANT.
      IF(XF(9).NE.0.) PFLUID=XF(9)
      IF(NPP.NE.0) CALL PRESS(PFLUID,temp)
      Plog = Dlog10(PFLUID)
*     HERE WE RECOMPUTE GAS AND MINERAL DATA WITH NEW P AND T
      DO 117 NM = 1, NMTT
      DO 111 I = 1, 5
  111 COE(I) = BPFCL(NM,I)
  114 CALL POWERK(COE,tempK,AKlog)
      LGKM(NM) = AKlog
  117 CONTINUE
*
*      BIT NOT EQUAL TO SINC ONLY IF BACK TITRATION
  160 IF(BIT.NE.SINC) GO TO 210
*     FRACTIONATION OPTIONS BELOW
*      NO FRACTIONATON
      IF(IFRAC.GT.0) GO TO 202
      WRITE(IOUT1,907)
      GO TO 210
*
*      FRACTIONATES SOLIDS ONLY
  202 DO 170 NS = 1, SAQTOT
  170 MTOT(NS) = MTOT(NS) - TMIN(NS)
      WRITE(IOUT1,908)
      WRITE(IOUT2,908)
      GO TO 210
*
  210 CONTINUE
*     NEXT TWO STATEMENTS CAUSE MINSAT TO USE THE LOWEST K/Q FOR
*     MINERAL SELECTION RATHER THAN T-CHANGE TO FIND SINGLE
*     SUPERSATURATED MINERAL.
      IF(DABS(BIT).LT..05.AND.(BIT*SINC).LT.0.) LOOC=0
      IF (DABS(BIT/SINC).LT..24)  LOOC=0
      RETURN
  625 FORMAT(///5X,'Temperature greater than 1227 C. temp =',F8.3,
     *         /5X,'EXECUTION STOPPED.'//)
  626 FORMAT(///5X,'Temperature reached given limit of ',F8.3,
     *         /5X,'Temp = ',F8.3,'  EXECUTION STOPPED.'//)
  627 FORMAT(///5X,'MIXING FRACTION REACHED GIVEN LIMIT OF ',G10.4,
     *         /5X,'TOTMIX = ',G10.4,'  EXECUTION STOPPED.'//)
  628 FORMAT(///5X,'PRESSURE REACHED GIVEN LIMIT OF ',F8.3,
     *         /5X,'PFLUID = ',F8.3,'  EXECUTION STOPPED.'//)
  630 FORMAT(///5X,'ENTHALPY REACHED GIVEN LIMIT OF ',F10.4,
     *         /5X,'ENTH  =',F10.4,'  EXECUTION STOPPED.'//)
*
  842 FORMAT(5X,' Next temperature change of ',F8.3,'deg.C.',///)
  846 FORMAT(5X,' Next pressure change of ',F8.3,' atm',///)
  847 FORMAT(5X,' Next enthalpy change of ',F10.4,'KCAL',/,'  TOTAL ',
     1' ENTHALPY ADDED (INCL. STARTING ENTHALPY)',F10.4,' KCAL'///)
  903 FORMAT ('  Temperature changed by ',F8.3,' DEGREES TO ',F8.3,' DE
     1GREES C.',/)
  904 FORMAT ('  PRESSURE CHANGED BY ',F8.3,' ATM TO ',F8.3,' ATM'/)
  905 FORMAT ('  ENTHALPY CHANGED BY ',F10.4,' KCAL TO ',F10.4,' KCAL'/,
     1'  TOTAL ENTHALPY ADDED (INCL. STARTING ENTHALPY)',F10.4,' KCAL'/)
  907 FORMAT (5X,'SOLIDS AND GASES CARRIED')
  908 FORMAT(5X,'SOLIDS ONLY LEFT BEHIND')
  909 FORMAT(5X,'GASES ONLY LEFT BEHIND')
  910 FORMAT(5X,'SOLIDS AND GASES LEFT BEHIND')
      END
      SUBROUTINE MIXER
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DOUBLE PRECISION MTOT,LGK,LGKM
      INTEGER BSTOT,SPEC,BSTP1
      COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8)
      COMMON/BULL/BSTOT,BSTP1,JGSTOT
      COMMON/FROID/ SINC,BIT,SLIM,AKlog,APFCL(BG2,5),
     1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4),
     2 XDH(4),ALK(5),LGK(BG2),temp,tempC, HEATC,COMTOT(BG3), HEATM,
     3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3),
     4 CVCCF(7,3), IFRAC
      COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,Plog
     1,GASTOT,NMTT,LOOC
      COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
*
*     TITRATE MIXER (COOL) SOLUTION INTO ORIGINAL (HOT) SOLUTION
      WRITE(IOUT1,925) TOTMIX
      WRITE(IOUT2,925) TOTMIX
      DO 101 I = 1, BSTOT
  101 MTOT(I) = MTOT(I) + BIT * COMTOT(I)
      IF (temp.EQ.tempC) GO TO 102
      temp = (temp*(1+TOTMIX)+tempC*BIT)/(1+TOTMIX+BIT)
*     Temperature is currently a simple weighted sum of the mixed
*     comonents.  NO ACCOUNT FOR NON-LINEARITIES IN HEAT CAPACITIES
*     IS TAKEN.
  102 TOTMIX=TOTMIX+BIT
*
      WRITE(IOUT2,920)  TOTMIX,temp
      WRITE(IOUT1,920)  TOTMIX,temp
      RETURN
*
  920 FORMAT(5X,'NEW MIXER FRACTION: ',G10.4,/,
     1 5X,'RESULTING tempERATURE: ',F8.3,' DEG.C.'/)
  925 FORMAT(/5X,'MIXER FRACTION: ',G10.4,/)
      END
      SUBROUTINE BULK
*     BULK COMPUTES COMPOSITION OF SOLID PHASE IN TERMS OF COMPONENT
*     SPECIES, GIVEN WT. PERCENTS OF OXIDES, ELEMENTS, OR MINERALS.
*     IF MINERAL REACTANTS ARE USED, THE MINERAL NAMES SHOULD BE WRITTEN
*     TO MATCH THE FORM USED IN GASTHERM IF YOU WANT MINSOL = 2
*     OPTION TO WORK
*
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DOUBLE PRECISION  MTOT,MINTRY,LGK
      INTEGER SPECM, SEQ, QES, SPEC, SPECOX, BSTOT,BSTP1
      CHARACTER*8 AMOLE,MIN,NOMOX,TYP
      DIMENSION  SOLMOL(BG3)
      COMMON/MORRIS/  MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW
     1 (BG1), tempK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1)
      COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8)
      COMMON/BULL/BSTOT,BSTP1,JGSTOT
      COMMON/FROID/ SINC,BIT,SLIM,AKlog,APFCL(BG2,5),
     1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4),
     2 XDH(4),ALK(5),LGK(BG2),temp,tempC, HEATC,COMTOT(BG3), HEATM,
     3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3),
     4 CVCCF(7,3), IFRAC
      COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
      COMMON/BESS/ WTPC(BG3), COEFOX(BG3,8), XMW(BG3),
     1 NOMOX(BG3),ITOX(BG3),NOOX, SPECOX(BG3,8),MINSOL
      COMMON /GEORGE/ AZERO(BG2),XI, NS, NST
*
      DATA TYP/'GRAMS   '/, AMOLE/'MOLES   '/
*
      OLDTOT = TOTMIX
      TOTMIX  = TOTMIX  + BIT
      DO 102 I = 1, BSTOT
  102 SOLMOL(I) = 0.0
      DO 105 I = 1, NOOX
      IF(XMW(I).EQ.0.) WRITE(IOUT1,600) NOMOX(I)
  600 FORMAT(//5X,'WARNING, MOLECULAR WEIGHT OF REACTANT ',A8,' IS 0.'/
     1  5X,'CHECK MINOX DATA FILE FOR ERRORS OR DOUBLE ENTRIES.'//)
      OXMOL =  BIT * WTPC(I) * .01/XMW(I)
      IF(MINSOL.GT.2) OXMOL = BIT*WTPC(I)* 0.01
*     MINSOL = 1  TITRATION IN GRAMS (SINC IN GRAMS, WTPC IN WT%)
*     MINSOL = 2  ACTIVATES NON-TITRATION OF EQUILIBRATED REACTANTS
*                 WITH GRAM TITRATION OPTION
*     MINSOL = 3 TITRATION OF IN MOLES  (SINC IN MOLES, WTPC IN MOLE%)
*     MINSOL = 4 ACTIVATES NON-TITRATION OF EQUIL. REACTANTS W/ MOLE
*                 TITRATION OPTION
      IF(MINSOL.NE.2.AND.MINSOL.NE.4) GO TO 95
      DO 90 J = 1, NMT
      NM = SEQ(J)
      IF(NOMOX(I).NE.MIN(NM)) GO TO 90
      OXMOL = 0.0
      WRITE(IOUT1,900) MIN(NM)
      WRITE(IOUT2,900) MIN(NM)
      GO TO 95
   90 CONTINUE
  900 FORMAT(2H  ,'REACTANT  ',A8,'  NOT TITRATED BECAUSE EQUILIBRATED')
*
   95 CONTINUE
      ITOT = ITOX(I)
      DO 103 IX = 1, ITOT
      CX = COEFOX(I,IX)
      JSPX = SPECOX(I,IX)
      IF(JSPX.EQ.0) GO TO 105
      SOLMOL(JSPX) = SOLMOL(JSPX) + CX * OXMOL
  103 CONTINUE
  105 CONTINUE
*   WE TITRATE SOLID INTO SOLUTION AND RECALCULATE MASS BALANCE
      DO 205 I = 1,BSTOT
      MTOT(I) = MTOT(I) + SOLMOL(I)
  205 CONTINUE
*
      IF(MINSOL.GT.2) TYP = AMOLE
      WRITE(IOUT1,910) BIT, TYP
      WRITE(IOUT2,910) BIT, TYP
      WRITE(IOUT2,920)  OLDTOT,TYP,TOTMIX,TYP
      WRITE(IOUT1,920)  OLDTOT,TYP,TOTMIX,TYP
      RETURN
*
  910 FORMAT(/5X,'NEXT SOLID TITRATION INCREMENT OF: ',F10.6,1X,A8)
  920 FORMAT(5X,'CURRENT TOTAL AMOUNT OF SOLID ADDED ',F11.6,1X,A8
     1 /5X,'NEW TOTAL AMOUNT OF SOLID ADDED ',F11.6,A8//)
      END
*
      SUBROUTINE POWER(COE,temp, AKlog)
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DIMENSION COE(5)
      COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
      COMMON/PRECIS/ RMIN, RlogMN, RlogMX
      AKlog = COE(1)
      DO 107 I = 1, 4
      J = I + 1
      AKlog = AKlog + COE(J) * temp**I
  107 CONTINUE
      RETURN
      END
*
      SUBROUTINE POWERK(COE,tempK,AKlog)
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DIMENSION COE(5)
      COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
      COMMON/PRECIS/ RMIN, RlogMN, RlogMX
      AKlog = COE(1)
      AKlog = AKlog + COE(2)/tempK
      AKlog = AKlog + COE(3)*tempK
      AKlog = AKlog + COE(4)/tempK**2
      AKlog = AKlog + COE(5)*Dlog10(tempK)
      RETURN
      END
      SUBROUTINE SLUD
      IMPLICIT DOUBLE PRECISION (A-G,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      COMMON/SLOUGH/ A(BG4,BG4),F(BG4),T(BG4,BG4),IV(BG4),N,ITREF
*
      DO 8100 J=1,N
      DO 8100 I=1,N
      T(I,J) = A(I,J)
 8100 CONTINUE
 8110 CONTINUE
*
      IV(N) = 1
      DO 8260 K=1,N
      PIV = DABS(T(K,K))
      IF(K.GE.N) GO TO 8260
*
      L = K
      KP1 = K+1
      DO 8200 I=KP1,N
      IF(PIV.GE.DABS(T(I,K))) GO TO 8200
      PIV = DABS(T(I,K))
      L = I
 8200 CONTINUE
      IV(K) = L
*
      TMP = T(L,K)
      IF(K.NE.L) GO TO 8210
      IF(PIV.GT.0.0) GO TO 8220
      IV(N) = 0
      GO TO 8260
 8210 CONTINUE
      IV(N) = -IV(N)
      T(L,K) = T(K,K)
      T(K,K) = TMP
 8220 CONTINUE
*
      TMP = -TMP
      DO 8230 I=KP1,N
      T(I,K) = T(I,K)/TMP
 8230 CONTINUE
*
      DO 8250 J=KP1,N
      TMP = T(L,J)
      T(L,J) = T(K,J)
      T(K,J) = TMP
      DO 8240 I=KP1,N
      T(I,J) = T(I,J) + T(I,K) * TMP
 8240 CONTINUE
 8250 CONTINUE
 8260 CONTINUE
      IF(PIV.EQ.0.0) IV(N) = 0
      RETURN
      END
      SUBROUTINE SIR(X)
      IMPLICIT DOUBLE PRECISION (A-G,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DIMENSION X(BG4), B(BG4), R(BG4)
      COMMON/SLOUGH/ A(BG4,BG4),F(BG4),T(BG4,BG4),IV(BG4),N,ITREF
      COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
*
*     B IS THE VECTOR CONTAINING THE RIGHT-HAND SIDE OF THE
*     SYSTEM OF LINEAR EQUATIONS WHICH IS  -F. NOTE THE MINUS SIGN WHICH
*     IS COMPATIBLE WITH XNEW = XOLD + X IN MAIN (NOTE PLUS SIGN)
*
      DO 150 I=1,N
      B(I)=-F(I)
  150 CONTINUE
*
      IER = 0
      DO 8100 I=1,N
      X(I) = B(I)
 8100 CONTINUE
      IF(IV(N).NE.0) GO TO 8500
*
      WRITE(IOUT1,915)
  915 FORMAT(/,5X,'******** WARNING ********',/,2X,
     1'COMPUTATIONAL SINGULARITY, IV(N) = 0.0',/)
      STOP
*
 8500 CONTINUE
      QUOT = 1.0
*
      CALL SBS(X)
      IF(ITREF.EQ.0) RETURN
 8200 CALL SAXMB(X,B,R)
      IF(ITREF.EQ.1) RETURN
      CALL SBS(R)
      IER = IER+1
      IF(ISTEP.EQ.1) WRITE(IOUT1,909)  IER
  909 FORMAT(5X,' RESIDUAL VECTOR (IER =',I2,')')
      IF(ISTEP.EQ.1) WRITE(IOUT1,973) (R(I), I=1,N)
*
      XN1 = 0.0
      RN1 = 0.0
      DO 8210 I=1,N
      XN1 = XN1 + DABS(X(I))
      RN1 = RN1 + DABS(R(I))
      X(I) = X(I) - R(I)
 8210 CONTINUE
*
      IF(XN1.GT.0.0) RN1 = RN1/XN1
      IF((RN1 + RN1).LE.QUOT) GO TO 8220
      IER = -IER
      WRITE(IOUT1,908)
  908 FORMAT(2X,'*************** WARNING *******************',/
     1,5X,'THE SYSTEM IS ILL-CONDITIONNED')
      RETURN
 8220 CONTINUE
      QUOT = RN1
      IF(1.0 + SNGL(RN1).NE.1.0) GO TO 8200
      RETURN
      ENTRY DUMP
      WRITE(IOUT1,923)
      WRITE(IOUT1,974) N
      DO 507 I = 1, N
      WRITE(IOUT1,973) (A(I,J), J = 1, N), F(I)
  507 CONTINUE
      RETURN
      ENTRY DUMP2
      WRITE(IOUT1,910) IER
      WRITE(IOUT1,922)
      RETURN
  910 FORMAT(5X,'NUMBER OF ITERATIVE REFINEMENTS : ',I4)
  922 FORMAT(//)
  923 FORMAT (///)
  973 FORMAT (2H  , 7D11.4)
  974 FORMAT (2H  , 'MATRIX OF A(I,I), F(I)  I=',I2//)
      END
      SUBROUTINE SBS(B)
      IMPLICIT DOUBLE PRECISION (A-G,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DIMENSION B(BG4)
      COMMON/SLOUGH/ A(BG4,BG4),F(BG4),T(BG4,BG4),IV(BG4),N,ITREF
*
      DO 8110 K=1,N
      IF(K.GE.N) GO TO 8110
      L = IV(K)
      TMP = B(L)
      B(L) = B(K)
      B(K) = TMP
      KP1 = K + 1
      DO 8100 I=KP1,N
      B(I) = B(I) + T(I,K)*TMP
 8100 CONTINUE
 8110 CONTINUE
*
      K = N
 8200 B(K) = B(K)/T(K,K)
      IF(K.LE.1) RETURN
      TMP = -B(K)
      KP1 = K
      K = K-1
      DO 8210 I=1,K
      B(I) = B(I) + T(I,KP1) *TMP
 8210 CONTINUE
      GO TO 8200
      END
      SUBROUTINE SAXMB(X,B,R)
      IMPLICIT DOUBLE PRECISION (A-G,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
      DIMENSION X(BG4), B(BG4), R(BG4)
      COMMON/SLOUGH/ A(BG4,BG4),F(BG4),T(BG4,BG4),IV(BG4),N,ITREF
*
      DO 8110 I=1,N
      SUM = -B(I)
      DO 8100 J = 1,N
      SUM = SUM + A(I,J)*X(J)
 8100 CONTINUE
      R(I) = SUM
 8110 CONTINUE
      RETURN
      END
      SUBROUTINE PRESS(P,T)
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      INTEGER BG1,BG2,BG3,BG4
      PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
*
*     THIS SUBROUTINE USES THE POINTS OF A GIVEN P-T CURVE, READ
*     IN INP3, AND INTERPOLATES LINEARLY BETWEEN THE TWO CLOSEST
*     GIVEN POINTS TO RETURN THE PRESSURE. N.S. ADDED OCT 83
*
      COMMON/PP/TS(200),PS(200),NPP,INCP
      COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH
*
*
      IF(T.GT.TS(NPP)) RETURN
*     GET THE APPROXIMATE POSITION OF POINT IN CURVE ARRAY
      ITEM=IDINT((T-TS(1))*NPP/(TS(NPP)-TS(1)))+1
      I=ITEM
      IF(I.GT.NPP) I=NPP
      IF(TS(I)-T) 40,30,50
   30 P=PS(I)
      RETURN
   40 K=I+1
      IF(TS(K).GE.T) GO TO 60
      I=K
      GO TO 40
   50 K=I-1
      IF(TS(K).LE.T) GO TO 60
      I=K
      GO TO 50
   60 CONTINUE
      P=((T-TS(I))*(PS(K)-PS(I))/(TS(K)-TS(I))) + PS(I)
      RETURN
      END
*     SUBROUTINE PHI(PF)
*     *****THIS SUBROUTINE IS NOT IN USE AT PRESENT..  IN ORDER TO USE
*     YOU MUST FIRST ADJUST ALL COMMON BLOCKS AND ARRAY SIZES TO CONFORM
*     WITH PRESENT GASWORKS VERSION
*
**    IMPLICIT DOUBLE PRECISION (A-H,O-Z)
**    INTEGER BG1,BG2,BG3,BG4
**    PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
**    CHARACTER*8 MIN
**    COMMON/MORRIS/  MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW
**   1 (BG1), tempK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1)
**    COMMON/MORT/MTRY(310),GAMMA(310),WATM,FO2,FO2log,KK
**   1,OFLAG
**    COMMON/FUGA/GASlog(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3,
**   1 NGAS,IMIX,IDEAL,ITGAS
**    COMMON/MORSEL/ GAMIN(BG1),B(BG1),SUMSOL(15),SOLTRY(15),
**   1 DGDN(BG3,20),SOLS(15,20),ISOLS(BG1),LIMSOL,JSOLS(BG1),IST(15)
**    COMMON/VAL/BVAL(20), CVAL(20), BMIX, CMIX
**    INTEGER SOLS, SPECM, SEQ, QES, OFLAG
**    DOUBLE PRECISION MINTRY
**    FUNC1(A1,A2,A3) = W1*A1 + W2*A2 + W3*A3
**    FUNC2(A1,A2,A3,A4,A5,A6) = W1*W1*A1 + W2*W2*A2 + W3*W3*A3
**   &    + 2.*( W1*W2*A4 + W1*W3*A5 + W2*W3*A6 )
**
**    THIS SUBROUTINE COMPUTES FUGACITY COEFFICIENTS FOR GASES.
**    WE ASSUME NON IDEAL MIXING FOR MIXTURES OF H2O, CO2 AND CH4
**    AND IDEAL MIXING FOR THE OTHER GASES.
**    SINCE H2O, CO2 AND CH4 GENERALLY ACCOUNT FO 99 PERCENT OF THE
**    THE GAS PHASE, THE RESULTING ERROR WHEN COMPUTING GAS COMPOSITIONS
**    IS ACCEPTABLE. THE FUGACITY COEFFICIENTS, PHI, FOR H2O, CO2 AND CH
**    ARE THEREFORE DEPENDENT ON T, P AND MOLE FRACTION Y OF THESE
**    THREE COMPONENTS. PHI FOR OTHER GASES ARE ONLY DEPENDENT ON T AND
**    P (PHI OF PURE GAS = PHI OF GAS IN MIXTURE)
**    ENTRIES FOR H2O, CO2 AND CH4 GAS HAVE TO BE AT THE TOP OF THE GAS
**    LIST IN DATA FILE SOLTHERM, H2O FIRST, CO2 SECOND AND CH4 THIRD.
**    FUGACITY COEFFICIENTS FOR H2O, CO2 AND CH4 ARE CALCULATED WITH
**    THEIR MOLE FRACTIONS NORMALIZED TO 1 (W1 W2 W3).
**    IMIX IS A FLAG WHICH INDICATES WHICH OF THESE THREE GASES ARE
**    PRESENT (H2O IS ALWAYS PRESENT) :
**       IMIX=3  H2O-CO2-CH4           NON IDEAL MIXING
**       IMIX=2  H2O-CO2  (NO CH4)     NON IDEAL MIXING
**       IMIX=1  H2O ONLY (NO CH4,CO2) IDEAL MIXING WITH OTHER GASES
**
**    PHI ARE CALCULATED FROM COMPRESSIBILITY DATA, USING A VIRIAL
**    EQUATION IN PRESSURE   Z = P*V/R/T = 1 + VB*PF + VC*PF**2.
**    FOR NON IDEAL MIXING WE USE Z = ZMIX. FOR A MULTICOMPONENT
**    MIXTURE ZMIX = 1 + BMIX*P + CMIX*P*P WHERE BMIX AND CMIX ARE
**    EXPRESSED IN TERMS OF CROSS VIRIAL COEFFICIENTS BIJ' AND CIJK'.
**    VALUES FOR BIJ'S AND CIJK'S WHERE OBTAINED FROM LEAST SQUARE
**    FIT OF PUBLISHED PVT OR SOLUBILITY DATA. PHI(I)'S ARE THEN
**    OBTAINED BY TAKING THE APPROPRIATE PARTIAL DERIVATIVES.
**    SEE SPYCHER AND REED, IN PREPARATION.
**    CROSS VIRIAL COEFFICIENTS ARE ASSIGNED AS FOLLOW (1=H2O, 2=CO2,
**    3=CH4):
**     CVB(1) : B12
**     CVB(2) : B13
**     CVB(3) : B23
**     CVC(1) : C112
**     CVC(2) : C122
**     CVC(3) : C113
**     CVC(4) : C133
**     CVC(5) : C223
**     CVC(6) : C233
**     CVC(7) : C123
**
********* THIS ROUTINE NEEDS SUBROUTINE MIXUP FOR NON IDEAL MIXING
**
********* AS OF FEB 86, BECAUSE OF LACK OF EXPERIMENTAL DATA TO
**        TO COMPUTE THESE VALUES, C113, C133 AND C123 ARE SET TO
**        ZERO.
**        DATA FOR H2O-CH4 IS EXTRAPOLATED ABOVE 100 DEG.C. DATA FOR
**        CO2-CH4 IS EXTRAPOLATED ABOVE 240 DEG.C.
********* WATCH ! THE CURRENT PRESSURE LIMIT FOR COMPUTING
**        FUGACITY COEFFICIENTS WITH ACCURACY AROUND 1-2% OR LESS
**        IS ABOUT 500 ATM FROM 150 TO 300 DEG.C, AND 200 ATM FROM
**        FROM 80 TO 150 DEG.C. IN THIS LOWER T RANGE, ERRORS CAN BE
**        UP TO 10% AT 500 BARS.  NOT RECOMMENDED BELOW
**        75 DEG.C. EXCEPT AT LOW P.  H2 DATA IS VALID UP TO
**        1300 BARS FROM 0 TO 600 DEG.C. A SET OF COEFFS VALID ABOVE
**        350 C UP TO 1000 BARS IS AVAILBALE FOR H2O-CO2.  SEE SPYCHER
**        AND REED, IN PREPARARION.
**
**    SEE SOLTHERM DATA FILE FOR SOURCE OF CURRENT DATA FOR PURE AND
**    MIXED GASES.
**
**    WRITTEN BY N.SPYCHER  FEB 86
**    **** LAST REVISED JAN 87. ADD DGDN TO EASE CONVERGENCE IN SOME
**         RARE CASES. N.S.
**
**    N1=1
**    IF(IDEAL.NE.0.OR.IMIX.EQ.0) GO TO 10
**    CALL MIXUP(VB,VC,CVB,CVC)
**    **********
**    N1 = IMIX + 1
**
**    IF IMIX IS 2, BVAL(3) AND CVAL(3) ARE OVERWRITTEN BELOW
** 10 CONTINUE
**    IDEAL MIXING
**    DO 15 I = N1,NGAS
**    BVAL(I) = VB(I)
** 15 CVAL(I) = VC(I)
**
**    NOW WE COMPUTE FUGACITY COEFFICIENTS
**    DO 20 I=1,NGAS
**    VV = BVAL(I)*PF + CVAL(I)*PF*PF/2.
**     SKIP IF VV IS TOO SMALL OR TOO LARGE AND KEEP OLD GAMIN VALUE
**    IF(DABS(VV).GT. 10.) GO TO 20
**    GAMIN(I) = DEXP(VV)
**    B(I) = 1.
** 20 CONTINUE
**
**    C-OUT BELOW AND FUNCTIONS ABOVE IF CONVERGENCE PROBLEMS ARISE,
**    PARTICULARTLY AT HIGH CO2 MOLE FRAC. IN GAS PHASE
**
**    NOW WE COMPUTE DERIVATIVES OF GAMIN W/ RESPECT TO Y (NON IDEAL
**    MIXING) TO EASE CONVERGENCE IN MINERA IN SOME RARE CASES.
**    SKIP IF: GASES ARE NOT SATURATED, OR IDEAL MIXING.
**    IF(ISOLS(1).EQ.0) RETURN
**    IF(IDEAL.NE.0.OR.IMIX.EQ.0) RETURN
**    SUMY = MINTRY(1) + MINTRY(2)
**    IF(IMIX.GT.2) SUMY = SUMY + MINTRY(3)
**    D1 = FUNC1(VB(1),CVB(1),CVB(2))
**    D2 = FUNC1(CVB(1),VB(2),CVB(3))
**    D3 = FUNC1(CVB(2),CVB(3),VB(3))
**    D12 = ( 2.*CVB(1) + 2.*BMIX - 2.*D1 - 2.*D2) * PF
**    D13 = ( 2.*CVB(2) + 2.*BMIX - 2.*D3 - 2.*D1) * PF
**    D23 = ( 2.*CVB(2) + 2.*BMIX - 2.*D1 - 2.*D3) * PF
**    E1 = FUNC2(VC(1),CVC(2),CVC(4),CVC(1),CVC(3),CVC(7))
**    E2 = FUNC2(CVC(1),VC(2),CVC(6),CVC(2),CVC(7),CVC(5))
**    E3 = FUNC2(CVC(3),CVC(5),VC(3),CVC(7),CVC(4),CVC(6))
**    E12 = (6.*FUNC1(CVC(1),CVC(2),CVC(7)) + 6.*CMIX - 6.*E1
**   &  - 6.*E2) * PF*PF/2.
**   E13 = (6.*FUNC1(CVC(3),CVC(7),CVC(4)) + 6.*CMIX - 6.*E1
**   &  - 6.*E3) * PF*PF/2.
**    E23 = (6.*FUNC1(CVC(7),CVC(5),CVC(6)) + 6.*CMIX - 6.*E2
**   &  - 6.*E3) * PF*PF/2.
**
**    DGDN(1,1) = GAMIN(1)/SUMY * (( 2.*VB(1) + 2.*BMIX - 4.*D1) * PF
**   &  + (6.*FUNC1(VC(1),CVC(1),CVC(3)) + 6.*CMIX - 12.*E1)*PF*PF/2.)
**    DGDN(1,2) = GAMIN(1)/SUMY * (D12 + E12)
**    DGDN(1,3) = GAMIN(1)/SUMY * (D13 + E13)
**
**    DGDN(2,1) = GAMIN(2)/SUMY * (D12 + E12)
**    DGDN(2,2) = GAMIN(2)/SUMY * (( 2.*VB(2) + 2.*BMIX -  4.*D2) * PF
**   &  + (6.*FUNC1(CVC(2),VC(2),CVC(5)) + 6.*CMIX - 12.*E2)*PF*PF/2.)
**    DGDN(2,3) = GAMIN(2)/SUMY * (D23 + E23)
**
**    IF(IMIX.LT.3) RETURN
**    DGDN(3,1) = GAMIN(3)/SUMY * (D13 + E13)
**    DGDN(3,2) = GAMIN(3)/SUMY * (D23 + E23)
**    DGDN(3,3) = GAMIN(3)/SUMY * (( 2.*VB(3) + 2.*BMIX - 4.*D3) * PF
**   &  + (6.*FUNC1(CVC(4),CVC(6),VC(3)) + 6.*CMIX - 12.*E3)*PF*PF/2.)
**    RETURN
**    END
**    SUBROUTINE MIXUP(B,C,CB,CC)
**    IMPLICIT DOUBLE PRECISION (A-H,O-Z)
**    INTEGER BG1,BG2,BG3,BG4
**    PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
**    COMMON/FUGA/GASlog(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3,
**   1 NGAS,IMIX,IDEAL,ITGAS
**    COMMON/VAL/BVAL(20),CVAL(20),BMIX,CMIX
**    DIMENSION B(20),C(20),CB(3),CC(7)
**
**    THIS ROUTINE COMPUTES BMIX, CMIX AND THEIR DERIVATIVE WITH
**    RESPECT TO WI AT T,P WJ NOT WI CONSTANT. THE RESULTING VALUES
**    BVAL AND CVAL ARE THEN USED IN PHI TO CALCULATE FUGACITY COEFFS.
**    SEE SUBROUTINE PHI FOR COMMENTS. FOR ENTHALPIES, WE PLUG IN
**    VALUES OF DB/DT, DC/DT, ETC. AND THEREFORE WE COMPUTES
**    DBMIX/DT, DCMIX/DT ETC. WE THEN USE RESULTING VALUES OF BVAL
**    AND CVAL TO COMPUTE ENTHALPIES OF GASES. SEE SUBROUTINE ENTHAL
**    B, C ARE 2ND AND 3RD VIRIAL COEFFS; CB, CC ARE 2ND AND 3RD CROSS
**    VIRIAL COEFFS AT GIVEN temp. (IF MIXUP IS CALLED FROM ENTHALP
**    THEN B ETC..  ARE -T**2 * DB/DT ETC...)
**    CROSS VIRIAL COEFFICIENTS ARE ASSIGNED AS FOLLOW AS IN SUBROUTINE
**    PHI.
**
**     WRITTEN BY NICOLAS SPYCHER, U.OF.O JAN 1986
**     LAST REVISED: FEB 86       BY: N.S
**
**    FUNC1(A1,A2,A3) = W1*A1 + W2*A2 + W3*A3
**    FUNC2(A1,A2,A3,A4,A5,A6) = W1*W1*A1 + W2*W2*A2 + W3*W3*A3
**   &    + 2.*( W1*W2*A4 + W1*W3*A5 + W2*W3*A6 )
**
**    WE HAVE TO NORMALIZE THE SUM OF THE MOLE FRACTIONS OF H2O,CO2,CH4
**    TO 1 BECAUSE WE IGNORE THE MIXING EFFECTS OF ALL OTHER GASES.
**
**    WTOT = Y(1) + Y(2)
**    W3=0.
**    IF(IMIX.LT.3) GO TO 5
**    WTOT=WTOT+Y(3)
**    W3=Y(3)/WTOT
**  5 W1=Y(1)/WTOT
**    W2=Y(2)/WTOT
**    IF(Y(1).LE.1.D-23) W1 = 0.
**    IF(Y(2).LE.1.D-23) W2 = 0.
**    IF(Y(3).LE.1.D-23) W3 = 0.
**    NON IDEAL MIXING
**
**    BMIX = FUNC2(B(1),B(2),B(3),CB(1),CB(2),CB(3))
**    CMIX = W1*W1*W1*C(1) + W2*W2*W2*C(2) + W3*W3*W3*C(3)
**   &     + 3.*( W1*W1*W2*CC(1)
**   &     + W1*W2*W2*CC(2) + W1*W1*W3*CC(3)  +  W1*W3*W3*CC(4)
**   &     + W2*W2*W3*CC(5) + W2*W3*W3*CC(6) + 2.*W1*W2*W3*CC(7) )
**    BVAL AND CVAL ALSO NEEDED IN PSAT TO SET UP PARTIAL DERIVATIVES
**    OF GAMIN WITH RESPECT TO PRESSURE  (TO SOLVE FOR SAT. PRESSURE)
**
**    Z   = 1 + BMIX*P + CMIX*P*P
**    BVAL(1) = 2.*FUNC1(B(1),CB(1),CB(2)) - BMIX
**    CVAL(1) = 3.*FUNC2(C(1),CC(2),CC(4),CC(1),CC(3),CC(7) )
**   &        - 2.*CMIX
**
**    BVAL(2) = 2.*FUNC1(CB(1),B(2),CB(3)) - BMIX
**    CVAL(2) = 3.*FUNC2(CC(1),C(2),CC(6),CC(2),CC(7),CC(5) )
**   &        - 2.*CMIX
**
**    BVAL(3) = 2.*FUNC1(CB(2),CB(3),B(3)) - BMIX
**    CVAL(3) = 3.*FUNC2(CC(3),CC(5),C(3),CC(7),CC(4),CC(6) )
**   &        - 2.*CMIX
**    IF IMIX IS 2, BVAL(3) AND CVAL(3) WILL BE OVERWRITTEN
**    RETURN
**    END
**    SUBROUTINE ENTHAL(TC)
**
**    IMPLICIT DOUBLE PRECISION (A-H,O-Z)
**    INTEGER BG1,BG2,BG3,BG4
**    PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85)
**    COMMON/FROID/ SINC,BIT,SLIM,AKlog,APFCL(BG2,5),
**   1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4),
**   2 XDH(4),ALK(5),LGK(BG2),temp,tempC, HEATC,COMTOT(BG3), HEATM,
**   3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3),
**   4 CVCCF(7,3), IFRAC
**    COMMON/FUGA/GASlog(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3,
**   1 NGAS,IMIX,IDEAL,ITGAS
**    COMMON/VAL/BVAL(20),CVAL(20)
**    DIMENSION WATCF1(4),WATCF2(5),DBDT(20),DCDT(20),DCBDT(3),DCCDT(7)
**    DOUBLE PRECISION LGK
**    DATA WATCF1/2.335628,-.00441813,2.30965D-5,-777.40446/,
**   1  WATCF2/240.912149,-.35584278,.00021415,-7.195396D+4,7.895231D-6/
**    DATA A1,B1,C1/7.17,2.56D-3,.08D+5/,
**   1     A2,B2,C2/10.55,2.16D-3,-2.04D+5/,
**   2     A3,B3,C3,D3/2.975,18.329D-3,.346D+5,-4.303D-6/
**    FUNCTION DHT FROM CP = A + B*T + C/T**2
**    DHT(A,B,C,T) = A*(T-298.15) + B/2.*(T*T-298.15*298.15) -
**   1               C*(1./T - 1./298.15)
**
**    THIS SUBROUTINE CALCULATES ENTHALPIES OF WATER, H2O GAS,
**    CO2 AND CH4. WATER ENTHALPY IS ASSUMED  TO BE PRESSURE
**    INDEPENDENT.
**      THE ENTHALPY OF WATER FOR A GIVEN   tempERATURE   IS
**      CALCULATED AT THE CORRESPONDING SATURATION PRESSURE. SINCE
**      WATER IS NEARLY INCOMPRESSIBLE, THE DIFFERENCE IN ENTHALPY
**      AT PRESSURES HIGHER THAN THE SATURATION PRESSURE IS NEGLIGIBLE.
**      THE ENTHALPY IS CALCULATED FROM A POWER FUNCTION OF  THE
**      tempERATURE.
**            H(TK) = A + B*TK + C*TK**2 + D/TK + E/TK**2
**      WHERE TK IS THE tempERATURE IN DEG. KELVIN. THE COEFFICIENTS
**      A,B,C,D,E WERE LEAST-SQUARE FITTED FROM HELGESON AND KIRKHAM
**      1974 PART I (AJS, 274, TABLE 37,P.1184). THE VALUES OF THESE
**      COEFFICIENTS ARE SUCH THAT THE CALCULATED ENTHALPY H(TK) IS
**      IN KCAL/MOLE AND IS ZERO AT 0 DEG.C.
**    ENTHALPY OF H2O, CO2 AND CH4 GAS:
**      THE ENTHALPIES OF THESE GASES ARE CALCULATED FROM  HEAT
**      CAPACITY DATA (BARIN&KNACKE,1973,1977,THERMOCHEMICAL PROPERTIES
**      OF INORGANIC SUBSTANCES, SPRINGER-VERLAG) AND FROM VIRIAL
**      COEFFICIENTS (SAME DATA THAT WE USE TO COMPUTE FUGACITY
**      COEFFICIENTS, SEE SUBROUTINE PHI). THE ENTHALPY IS THEN:
**          H GAS MIXTURE = SUM OF (MOLES(I) * (DH/DT (I) + DH/DP (I)))
**      IF IDEAL MIXING, THEN DH/DP IS FOR PURE GAS. IF NO-IDEAL MIXING
**      THEN DH/DP IS FOR THE GAS IN THE MIXTURE.
**      DHT IS THE CHANGE OF ENTHALPY WITH tempERATURE, CALCULATED FROM
**      A MAIER-KELLY HEAT CAPACITY POWER FUNCTION, AND DHP IS THE
**      CHANGE WITH PRESSURE, CALCULATED FROM THE FOLLOWING EQUATION:
**          DHP = R*TK**2 * ( DB/DT*P + DC/DT*P**2/2 )
**      WHERE B AND C ARE THE SECOND AND THIRD VIRIAL COEFFICIENTS
**      OF THE VIRIAL EQ. OF STATE IN PRESSURE. FOR NON IDEAL MIXTURES
**      WE TAKE THE APPROPRIATE PARTIAL DERIVATIVES.
**      ****** WARNING:  DHP MAY BE UP TO  8% IN ERROR AT HIGHER T,P
**      HOWEVER THE ERROR IN THE TOTAL ENTHALPY OF THE GASES DOES NOT
**      EXEED A FEW PRECENT. MAKE SURE TO RESPECT T,P LIMITS (SEE
**      SUBROUTINE PHI).
**      ENTHAPIES OF CO2 AND CH4 ARE ZERO AT 0 DEG.C. AND 0 BAR.
**      FOR H2O GAS, AT 0 DEG.C. 0.006 BAR, H = ENTH. OF VAPORIZATION.
**      (10.77 KCAL/MOLE)
**      WRITTEN BY NICOLAS SPYCHER, U.OF.O FEB 86
**      LAST REVISED:           BY:
**
**    TK=TC + 273.15
**    GETS HEAT OF WATER IF LESS THAN 100 DEG.C. COEFFICIENTS
**    WATCF1 VALID ONLY FROM 0 TO 99 DEG.C.
**    IF(temp.GE.100.) GO TO 20
**    ENTHW = WATCF1(1) + WATCF1(2)*TK + WATCF1(3)*TK*TK + WATCF1(4)/TK
**    GETS HEAT OF WATER IF 100 DEG.C.OR ABOVE. WATCF2 COEFFICIENTS
**    VALID FROM 100 TO 350 DEG.C.
** 20 CONTINUE
**    ENTHW = WATCF2(1) + WATCF2(2)*TK + WATCF2(3)*TK*TK + WATCF2(4)/TK
**   & + WATCF2(5)/TK/TK
**    RETURN
**
**    ENTRY PFENTH(TC,P)
**    ***************
**    TK=TC + 273.15
**    GET HEAT OF GASES. 1 H2O; 2 CO2; 3 CH4;
**    ENTHM(1)= DHT(A1,B1,C1,TK)/1000. +  10.970
**    ENTHM(2)= DHT(A2,B2,C2,TK)/1000. + .210
**    ENTHM(3)= DHT(A3,B3,C3,TK)/1000. + D3*(TK**3.-298.15**3.)/3000.
**   1    + .210
**
**    IF(IDEAL-1) 53,52,51
** 51 RETURN
**    PREF IS ONE ATM SINCE CP DATA IS AT ONE ATM
**    IDEAL MIXING/REAL GAS WE GET DHP FOR PURE GASES THEN ADD TO DHT
** 52 CONTINUE
**    DO 50 I=1,IMIX
**     DB = 2.*VBCOE(I,1)/TK + VBCOE(I,2)
**     DC = 2.*VCCOE(I,1)/TK + VCCOE(I,2)
**     DHDP(I) = 1.987 * (DB + DC*P)/1000.
**     DHP = 1.987 * (DB*(P-1.)  + DC*(P*P-1.)/2.)/1000.
** 50  ENTHM(I) = ENTHM(I) + DHP
**    RETURN
**
** 53 CONTINUE
**    IF(Y(1).EQ.0.) GO TO 52
**    NON IDEAL MIXING/REAL GAS. WE CALCULATE DHP FOR GAS IN MIXTURE
**    THEN ADD TO DHDT. PREF IS ONE ATM SINCE CP DATA IS AT 1 ATM.
**     BELOW WE GET VALUES OF T**2 * DB/DT (ALSO FOR C AND CROSS B,C)
**    DO 60 I = 1,IMIX
**     DBDT(I) = 2.*VBCOE(I,1)/TK + VBCOE(I,2)
** 60  DCDT(I) = 2.*VCCOE(I,1)/TK + VCCOE(I,2)
**    DO 62 I = 1,3
** 62  DCBDT(I) = 2.*CVBCF(I,1)/TK + CVBCF(I,2)
**    DO 64 I = 1,7
** 64  DCCDT(I) = 2.*CVCCF(I,1)/TK + CVCCF(I,2)
**
**    CALL MIXUP(DBDT,DCDT,DCBDT,DCCDT)
**    **********
**    DO 66 I=1,IMIX
**     DHDP(I) = 1.984* (BVAL(I) + CVAL(I)*P)/1000.
**     DHP= 1.984* (BVAL(I)*(P-1.) + CVAL(I)*(P*P-1.)/2.)/1000.
** 66  ENTHM(I) = ENTHM(I) + DHP
**    RETURN
**    END
