C    ******************IBM PC OR VAX VERSION****************************
C    *******************************************************************
C    **                                                               **
C    **                 SOLVGAS                                       **
C    **           (IBM PC VERSION 1.0)                                **
C    **                                                               **
C    **       (1) MARK H REED (UNIVERSITY OF OREOGN)                  **
C    **       (2) ROBERT B. SYMONDS (MICHIGAN TECH.)                  **
C    **  THIS COPY SENT TO JIM WOOD                                   **
C    *******************************************************************
C    *******************************************************************
C     THIS PROGRAM WAS MODIFIED FROM SOLVEQ(M. REED, U.C. BERKELEY,1976)
C     MODIFICATION OF SOLVEQ BEGINNING FEB. 1986 BY M. REED (UNIV. OF
C     OREGON) AND R. SYMONDS (MICHIGAN TECH.)
C    *******************************************************************
C     THIS PROGRAM NEEDS DATA FILE GASTHERM TO RUN
C
C     THIS VERSION IS DESIGNED TO RUN WITH THE GASTHERM DATA FILE,
C     WHICH CONTAINS LOG K'S AS WELL AS REGRESSION COEFFICIENTS, WHICH
C     ALLOW CALCULATIONS AT ANY GIVEN TEMPERATURE WHEN NTEMP IS SET
C     TO ZERO.
C     ***************REVISIONS OF SOLVGASQ FORTRAN**********************
C
C      1). AUGUST 87: N. SPYCHER ELIMINATED SOME BUGS, SIMPLIFIED SOME
C          NUMERICAL EXPRESSIONS AND CHANGED CONVERGENCE CRITERION TO
C          F(I) VALUES.
C      2). NOV. 87: R. SYMONDS ADDED PLOTTING FILES AT OUTPUT, MADE
C          UNACT AN INTERACTIVE VARIABLE AND FIXED A CONVERGENCE ERROR.
C      3). JUL. 88: MH. REED AND M. HIRSCHMANN SWAPPED O2 FOR H2O AS
C          DERIVED COMPONENTS, CHANGED CONVERGENCE TEST, AND ALTERED
C          FORMAT FOR READING IN GASRUN  FILES.
C      4). SEPT 88: M. REED SWATTED GNATS IN OXY AND WORKED ON THE
C          "DO 358" ZONE
C      5). OCT - NOV 88: R SYMONDS WOKED ON CONVERGENCE CRITERIA.  ADDED
C          MANY IDEAS FROM SUPER SOLVEQ TO AID CONVERGENCE.  ALSO ADDED
C          ALKALN AS A FUDGE FACTOR TO STOP MTRYS FOR COMPONENT AND
C          DERIVED SPECIES FROM GETTING TOO LARGE.  ALSO INITIALIZED
C          SOME ARRAYS AND SWATTED OTHER GNATS.
C      6). DEC 88: M REED, REORDERED TYPE STATEMENTS, CHANGED TO FO2LOG
C          INPUT, CORRECTED A 358 ZONE BUG.
C      7). FEB 89: R SYMONDS ADDED PRESSURE ROUTINES TO CALCULATE MOLE
C          FRACTIONS AS A FUNCTION OF PRESSURE.  ALSO CORRECTED
C          AN ERROR IN RESETTING UNACT.
C       8) NOV 90: R SYMONDS CONVERTED PROGRAM TO RUN WITH NEW
C          GASTHERM.
C       9) NOV 90: R SYMONDS CONVERTED PROGRAM TO FORTRAN 77
C      10) MARCH 90: R SYMONDS CONVERTED PROGRAM TO RUN ON 386 PC'S WITH
C          LAHEY FORTRAN.
c      11) May 2002, Changed I/O device numbers (MR)
C     ******************************************************************
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      DOUBLE PRECISION LGK, MTRY, K, MTOT
      INTEGER SAQ, QAS, SAQTOT, ERLOOP, SPEC, BSTOT, BSTP1, BBST,
     1 BSTP2
      CHARACTER*1 ANS, POS, pos2, STOP,stop2, SUP, BLNK
      CHARACTER*2 STA1,STA2,STA3,STAR(50)
      CHARACTER*8 ETOIL,TEXT, TITLE1, TITLE2, NAME
      DIMENSION TEMP(13),CONV(50)
      DIMENSION NAME(999), TITLE1(10),TITLE2(10),W(50),
     1  ZILCH(13), MTOT(50), X(50)
      COMMON/SLOUGH/A(50,51), T(50,51), IV(50),N
      COMMON/MINNIE/ ACT(999), TEMPC, NTEMP, NMTT
      COMMON/DRUNK/COEF(999,7), K(999),LGK(999),SPEC(999,7),IT(999),
     1 BSTP1,BSTOT
      COMMON/HUH/ UNACT, TERMLN,NSUNK
      COMMON/IO/IPRINT,IN1,IN2,IN3,IOUT1,IOUT2,IOUT3,IOUT4,ITREF
      COMMON/MORT/ MTRY(999),GAMMA(999),FO2,F(50),KK
      COMMON/EREBUS/ PF,GASTOT,UNALOG,ALA1
      COMMON/GEORGE/ AZERO(999),
     1NS,NST
      COMMON/TURNIP/ SAQ(50), QAS(50), SAQTOT
      DATA BLNK/' '/
      DATA   STA1  /'  '/
      DATA  STA2 /'CM'/
      DATA   STA3  /'CM'/
      DATA TEMP/25.,100.,200.,300.,400.,500.,600.,700.,800.,900.,1000.,
     11100.,1200./
      DATA POS/'Y'/,pos2/'y'/,STOP/'S'/,stop2/'s'/,ETOIL/'********'/
C
C     THESE ASSIGNMENTS DEFINE WHAT DEVICE NUMBERS ARE USED IN THE
C     READ AND WRITE STATEMENTS:
C     IN1 IS GASRUN, IN2 IS GASTHERM, IOUT1 IS GASOUT
C       INPUT: - IN1    DATA SPECIFIC TO THE CALCULATION (FO2,MTOT,ETC.)
C              - IN2    COMPLEXES, GASES AND MINERAL DATA
C              - IN3    TERMINAL DISPLAY FOR INTERACTIVE RUNS
C       OUTPUT: - IOUT1    MAIN PROGRAM OUTPUT
C               - IOUT2    TERMINAL DISPLAY FOR INTERACTIVE RUNS
C               - IOUT3    FILE TO PLOT LOG(FUGACITIES)
C               - IOUT4    FILE TO PLOT LOG(Q/K)'S
C
      IN1=4
      IN2=1
      IN3=5
      IOUT1=9
      IOUT2=6
      IOUT3=2
      IOUT4=7
C     THESE OPEN STATEMENTS ASSIGN EXPLICIT, 'DOS-APPROVED', FILENAMES
C     TO THE DEVICE NUMBERS DEFINED ABOVE. USERS WORKING IN MAINFRAME
C     ENVIRONMENTS, WHO WOULD LIKE TO HAVE MORE FLEXIBILTY IN FILE
C     MANAGEMENT, CAN SIMPLY DISABLE THESE OPEN STATEMENTS AND MAKE
C     THEIR DEVICE ASSIGNMENTS AT RUN-TIME.
      OPEN (UNIT = 1, FILE = 'GASTHERM.DAT', STATUS = 'OLD')
      OPEN (UNIT = 2, FILE = 'GASPLOT.DAT',  STATUS = 'UNKNOWN')
      OPEN (UNIT = 4, FILE = 'GASRUN.DAT',   STATUS = 'OLD')
      OPEN (UNIT = 9, FILE = 'GASOUT.DAT',   STATUS = 'UNKNOWN')
      OPEN (UNIT = 7, FILE = 'SUBPLOT.DAT',  STATUS = 'UNKNOWN')
C
C     ARRAY MAXIMUM SIZES: MAXSAQ = MAXIMUM NUMBER OF COMPONENT SPECIES;
C     MAXNS = MAXIMUM NUMBER OF GAS SPECIES ALLOWED;
C     MAXSTO = MAXIUM NUMBER OF STOICHIOMETRIC COMPONENTS
C     MAXNM = MAXIMUM NUMBER OF SOLIDS AND LIQUIDS
      MAXSAQ = 50
      MAXNS = 999
      MAXNM = 900
      MAXSTO = 8
C
C     COMPUTER PRECISION LIMITS TO AVOID OVERFLOWS, UNDERFLOWS, ETC.
C     FOR IBM 3281, SMALLEST DOUBLE PRECISION NUMBER IS ABOUT 10**-70,
C     LARGEST IS ABOUT 10**70.  FOR VAX AND 386 PERSONAL COMPUTER
C     SYSTEMS, THESE are 10**-300 and 10**300.  RMIN: MININMUM REAL #;
C     RLOGMN: MINIMUM LOG10 REAL #; RLGMN2: MINIMUM LOG10 REAL # -2
      RMIN = 1.D-300
      RLOGMN = -300.
      RLGMN2 = -302.
C
      FO2=0.
      FO2LOG = 0.
      UNACT = 0.
      UNALOG = 0.
      GASTOT=0.
      GSTOLD=0.
      BASE = 10.
      TERMLN = 1./DLOG(BASE)
C
      READ(IN1,870) TITLE1
      READ(IN1,870) TITLE2
C     **************
C     TITLE1 AND TITLE2 CAN BE ANY TEXT STRINGS USED TO IDENTIFY THE
C     GAS ANALYSIS
      WRITE(IOUT3,871) TITLE1,TITLE2
      WRITE(IOUT4,871) TITLE1,TITLE2
C
      READ(IN1,943)
      READ(IN1,833) ERPC, FO2LOG, UNALOG, ALKALN, TEMPC, PF,RAT
C     **************
      IF(FO2LOG.NE.0.) FO2 = 10.**FO2LOG
      IF(UNALOG.NE.0.) UNACT = 10.**UNALOG
C     ERPC IS THE FACTOR FOR DETERMINING NEWTON-RAPHSON CONVERGENCE.
C     USE E-12 OR SMALLER
C     UNACT CONTAINS THE FUGACITY OF ANOTHER COMPONENT SPECIES, SUCH
C       AS H2S, HCL OR HF.
C     ALKALN IS A FUDGE FACTOR TO AID CONVERGENCE BY PREVENTING MTRYS
C        OF COMPONENT AND DERIVED SPECIES FROM GETTING TOO LARGE
C        NORAMALLY, ONE SHOULD SET THIS EQUAL TO THE APPOXIMATE GASTOT.
C     TEMPC :     THE USER CAN ENTER ANY TEMPERATURE. THE LOG K'S WILL
C       BE CALCULATED FROM THE REGRESSION COEFFICIENTS (SLOWER PROCESS
C       THAN WHEN READING LOG K'S DIRECTLY. SEE BELOW)
C     PF = THE FLUID PRESSURE IN ATMOSPHERES
C     FO2LOG = LOG OF OXYGEN FUGACITY.  USUALLY FO2 WILL BE ENTERED FOR
C     GIVEN GAS ANALYSIS AND HELD CONSTANT.  SET FO2LOG = 0. IF YOU WANT
C     THE PROGRAM TO CALCULATE IT.
C     RAT: IS THE LOG OF THE FUGACITY BELOW WHICH SPECIES ARE NOT
C          PRINTED IN THE OUTPUT
      READ(IN1,943)
      READ(IN1,800) BSTOT, NTEMP, NMTT, NMOLF, IGAMP, NLOOP, NSUNK,
     1IPRINT,ITREF
C     ******************************************************************
C     ITREF = 0    : NO ITERATIVE REFINEMENT
C           = 1    : NO ITERATIVE REF. BUT CALCULATION OF RESIDUAL
C           > 1    : 1 OR MORE ITERATIVE REFINEMENTS
C     NSUNK IS THE SAQ-TYPE NUMBER CORRESPONDING TO UNACT
C     NTEMP= 1 TO 13 IS FOR CALCULATIONS WITH LOG K'S LISTED IN
C       POSITIONS 1 TO 13. THESE ARE CURRENTLY LOG K (25) (100)
C       (200) (300) (400) (500) (600) (700) (800) (900) (1000)
C       (1100) (1200) DEG. C.  TEMPC NEEDS TO BE ZERO FOR THIS OPTION.
C     BSTOT IS THE NUMBER OF COMPONENT SPECIES IN THE INPUT FILE
C      (GASTHERM)
C     NMTT SHOULD BE NOT ZERO IF MINERALS ARE TO BE
C       READ FROM INPUT FILE (SOLTHERM). IF CALCULATION OF Q/K FOR
C       MINERALS IS NOT DESIRED, SET NMTT=0.
C     NMOLF IS USED FOR CHANGING PRESSURE.  IF YOU WANT TO INVOKE THIS
C        OPTION, SET NMOLF = 1 (OTHERWISE SET TO 0).  NOTE THAT
C        THIS WILL CAUSE LOG MOLE FRACTIONS (NOT LOG FUGACITIES) TO
C        BE PRINTED IN THE GASPLOT DATA FILE.
C     IGAMP NOT ZERO MEANS FUGACITY COEFFICIENTS ARE TO BE READ IN
C      BELOW FOR ALL SPECIES.  IF YOU'VE NEVER RUN THIS GAS BEFORE
C      IT IS EASIEST TO SET IGAMP = 0
C     NLOOP SETS A LIMIT ON THE NUMBER OF NEWTON ITERATIONS THAT
C        WILL BE TRIED.  40 IS A SAFE LIMIT.
C     BE SURE TO CHECK THAT THE ACTUAL NUMBER OF LOOPS USED WAS LESS
C      THAN NLOOP; IF NOT, INCREASE NLOOP (E.G. TO 60 OR 80) AND RERUN
C     IPRINT NOT ZERO RESULTS IN EXTRA PRINTING OF RESULTS
C
C      INITIALIZE ARRAYS BEFORE USING THEM
      I=0
      DO 141   I = 1,MAXSAQ
      SAQ(I) =   0
      MTOT(I) =  0.
      GAMMA(I) = 0.
  141 QAS(I)   = 0
  142 CONTINUE
C
C     INPUT FORMAT THAT FOLLOWS HAS BEEN ADJUSTED TO READ NEW
C     (SOLRUN-LIKE) GASRUN FILES, MH JULY, 1988
      I=1
      READ(IN1,943)
  100 READ(IN1,805) SAQ(I), MTOT(I), MTRY(I),GAMMA(I)
  805 FORMAT (I5,13X,2E16.9,F8.4)
C     *********************
C
C     SAQ(I) SPECIFIES WHICH OF THE COMPONENT SPECIES ARE TO BE INCLUDED
C     E.G. FOR H2, H20, CO2, H2S, HCL, HF AND HBR, PUT 1 2 3 4 5 6 7
C     THESE NUMBERS ARE THE SEQUENCE POSITIONS OF INPUT CARDS FOR
C     BASIS SPECIES AS READ IN DO 110 BELOW (FILE GASTHERM)
C     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
      IF(J.EQ.NSUNK) NSUNK=I
      IF (MTRY(I).LE.0.)  MTRY(I) = 1.E-10
      IF (GAMMA(I).LE.0.)  GAMMA(I) = 1.
      I = I + 1
      IF(I.LE.MAXSAQ) GO TO 100
      READ(IN1,805) IDUM
C      (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 GASRUN
     1 EXCEEDS THE '/5X, 'CURRENT ARRAY DIMENSION (',I2,').  CHECK',
     2 ' GASRUN 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
C     THE DO 123 LOOP THAT FOLLOWS IS TO ESTIMATE THE INITIAL NUMBER
C     OF MOLES OF GAS.  GASTOT IS THE TOTAL NUMBER OF MOLES OF GAS.
  144 SAQTOT = I - 1
      DO 123 I = 1, SAQTOT
C     **************
  123 GASTOT = GASTOT + MTOT(I)
      GSTOLD = GASTOT
C     ********************
C     MTRY(I) ARE TRIAL MOLE NUMBERS FOR ALL COMPONENT SPECIES. IF THIS
C       IS THE FIRST RUN FOR THIS GAS ANALYSIS, THEN USE THE SAME
C       VALUES AS THE GAS ANALYSIS AND ESTIMATE FO2 BY USING THE QFM
C       (WONES AND GILBERT, 1969) OR NI-NIO BUFFER AT THE TEMPERATURE
C       OF INTEREST.
      BBST=BSTOT
      IF(TEMPC.EQ.0.) WRITE(IOUT2,780) TEMP(NTEMP),FO2LOG,PF
      IF(TEMPC.NE.0.) WRITE(IOUT2,780) TEMPC,FO2LOG,PF
      WRITE(IOUT2,782)
      GO TO 30
C
C       SOLVGAS INTERACTIVE SECTION FROM HERE TO FORTRAN LINE 30
C       *******************************************************
C
C     THE FO2LOG IS SET EQUAL TO ZERO AND MTOT(H2) STAYS THE SAME
C     AS THE ONE CALCULATED IN THE FIRST RUN WITH GIVEN FO2.
C     MTOT(H2) BECOMES THEREFORE AN INPUT DATUM AND THE FO2LOG
C     AN OUTPUT RESULT
C
C     IF YOU INTEND TO RUN THIS PROGRAM IN BATCH MODE, YOU
C     WILL HAVE TO INSERT A "STOP" STATEMENT AFTER THE
C     "50 CONTINUE" LINE BELOW. IF YOU FORGET, THE COMPUTER
C     WILL COME BACK TO THIS STATEMENT AFTER THE FIRST
C     EXECUTION AND WAIT FOR AN ANSWER TO CONTINUE OR NOT.
C
   50 CONTINUE
C     STOP CARD HERE FOR BATCH JOBS !!!!!!!!!!!!!!!!!!!!!!!!!!
      REWIND 1
      BSTOT=BBST
      IF (FO2LOG.EQ.0.) FO2LGP = DLOG10(ACT(BSTP1))
      IF (FO2LOG.NE.0.) FO2LGP = FO2LOG
      WRITE(IOUT2,701) FO2LGP
      IF (UNACT.NE.0.) WRITE(IOUT2,702) NAME(NSUNK), UNALOG
      WRITE(IOUT2,965) ERLOOP, NLOOP
      WRITE(IOUT2,700)
      READ(IN3,710) ANS
      IF(ANS.NE.POS.and.ans.ne.pos2) GO TO 19
C     RESETTING FO2 = 0 AND UNACT= 0 SO PROGRAM WILL CALCULATE THEM AS
C     TEMPERATURE IS CHANGED
      FO2 = 0.
      FO2LOG = FO2
      UNACT = 0.
      UNALOG = UNACT
      NSUNK = 0
C
      IF(TEMPC.NE.0.) GO TO 15
   14 IF(NTEMP.NE.0) DEGC=TEMP(NTEMP)
      IF(TEMPC.NE.0.) DEGC=TEMPC
      WRITE(IOUT2,720) DEGC
      READ(IN3,730) NTEMP
      GO TO 17
C
   15 WRITE(IOUT2,500) TEMPC
  500 FORMAT(///5X,'THE PREVIOUS TEMPERATURE WAS: ',F6.0,/
     1 5X,'INPUT  THE NEW TEMPERATURE IN DEG. C.',/5X,'OR ZERO TO '
     2, 'CHOOSE FROM PRESELECTED TEMPERATURES'/)
      READ(IN3,*) TEMPC
      IF(TEMPC.EQ.0.) GO TO 14
C
   17 CONTINUE
      WRITE(IOUT2,740)
      WRITE(IOUT3,942)
      WRITE(IOUT4,942)
      GO TO 30
   19 IF(NMOLF.NE.1) GO TO 20
      WRITE(IOUT2,745)
      READ(IN3,710) ANS
      IF(ANS.NE.POS.and.ans.ne.pos2) GO TO 20
      WRITE(IOUT2,746) PF
      READ(IN3,771) PF
C     RESET FO2 AND UNACT HERE
      FO2 = 0.
      FO2LOG = FO2
      UNACT = 0.
      UNALOG = UNACT
      NSUNK = 0
      GO TO 17
   20 IF(UNACT.EQ.0.) GO TO 21
      WRITE(IOUT2,749) NAME(NSUNK)
      READ(IN3,710) ANS
      IF(ANS.NE.POS.and.ans.ne.pos2) GO TO 21
      WRITE(IOUT2,759) NAME(NSUNK), UNALOG
      READ(IN3,770) UNALOG
      UNACT = 10.**(UNALOG)
C     RESET FO2 HERE
      FO2 = 0.
      FO2LOG = FO2
      GO TO 17
   21 WRITE(IOUT2,750)
      READ(IN3,710) ANS
      IF(ANS.EQ.STOP.or.ans.eq.stop2) STOP
C     RESET UNACT HERE
      UNACT = 0.
      UNALOG = UNACT
      NSUNK = 0
      IF(ANS.NE.POS.and.ans.ne.pos2) GO TO 17
      WRITE(IOUT2,760) FO2LGP
      READ(IN3,770) FO2LOG
      FO2 = 10.**(FO2LOG)
      GO TO 17
  700 FORMAT (////5X,'DO YOU WANT TO CHANGE THE TEMPERATURE ?',
     11X,'(YES/NO)')
  701 FORMAT(////5X,'THE LOG OF THE OXYGEN FUGACITY IS:  ',F12.8)
  702 FORMAT(5X,'THE LOG OF THE ',A8,'FUGACITY IS:  ',F12.8)
  710 FORMAT(A1)
  720 FORMAT(5X,'THE OLDER TEMPERATURE WAS: ',F5.0,' DEGREES C.',//,5X,
     1'TO INPUT THE NEW TEMPERATURE TYPE (I2):',/,7X,
     2' 1    2     3     4     5     6      7     ',/,24X,'FOR',/,7X,
     3'(25)(100) (200) (300) (400) (500)  (600)  DEG. C',//,7X,
     4'  8     9     10    11     12     13       ',/,24X,'FOR',/,7X,
     5'(700) (800) (900) (1000) (1100) (1200)    DEG. C',//,7X,
     6'INPUT TWO DIGITS, E.G. "07" FOR "7"',/)
  730 FORMAT(I2)
  740 FORMAT(5X,'OKAY. THE SPECIES DISTRIBUTION IS GOING TO BE',
     1'CALCULATED')
  745 FORMAT(//5X,'DO YOU WANT TO CHANGE THE PRESSURE (YES/NO)?',///)
  746 FORMAT(//5X,'THE OLD PRESSURE WAS: ',F12.3,' ATMOSPHERES',//,5X,
     1'TYPE THE NEW ONE (F12.3)',///)
  749 FORMAT(//5X,'DO YOU WANT TO CHANGE THE ',A8,'FUGACITY (YES/NO)?')
  750 FORMAT(//5X,'DO YOU WANT TO CHANGE THE FO2 (YES), OR STOP?',
     1'(STOP)'///)
  759 FORMAT(//5X,'THE OLD LOG ',A8,'WAS: ',F12.8,//,5X,
     1'TYPE THE NEW ONE (F12.8)',///)
  760 FORMAT(//5X,'THE OLD LOG FO2 WAS: ',F12.8,//,5X,
     1'TYPE THE NEW ONE (F12.8)',///)
  770 FORMAT(F12.8)
  771 FORMAT(F12.8)
  780 FORMAT(//,5X,'FIRST RUN STARTS...',//,8X,'ORIGINAL TEMPERATURE:',
     1F5.0,'  DEGREES C.',/,8X,'ORIGINAL LOG FO2: ',F12.8,/,
     18X,'ORIGINAL PRESSURE (ATMOSPHERES): ',F12.3,//)
  782 FORMAT(/,5X,'OOOOOWEEEE...HOT GASES...',/,5X,
     1'  SYMONDS SAYS THEY ARE GREAT!!!',/)
C
C
   30 CONTINUE
C     STARTS READING GASTHERM DATA FILE. EITHER STRAIGHT LOG K'S OR
C     REGRESSION COEFFICIENTS FOR LOG K AS F(TEMP) ARE READ, DEPENDING
C     ON THE SELECTED OPTION.
C
C     THE FOLLOWING SKIPS THE COMMENTS AT THE TOP OF THE DATA FILE
C     UNTIL A LINE OF STARS IS REACHED (8 STARS MINIMUM, FROM 1ST COL.)
   60 READ(IN2,820) TEXT
      IF(TEXT.NE.ETOIL) GO TO 60
C
C
C
C
C
      BSTP1 = BSTOT + 1
      IF(ERPC.EQ.0.) ERPC = 1.E-12
      DO 110 NJ = 1, BSTOT
      IF(QAS(NJ).NE.0)GO TO 103
      READ(IN2,820) XILCH
C     BYPASS SPECIES NOT SPECIFIED BY SAQ
      GO TO 110
  103 CONTINUE
      NS = QAS(NJ)
      READ(IN2,822) NAME(NS), W(NS)
C     **************
C     READ CHARACTERISTICS OF COMPONENT SPECIES
C     W(NS) CONTAINS MOLECULAR WEIGHT OF SPECIES
  110 CONTINUE
      NS = SAQTOT
  113 CONTINUE
      NS = NS + 1
      IF (NS.LE.MAXNS) GO TO 114
      NS = NS - 1
      WRITE (IOUT1,660) MAXNS, NAME(NS), NS
      WRITE (IOUT2,660) MAXNS, NAME(NS), NS
  660 FORMAT(/5X,'NUMBER OF SPECIES READ IN GASTHERM EXCEEDS MAXIMUM',
     1' OF ',I3,/5X,'LAST NAME READ IS:  ',A8,' (NS=',I3,')'//)
      STOP
  114 READ(IN2,820) NAME(NS), SUP, ITOT, (COEF(NS,I),
     1SPEC(NS,I), I = 1,ITOT)
C     **************
C
C     COEF(NS,I) CONTAINS THE STOICHIOMETRIC COEFFICIENT FOR SPECIES
C       "SPEC(NS,I)" IN THE DISSOCIATION EQUATION FOR THE DERIVED GAS
C     SPEC(NS,I) IS A SAQ-TYPE NUMBER IDENTIFYING THE COMPONENT SPECIES
C       THAT GOES WITH COEF(NS,I)
C     THE I = 1 POSITION IS LEFT BLANK.
C     IF H2 IS IN THE DERIVED GAS IT SHOULD BE IN THE I=2 POSITION.
      IF(ITOT.NE.0) GO TO 115
      NST = NS - 1
      GO TO 120
  115 CONTINUE
      IT(NS) = ITOT
      IF(ITOT.LE.MAXSTO) GO TO 111
      WRITE (IOUT1,662) NAME(NS), ITOT, MAXSTO
      WRITE (IOUT2,662) NAME(NS), ITOT, MAXSTO
  662 FORMAT(/5X,'SPECIES ',A8,' IN GASTHERM HAS TOO MANY COMPONENTS',
     1/5X,'ITOT FOR THAT SPECIES IS ',I3,' ; MAXIMUM IS ',I3//)
      STOP
C
  111 IF(TEMPC.NE.0.) GO TO 150
      IF(NTEMP.GT.8) GO TO 148
      READ(IN2,860) (ZILCH(I), I = 1,NTEMP)
      READ(IN2,860) DUM
      GO TO 149
  148 READ(IN2,860) DUM
      READ(IN2,860) (ZILCH(I), I = 9,NTEMP)
C     *************
  149 READ(IN2,855) DUM
C
C*    INPUT LOG K(DISSOC) FOR DERIVED GAS SPECIES
      LGK(NS) = ZILCH(NTEMP)
      GO TO 155
C
  150 CONTINUE
      READ(IN2,860) DUM
      READ(IN2,860) DUM
      READ(IN2,855) (ZILCH(I), I = 1, 5)
C     *************
  154 CALL POWER(ZILCH,TEMPC,VAL)
      LGK(NS)=VAL
  155 CONTINUE
C
C     THE ABOVE INPUT LOG K(DISSOC) FOR DERIVED GASES (EITHER LOG K'S OR
C       REGRESSION COEFFICIENTS  LOG(K) = F(TEMPC)
C     ALL COMPLEXES WITH A NON BLANK CHARACTER IN COLUMN 16 WILL
C     BE SUPPRESSED AND DISPLAYED ON THE SCREEN
C     LGK(NS) = 566 MEANS LOG K FOR THAT TEMPERATURE, (NTEMP) IS UNKNOWN
      IF(SUP.EQ.BLNK.AND.LGK(NS).NE.566.) GO TO 116
      WRITE(IOUT2,885) NAME(NS)
      GO TO 114
  116 CONTINUE
C     SORT OUT COMPLEXES CONTAINING COMPONENTS BEING USED
C     AND REDEFINE COMPONENT INDECIES FOR DERIVED GAS SPECIES. THE
C     LOOP THAT FOLLOWS SUMS UP ALL THE GAS SPECIES COEFFICIENTS.
      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
      AZERO(NS) = COSUM - 1.
C    LGK IS THE LOG OF THE DISSOCIATION CONSTANT OF SPECIES AT TEMP
C     ITCOEF(NS)=ITOT
C*    GAMMA IS USED FOR FUGACITY COEFFICIENT IN SOLVGAS
      GAMMA(NS) = 1.0
      GO TO 113
  120 CONTINUE
C*    CALL HELHUC
C
      BSTOT = SAQTOT
      BSTP1 = SAQTOT + 1
      IF(NMTT.NE.0) CALL MINSAT(MAXSTO,MAXNM)
C     ************************
C
C     SUBROUTINE MINSAT CONTAINS READ STATEMENTS FOR MINERAL LOG K DATA
C       AND STOICHIOMETRY
      JSTOP = 0
  175 CONTINUE
      ERLOOP = 0
      GO TO 313
  176 CONTINUE
      IF(IPRINT.EQ.0) GO TO 178
      WRITE(IOUT1,999)
C     ******************************************************************
      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,905) ERLOOP, GASTOT
      WRITE(IOUT1,923)
      WRITE(IOUT1,940)
      WRITE(IOUT1,945) (NAME(I),I,MTOT(I),MTRY(I),GAMMA(I), I = 1,BSTOT)
      WRITE(IOUT1,923)
      WRITE(IOUT1,950)
      WRITE(IOUT1,955)(NAME(I),I,LGK(I),MTRY(I),GAMMA(I),I = BSTP1, NST)
      WRITE(IOUT1,923)
C     ******************************************************************
      IF(JSTOP.EQ.1) GO TO 388
  178 CONTINUE
      NDUMP = 0
      IF(FO2.NE.0.) ACT(BSTP1) = FO2
C
C     ******************************************************************
C
      DO 179 I=1,SAQTOT
  179 CONV(I) = ERPC * DABS(MTOT(I))
      IF (FO2.NE.0.) CONV(1) = ERPC
      IF (NSUNK.NE.0) CONV(NSUNK) = ERPC
C**   IF (FO2.NE.0.) CONV(1) = ERPC*MTRY(1)
C**   IF (NSUNK.NE.0) CONV(NSUNK) = ERPC*UNACT
C
C
  200 CONTINUE
      ERLOOP = ERLOOP + 1
      WRITE(IOUT2,790) ERLOOP, GASTOT
  790 FORMAT(5X,'LOOP NO. ',I4,'  GASTOT= ',E15.7)
      N = SAQTOT
C
C
C
C     ******************************************************************
C     THE DO 250 LOOP INITIALIZES THE VALUES FOR THE DERIVATIVE MATRIX
C     TO ZEROS EXCEPT FOR THE DIAGONALS WHICH ARE 1'S.
C
      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
      DO 300 KK = 1, BSTOT
C     THE DO 300 LOOP GOES THROUGH ALL COMPONENTS TO SET UP PARTIAL
C      DERIVATIVE VECTORS OF MASS BALANCE EQUATIONS
C
C
C     THE FOLLOWING IF'S CONTROL PRODUCTION OF MATRIX ENTRIES
C     FOR FO2 AND UNACT SPECIES BY APPROPRIATE CALLING OF SUBROUTINES.
      IF(KK.EQ.1.AND.FO2.NE.0.)GO TO 252
      GO TO 254
  252 CALL OXY
      GO TO 300
  254 CONTINUE
      IF(NSUNK.NE.KK) GO TO 256
      CALL OXY1
      GO TO 300
  256 CONTINUE
C
C
C
      F(KK) = -MTOT(KK) + MTRY(KK)
      DO 290 NS = BSTP1, NST
      IF(SPEC(NS,1).EQ.99) GO TO 290
C     THE DO 290 LOOP GOES THROUGH ALL DERIVED SPECIES FOR INCLUSION IN
C     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)
      DO 288 II = 2,ITOT
      J = SPEC(NS,II)
C     DTERM IS THE DERIVATIVE OF MOLALITY EXPRESSION FOR DERIVED SPECIES
C     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
      NP1 = N + 1
      DO 305 I = 1,N
      A(I,NP1) = -F(I)
  305 CONTINUE
      IF(ERLOOP.EQ.2.AND.IPRINT.NE.0) GO TO   306
      GO TO  308
  306 NSTOP = 1
      GO TO 388
  307 NSTOP = 0
      NDUMP = 1
      CALL DUMP
  308 CONTINUE
      IF(ERLOOP.EQ.1.AND.IPRINT.NE.0) CALL DUMP
      CALL SLUD(N,N)
      CALL SIR(N,N,X)
      IF(IPRINT.NE.0) CALL DUMP2
      IF(NDUMP.EQ.1) CALL DUMP2
      NDUMP = 0
      GASTOT = 0.
      DO 310 KK = 1, N
      I = KK
      COLD = MTRY(I)
      MTRY(I) = MTRY(I) + X(KK)
      XX = X(KK)
C     BRACKETS MTRYS OF COMPONENT SPECIES TO EASE CONVERGENCE
      IF(MTRY(I).GT.ALKALN) MTRY(I) = COLD
  309 CONTINUE
      IF(MTRY(I).GT.0.) GO TO 311
      XX = XX/2.
      MTRY(I) = COLD + XX
      GO TO 309
  311 IF (MTRY(KK).GT.RMIN) GO TO 310
      WRITE(IOUT2,676) NAME(KK), RMIN
  676 FORMAT(5X,'** WARNING! MOLE NUMBER OF ',A8,' LESS THAN ',D10.3)
      MTRY(KK) = RMIN
  310 CONTINUE
C
  313 CONTINUE
      BSTP2=BSTP1
      IF(FO2.EQ.0.) GO TO 314
      MTRY(BSTP1) = FO2/GAMMA(BSTP1)/(PF/GSTOLD)
      BSTP2=BSTP2 + 1
  314 DO 315 NS= BSTP2, NST
      CSAVE= -DLOG10(GAMMA(NS)) - LGK(NS)
      ITOT = IT(NS)
      DO 312 I = 2, ITOT
      J = SPEC(NS,I)
      CSAVE2 = COEF(NS,I)*(DLOG10(GAMMA(J) * MTRY(J)))
      CSAVE = CSAVE + CSAVE2
  312 CONTINUE
      PTERM = AZERO(NS)*DLOG10(PF/GSTOLD)
      RSL = CSAVE + PTERM
C     BRACKETS MTRYS OF DERIVED SPECIES
      IF(RSL.LT.RLGMN2) RSL = RLGMN2
      MTRY(NS) = 10.**RSL
      IF(MTRY(NS).GT.ALKALN) MTRY(NS) = ALKALN
  315 CONTINUE
      GASTOT = 0.
      DO 320 NS = 1,NST
      GASTOT = GASTOT + MTRY(NS)
  320 CONTINUE
      IF(ERLOOP.EQ.0) GO TO 176
C
C     ******************************************************************
C
      IF(ERLOOP.EQ.NLOOP) GO TO 410
      JK = 1
C*    JUMP = 0
      IF(FO2.EQ.0.) GO TO 352
      ACT(1) = GAMMA(1) * MTRY(1) * PF/GASTOT
      JK = 2
  352 DO 358 NS = JK, NST
C*    IF(MTRY(NS).GT.8.) JUMP = 1
C     ACT = FUGACITY
      FTERM = DLOG10(MTRY(NS))
      IF(NS.GT.BSTOT) FTERM = FTERM + AZERO(NS)*DLOG10(GSTOLD/GASTOT)
      IF(FTERM.LE.RLOGMN) GO TO 356
      ATERM = FTERM + DLOG10(GAMMA(NS) * PF / GASTOT)
      IF(ATERM.LT.RLOGMN) GO TO 356
      MTRY(NS) = 10.**FTERM
      ACT(NS) = 10.**ATERM
      GO TO 358
  356 MTRY(NS) = RMIN
      ACT(NS) = RMIN
  358 CONTINUE
C%    IF(JK.NE.2) GO TO 359
      IF(FO2.EQ.0.) GO TO 359
C     SET MTRY OF O2 IF FO2 IS INPUT, OTHERWISE, ALREADY SET, ABOVE
      MTRY(BSTP1) = FO2 * GASTOT/(GAMMA(BSTP1) * PF)
      ACT(BSTP1) = GAMMA(BSTP1) * MTRY(BSTP1) * PF/GASTOT
  359 CONTINUE
C     TESTS THAT ALL F(I) ARE LESS OR EQUAL TO ERPC
C
      IF(IPRINT.NE.0) WRITE(IOUT1,905) ERLOOP, GASTOT
C     THE FOLLOWING IS A TEST FOR CONVERGENCE.
C     TESTS THAT ALL F(I) ARE LESS OR EQUAL TO CONV(I)
      IF(ERLOOP.EQ.1) GO TO 365
      DO 362 I = 1,N
      IF(DABS(F(I)).GT.CONV(I)) GO TO 365
  362 CONTINUE
      GO TO 410
C
  365 GSTOLD = GASTOT
C     IF(JUMP.EQ.1) GO TO 367
C*    CALL GAMCAL
  367 CONTINUE
C*    IF(JUMP.EQ.1) XISAVE = XISAVE + 10.
C
C
C
      GO TO 200
  410 CONTINUE
C     SUM H2-BEARING SPECIES IF FO2 IS INPUT
      IF (FO2.EQ.0.) GO TO 413
      H2TOT = MTRY(1)
      ACT(1) = MTRY(1)*GAMMA(1)*PF/GASTOT
      DO 411 NS=BSTP1,NST
      IF(SPEC(NS,2).NE.1) GO TO 411
      H2TOT = H2TOT + COEF(NS,2)*MTRY(NS)
  411 CONTINUE
      MTOT(1) = H2TOT
  413 CONTINUE
C     SUM UNACT-BEARING SPECIES IF UNACT IS INPUT
      IF (UNACT.EQ.0.) GO TO 417
      UNTOT = MTRY(NSUNK)
      DO 416 NS=BSTP1,NST
      ITOT=IT(NS)
      DO 415 I=2,ITOT
      IF(SPEC(NS,I).NE.NSUNK) GO TO 415
      UNTOT = UNTOT + COEF(NS,I)*MTRY(NS)
  415 CONTINUE
  416 CONTINUE
      MTOT(NSUNK) = UNTOT
  417 CONTINUE
      JSTOP = 1
      IF(IPRINT.EQ.0) GO TO 388
      GO TO 176
  388 CONTINUE
C
      GRT=0.
C     COMPUTE GRT WHICH IS THE TOTAL MASS OF GAS (IN GRAMS)
      DO 406 NS=1,BSTOT
      GRT = GRT + MTOT(NS)*W(NS)
  406 CONTINUE
      WRITE(IOUT1,999)
      WRITE(IOUT1,960) TITLE1
      WRITE(IOUT1,960) TITLE2
      WRITE(IOUT1,970)
      WRITE(IOUT1,965) ERLOOP, NLOOP
      IF(TEMPC.EQ.0.) WRITE(IOUT1,961) TEMP(NTEMP)
      IF(TEMPC.NE.0.) WRITE(IOUT1,961) TEMPC
      WRITE(IOUT1,962) PF
      WRITE(IOUT1,906) GASTOT
      WRITE(IOUT1,943)
  390 WRITE(IOUT1,915)
      IF(NTEMP.NE.0) DEGC=TEMP(NTEMP)
      IF(TEMPC.NE.0.) DEGC=TEMPC
      ALA1 = DLOG10(ACT(BSTP1))
      DO 400 NS = 1, NST
      ALA = DLOG10(ACT(NS))
      IF(ALA.LT.RAT.AND.NS.NE.BSTP1) GO TO 400
C     ALM = PARTIAL PRESSURE
      ALM = ACT(NS)/GAMMA(NS)
C     ALMF = LOG OF THE MOLE FRACTION
      ALMF = ALM/PF
      ALMF = DLOG10(ALMF)
      IF(NS.GT.SAQTOT) GO TO 396
      STAR(NS)  = STA1
      IF(FO2.NE.0.) STAR(1) = STA2
      IF (NS.EQ.NSUNK) STAR (NSUNK) = STA3
      PPM = 1.E6*MTOT(NS)*W(NS)/GRT
      WRITE(IOUT1,920) NS,NAME(NS),MTRY(NS),ALM,ACT(NS), ALA, ALMF,
     1GAMMA(NS),MTOT(NS), STAR(NS), PPM
      IF(NMOLF.EQ.1) GO TO 395
      WRITE(IOUT3,921) NAME(NS),ALA,DEGC,PF,ALA1,UNALOG
      GO TO 400
  395 WRITE(IOUT3,921) NAME(NS),ALMF,DEGC,PF,ALA1,UNALOG
      GO TO 400
  396 WRITE(IOUT1,920) NS,NAME(NS),MTRY(NS),ALM,ACT(NS),ALA,ALMF,
     1GAMMA(NS)
      IF(NMOLF.EQ.1) GO TO 397
      WRITE(IOUT3,921) NAME(NS),ALA,DEGC,PF,ALA1,UNALOG
      GO TO 400
  397 WRITE(IOUT3,921) NAME(NS),ALMF,DEGC,PF,ALA1,UNALOG
  400 CONTINUE
      WRITE(IOUT1,923)
      IF(FO2.NE.0..OR.NSUNK.NE.0) WRITE(IOUT1,947)
      WRITE(IOUT1,948)
      WRITE(IOUT1,949) RAT
      IF(ERLOOP.LT.NLOOP) GO TO 405
      WRITE(IOUT1,670) NLOOP
      WRITE(IOUT2,670) NLOOP
  670 FORMAT(//5X,'***** NUMBER OF LOOPS REACHED GIVEN MAXIMUM OF ',I3,
     1 /5X,'    NON-CONVERGENCE ASSUMED, EXECUTION ABORTED'//)
      STOP
  405 IF(NSTOP.EQ.1) GO TO 307
      IF(NMTT.NE.0) CALL MINSA1
C
      GO TO 50
C
C
C
  800 FORMAT (16I5)
  820 FORMAT (A8,7X,A1,I1,7(F6.3,I2))
  822 FORMAT (A8,6X,F8.4)
  833 FORMAT (8E10.5)
  835 FORMAT (17X,E13.5)
  853 FORMAT (10X,7F10.5)
  855 FORMAT (8X,5E14.7)
  860 FORMAT (8X,8(F8.3))
  870 FORMAT (10A8)
  871 FORMAT(10A8,/,10A8,/)
  885 FORMAT(5X,A8,'  IS SUPPRESSED FROM THE CALCULATION')
  905 FORMAT (2H  ,'NO. OF LOOPS = ',I3,5X,'MOLES OF GAS =',E12.7)
  906 FORMAT (2H  , 'TOTAL MOLES OF GAS = ',E12.7)
  915 FORMAT (2H  ,'   SPECIES     MOLES      PARTIAL PRESS.    FUGACITY
     1     LOG(F)       LOG(MF)   GAMMA ',
     2'       TOT MOLES          PPM ',/)
  920 FORMAT (1X,I3,1X,A8,3X,3(E12.5,1X),2(F12.5,1X), 2(E12.5, 1X),
     1A2,1X,E13.6)
  921 FORMAT (A8,5F10.5)
  923 FORMAT (///)
  924 FORMAT (1X,I3, 1X, A8,3X,26X,E12.4,1X,F12.4,1X , 3(E12.4, 1X))
  940 FORMAT (2H  ,'SPECIES     NO.  TOT.MOLES   MTRY        GAMMA',//)
  942 FORMAT('       ')
  943 FORMAT (/)
  944 FORMAT (//)
  945 FORMAT (2H  , A8, I5, 3X, 3E12.5)
  947 FORMAT (2H  , 'CM MEANS TOTAL MOLES OF THAT SPECIES WAS COMPUTED',
     1' FROM INPUT OF ITS FUGACITY')
  948 FORMAT (2H  , 'ANY SPECIES WHOSE LOG(MOLES) = -300.00 IS ACUALLY',
     1' LESS ABUNDANT THAN THAN THAT')
  949 FORMAT (2H  , 'ANY SPECIES WITH LOG FUGACITY LESS THAN ', F6.1,
     1' IS NOT PRINTED')
  950 FORMAT (2H  ,'SPECIES     NO.    LOG K(DISS)   TRIAL MOLES  GAMMA'
     1 //)
  955 FORMAT (2H  , A8, I5, 3X, F10.4,5X,2E12.4)
  960 FORMAT (2H  , 10A8)
  961 FORMAT (/5X,'TEMPERATURE: ',F5.0,'  DEGREES CENTIGRADE')
  962 FORMAT (5X,'PRESSURE: ',F10.3,'  ATM')
  965 FORMAT (2H  ,'NO. OF LOOPS USED = ',I3,'   LOOP LIMIT = ',I3)
  970 FORMAT (2H  , 100('+'))
  999 FORMAT (1H1)
      END
      SUBROUTINE MINSAT(MAXSTO,MAXNM)
C     THIS SUBROUTINE DETERMINES WHETHER MINERALS ARE UNDER- OR SUPER-
C       SATURATED BY CALCULATING LOG(Q/K) FOR EACH MINERAL
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      DOUBLE PRECISION  KM
      CHARACTER*8 MIN,NMIN
      CHARACTER*1 SUP
      INTEGER SAQ, QAS, SAQTOT
      INTEGER SPECM, SPMM
      DIMENSION MIN(900), COEFM(900,8), SPECM(900,8),ZILCH(13),KM(900)
     1 , NMSAV(900), XLGQM(900), ISOLS(900), B(900),IT(900),XMOLEW(900),
     2 TEMP(13)
      COMMON/MINNIE/  ACT(999), TEMPC,NTEMP,NMTT
      COMMON/TURNIP/ SAQ(50), QAS(50), SAQTOT
      COMMON/IO/IPRINT,IN1,IN2,IN3,IOUT1,IOUT2,IOUT3,IOUT4,ITREF
      COMMON/EREBUS/ PF,GASTOT,UNALOG,ALA1
      DATA TEMP/25.,100.,200.,300.,400.,500.,600.,700.,800.,900.,1000.,
     11100.,1200./
      IF(NTEMP.NE.0) DEGC=TEMP(NTEMP)
      IF(TEMPC.NE.0.) DEGC=TEMPC
      NM = 0
   33 CONTINUE
      NM = NM + 1
      IF(NM.LE.MAXNM) GO TO 30
      NM = NM - 1
      WRITE(IOUT1,664) MAXNM, MIN(NM)
      WRITE(IOUT2,664) MAXNM, MIN(NM)
  664 FORMAT(/5X,'NUMBER OF SOLIDS AND LIQUIDS READ IN GASTHERM EXCEEDS'
     1 ,' MAXIMUM OF ',I3,/5X,'LAST SPECIES READ IS:  ',A8//)
      STOP
   30 READ(IN2,854) NMIN,SUP,XMOLEW(NM),ITOT,(COEFM(NM,I),SPECM(NM,I)
     1 ,I=1,ITOT),B(NM), ISOLS(NM)
      MIN(NM) = NMIN
C     **********************
C
C     COEFM AND SPECM CONTAIN THE STOICHIOMETRIC COEFFICIENTS AND
C       SPECIES FOR THE CHARACTERISTIC REACTION THAT RELATES THE MINERAL
C       TO THE COMPONENT SPECIES
C     B IS USED FOR SOLID SOLUTION CALCULATIONS AND IS GENERALLY 1.
C     ISOLS IDENTIFIES MINERALS AS ENDMEMBERS IN PARTICULAR SOLID
C       SOLUTIONS.  MINERALS WITH THE SAME ISOLS ARE IN THE SAME SOLID
C       SOLUTION.
C     ALL OF THIS IS READ FROM FILE GASTHERM
      IF(ITOT.EQ.0) GO TO 165
C     CHECK COMPATIBILITY WITH ARRAY SIZE
      IF(ITOT.LE.MAXSTO) GO TO 25
      WRITE(IOUT1,666) MIN(NM), ITOT, MAXSTO
      WRITE(IOUT2,666) MIN(NM), ITOT, MAXSTO
  666 FORMAT(/5X,'SPECIES ',A8,' IN GASTHERM HAS TOO MANY COMPONENTS',
     1 /5X,'ITOT FOR THAT SOLID OR LIQUID IS ',I3,' ; MAXIMUM IS ',I3//)
      STOP
   25 CONTINUE
      IF(TEMPC.NE.0.) GO TO 150
      IF(NTEMP.GT.8) GO TO 146
      READ(IN2,865) (ZILCH(I), I = 1,NTEMP)
      READ(IN2,865) DUM
      GO TO 147
  146 READ(IN2,865) (ZILCH(I), I = 1,8)
      READ(IN2,865) (ZILCH(I), I = 9,NTEMP)
C     *******************
  147 READ(IN2,855) DUM
C
C     ZILCH CONTAINS LOG K FOR MINERAL REACTIONS
      VAL=ZILCH(NTEMP)
      KM(NM) = VAL
      GO TO 155
C
  150 CONTINUE
      READ(IN2,865) DUM
      READ(IN2,865) DUM
      READ(IN2,855) (ZILCH(I), I = 1, 5)
C     *******************
  154 CALL POWER(ZILCH,TEMPC,VAL)
      KM(NM) = VAL
C
  155 CONTINUE
      DO 40 I = 1, ITOT
      JS = SPECM(NM,I)
      IF(QAS(JS).EQ.0) GO TO 30
      SPECM(NM,I) = QAS(JS)
   40 CONTINUE
      IT(NM) = ITOT
      IF(IPRINT.NE.0) WRITE(IOUT1,977) MIN(NM), VAL
      GO TO 33
  165 CONTINUE
      NMTT = NM - 1
      RETURN
      ENTRY MINSA1
      MCNT = 0
      DO 101 NM = 1, NMTT
      CSAVE = 0.
      ITOT = IT(NM)
      DO 90 I = 1,ITOT
      SPMM = SPECM(NM,I)
C     COMPUTE ACTIVITY QUOTIENT Q FOR MINERAL
      CSAVE1=COEFM(NM,I)*DLOG10(ACT(SPMM))
      CSAVE= CSAVE + CSAVE1
   90 CONTINUE
      XLGQM(NM) = CSAVE
   96 IF(CSAVE.GT.KM(NM)) GO TO 98
      GO TO 101
   98 MCNT = MCNT + 1
      NMSAV(MCNT) = NM
  101 CONTINUE
C     WRITE SUBPLOT FILE TO PLOT LOG (Q/K)
      DO 108 NM = 1, NMTT
      RATLG = XLGQM(NM) - KM(NM)
      WRITE(IOUT4,990) MIN(NM),RATLG,DEGC,PF,ALA1,UNALOG
  108 CONTINUE
      WRITE(IOUT1,923)
      WRITE(IOUT1,970)
      WRITE(IOUT1,923)
      WRITE(IOUT1,930)
      WRITE(IOUT1,923)
      IF(MCNT.EQ.0) GO TO 110
      WRITE(IOUT1,975)
      WRITE(IOUT1,980)
C     WRITE LOG (Q/K) DATA ON GASOUT FILE
      DO 109 I = 1, MCNT
      NM = NMSAV(I)
      RATLG = XLGQM(NM) - KM(NM)
      WRITE(IOUT1,977) MIN(NM), KM(NM), XLGQM(NM), RATLG
      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
      RATLG =  XLGQM(NM) - KM(NM)
      WRITE(IOUT1,977) MIN(NM), KM(NM), XLGQM(NM), RATLG
  112 CONTINUE
      WRITE(IOUT1,923)
  107 CONTINUE
      RETURN
  854 FORMAT(A8,A1,1X,F7.2,I2,1X, 7(F6.3, I2),F4.3,/,I2,F6.3,I2,F6.3)
  855 FORMAT(8X,5E14.7)
  865 FORMAT  (8X,8F8.3)
  923 FORMAT (///)
  930 FORMAT (2H  ,' NOTE: FOR GASES LISTED BELOW, ',
     1'LOG (FUGACITY OF GAS =  LOG (Q/K)')
  970 FORMAT (2H  , 100('+'))
  975 FORMAT (2H  , 'THE GAS MIXTURE IS SUPERSATURATED WITH THE',
     1' FOLLOWING CRYSTALLINE SOLIDS OR LIQUIDS', //)
  977 FORMAT(2H  ,A8,2X,3(2X, F7.2))
  978 FORMAT(2H  ,A8,2X,2(2X, E7.2))
  979 FORMAT (2H  , ' THIS GAS MIXTURE IS UNDERSATURATED WITH THE',
     1' FOLLOWING CRYSTALLINE SOLIDS OR LIQUIDS',//)
  980 FORMAT        ('  MINERAL       LOG K    LOG Q   LOG Q/K',/)
  985 FORMAT (2H  ,'SOLID SOLUTIONS',//,'  MINERAL        X      SUM',/)
  990 FORMAT(A8,5F10.5)
      END
      SUBROUTINE OXY
      IMPLICIT  DOUBLE PRECISION  (A-H,O-Z)
      DOUBLE PRECISION  MTRY,K,LGK
      INTEGER BSTP1,BSTOT,SPEC
      COMMON/SLOUGH/  A(50,51), T(50,51),IV(50),N
      COMMON/MORT/ MTRY(999),GAMMA(999),FO2,F(50),KK
      COMMON/DRUNK/COEF(999,7),K(999),LGK(999),SPEC(999,7),IT(999),
     1 BSTP1,BSTOT
      COMMON/EREBUS/ PF,GASTOT,UNALOG,ALA1
      COMMON/HUH/ UNACT, TERMLN,NSUNK
      COMMON/IO/IPRINT,IN1,IN2,IN3,IOUT1,IOUT2,IOUT3,IOUT4,ITREF
      IF(IPRINT.NE.0) WRITE(IOUT1,911)
C
C
C
C     PSEUDO SUBROUTINE FH2CALC
C     O2 EQUATION IN LOG FORM.  NO PF/GASTOT TERM
      F(1)=-0.5*LGK(BSTP1)+DLOG10(GAMMA(2)*MTRY(2))
     1 -DLOG10(GAMMA(1)*MTRY(1))-0.5*DLOG10(FO2)
C     DERIVATIVES WITH RESPECT TO H2(1) AND H2O(2)
      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
      RETURN
  911 FORMAT (2H  , 40X, 'SUBROUTINE OXY WAS USED')
C
C
      ENTRY OXY1
      F(NSUNK) = -MTRY(NSUNK)/GASTOT*PF + UNACT/GAMMA(NSUNK)
      DO 67 J = 1, N
   67 A(NSUNK,J) = 0.0
      A(NSUNK,NSUNK) = -PF/GASTOT
      RETURN
      END
      SUBROUTINE SLUD (ADIM,TDIM)
      IMPLICIT DOUBLE PRECISION (A-G,O-Z)
      INTEGER ADIM,TDIM
      COMMON/SLOUGH/ A(50,51),T(50,51),IV(50),N
C
      IF(ADIM.LT.N) GO TO 8110
      DO 8100 J=1,N
      DO 8100 I=1,N
      T(I,J) = A(I,J)
 8100 CONTINUE
 8110 CONTINUE
C
      IV(N) = 1
      DO 8260 K=1,N
      PIV = DABS(T(K,K))
      IF(K.GE.N) GO TO 8260
C
      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
C
      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
C
      TMP = -TMP
      DO 8230 I=KP1,N
      T(I,K) = T(I,K)/TMP
 8230 CONTINUE
C
      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 (ADIM,TDIM,X)
      IMPLICIT DOUBLE PRECISION (A-G,O-Z)
      INTEGER ADIM,TDIM
      DIMENSION X(50),B(50),R(50)
      COMMON/IO/IPRINT,IN1,IN2,IN3,IOUT1,IOUT2,IOUT3,IOUT4,ITREF
      COMMON/SLOUGH/ A(50,51),T(50,51),IV(50),N
C
C     B, THE VECTOR CONTAINING THE RIGHT-HAND SIDE OF THE
C     SYSTEM OF LINEAR EQUATIONS, IS CONTAINED IN THE LAST
C     COLUMN OF THE MATRIX A(N,N+1)
C
      DO 150 I=1,N
      B(I)=A(I,N+1)
  150 CONTINUE
C
      IER = 0
      DO 8100 I=1,N
      X(I) = B(I)
 8100 CONTINUE
      IF(IV(N).NE.0)  GO TO 8500
C
      WRITE(IOUT1,915)
  915 FORMAT(/,5X,'************** WARNING **************',/,2X,
     1'COMPUTATIONAL SINGULARITY, IER = 0.0',/)
      STOP
C
 8500 CONTINUE
      QUOT = 1.0
C
      CALL SBS (TDIM,X)
      IF(ITREF.EQ.0) RETURN
 8200 CALL SAXMB (ADIM,X,B,R)
      IF(ITREF.EQ.1) RETURN
      CALL SBS  (TDIM,R)
      IER = IER+1
C
      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
C
      IF(XN1.GT.0.0) RN1 = RN1/XN1
      IF((RN1+RN1).LE.QUOT) GO TO 8220
      IER = -IER
      RETURN
 8220 CONTINUE
      QUOT = RN1
      IF((1.0 + SNGL(RN1)).NE.1.0) GO TO 8200
      RETURN
      ENTRY DUMP
      NP1 = N + 1
      WRITE(IOUT1,923)
      WRITE(IOUT1,974) N,NP1
      DO 507 I = 1, N
      WRITE(IOUT1,973) (A(I,J), J = 1, NP1)
  507 CONTINUE
      RETURN
      ENTRY DUMP2
      WRITE(IOUT1,973) (X(I), I = 1, N)
      WRITE(IOUT1,910) IER
      IF(ITREF.NE.0) WRITE(IOUT1,973) (R(I), I=1,N)
      RETURN
  910 FORMAT(8X,'NUMBER OF ITERATIVE REFINEMENTS : = ',I4)
  923 FORMAT (///)
  973 FORMAT (2H  , 10E11.4)
  974 FORMAT (2H  , 'MATRIX OF A(I,J), I=1,',I2,' J=1,',I2,//)
  922 FORMAT (//)
      END
      SUBROUTINE SBS (TDIM,B)
      IMPLICIT DOUBLE PRECISION (A-G,O-Z)
      INTEGER TDIM
      COMMON/SLOUGH/ A(50,51),T(50,51),IV(50),N
C
      DIMENSION B(50)
      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
C
      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 (ADIM,X,B,R)
      IMPLICIT DOUBLE PRECISION (A-G,O-Z)
      INTEGER ADIM
      DIMENSION X(50),B(50),R(50)
      COMMON/SLOUGH/ A(50,51),T(50,51),IV(50),N
C
      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 POWER(COEF,TEMP,VAL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION COEF(8)
      COMMON/IO/IPRINT,IN1,IN2,IN3,IOUT1,IOUT2,IOUT3,IOUT4,ITREF
      TEMPK=TEMP + 273.15
      VAL = COEF(1)
      VAL = VAL + COEF(2)/TEMPK
      VAL = VAL + COEF(3)*TEMPK
      VAL = VAL + COEF(4)/TEMPK**2
      VAL = VAL + COEF(5)*DLOG10(TEMPK)
      RETURN
      END