C**********************************************************************
C
C   PROGRAM GASCAL IS DESIGNED TO COMPUTE STANDARD STATE APPARANT
C   ENERGIES FOR  MINERALS AND  GASES AND WATER  AT VARIOUS TEMPERATURES
C   AND PRESSURES.  THESE VALUES ARE USED TO COMPUTE THE LOG OF THE
C   LOG OF THE EQUILIBRIUM CONSTANT FOR ANY REACTION  INVOLVING
C   GASES AND MINERALS.
C        THIS PROGRAM ALLOWS THE USER TO INTERACTIVELY TYPE IN
C   STOCHIOMETRIES OF REACTIONS. IT CAN ALSO READ THE STOCHIOMETRIES
C   OF GASES OR MINERALS  DIRECTLY FROM THE GASTHERM DATA FILE.
C   WHEN DOING SO, THE USER CAN CHOSE TO CREATE A NEW GASTHERM DATA
C   FILE CONTAINING THE NEW LOG K VALUES. THE DIFFERENCE BETWEEN THE
C   EXISTING VALUES AND THE NEW VALUES IS ALSO CALCULATED.
C        GASCAL WAS WRITTEN AFTER PROGRAM KCAL, WHICH WAS WRITTEN BY
C   NIC SPYCHER AT THE UNIVERSITY OF OREGON
C        THIS PROGRAM NEEDS THE FOLLOWING INPUT FILES TO RUN CORRECTLY
C     - A MINERAL DATA FILE
C     - A GAS DATA FILE
C     - A P-T GRID CONTAINING THE PRESSURES AND TEMPERATURES OF THE RUN
C     - (OPTIONAL) GASTHERM DATA FILE CONTAINING STOCHIOMETRIES OF
C       DISSOCIATION OF GASES AND MINERALS.
C
C             WRITTEN BY R. SYMONDS, MICHINGAN TECH, FEBRUARY  1989
C
C**********************************************************************
C**********************************************************************
C
C           LAST REVISED :  FEF 89         BY: R. SYMONDS
C
C**********************************************************************
C
      DOUBLE PRECISION RENAM1,RENAM2,ION,GAS,MIN,MINAM,BLANK,SPECM,NMIN
      INTEGER SPEC,SPECG
      COMMON/DATA/TREF,PREF,T(24),P(24,10),TC(24)
      COMMON/MAX/JMAX,KMAX,INDP(10),INDT(24),ITOTT,ITOTP,MAXG
      COMMON/GIBBS/G(40,21,10),GG(999,21),GM(8,21,10)
      COMMON/DATMIN/A,B,C,VREF,HRF,SRF,GRF,A5,B5,C5,A6,B6,C6,TTRAN1,
     1TTRAN2,VRE1,VRE2,HTRAN1,HTRAN2,SLOPE1,SLOPE2,TTRAN3
      COMMON/DATGAS/SGREF(999),GGREF(999),AG(999),BG(999),CG(999),
     1 DG(999),HGREF(999)
      COMMON/STOCH/ ACTI(10),ACTM(8),XMOLG(8),XRI(10),XRG(8),XRM(8),
     1 XRW,ACTW,NSION,NSGAS,NSMIN,SPEC(10),SPECG(8)
      COMMON/NOMS/GAS(999,2),ION(40),MINAM(8)
      COMMON/IO/INP1,INP2,INP3,INP4,INP5,IOUT1,IOUT2,IOUT3,ITTY
C
      DIMENSION RENAM1(10),RENAM2(10),MIN(2),SPECM(8)
      DATA BLANK/'        '/, POS/'Y'/
C
C     FORTRAN I/O VARIABLES:
C       INP1: NOT USED NOW
C       INP2: INPUT FILE W/ P-T AND DATA FOR WATER
C       INP3: INPUT FILE W/ GASES THERMO DATA
C       INP4: INPUT FILE W/ MINERALS THERMO DATA
C       INP5: GASTHERM DATA FILE  (ORIGINAL)
C       IOUT1: MAIN OUTPUT
C       IOUT2: GASTHERM FORMATTED OUTPUT WITH NEW LOG K'S
C       IOUT3: SCRATCH FILE WITH DIFFERENCES BETWEEN NEW AND OLD LOG K'S
C       ITTY:  TTY INPUT OR OUTPUT
C
C     THESE VARIALBLES CORRESPOND TO THE FOLLOWING NUMBERS
      INP1=10
      INP2=11
      INP3=12
      IOUT1=6
      IOUT2=15
      IOUT3=16
      ITTY=5
      INP4=13
      INP5=14
C
C     OPEN STATEMENTS FOR IBM DOS MACHINES
C
      OPEN (UNIT = 11, FILE = 'PT.DAT',  STATUS = 'OLD')
      OPEN (UNIT = 12, FILE = 'GAS.DAT',   STATUS = 'OLD')
      OPEN (UNIT = 13, FILE = 'MINERAL.DAT',   STATUS = 'OLD')
      OPEN (UNIT = 14, FILE = 'GASTHERM.DAT', STATUS = 'OLD')
      OPEN (UNIT = 6, FILE = 'GCALOUT.DAT',   STATUS = 'UNKNOWN')
      OPEN (UNIT = 15, FILE = 'GASTHERM.NEW', STATUS = 'UNKNOWN')
      OPEN (UNIT = 16, FILE = 'KDIFF.DAT',  STATUS = 'UNKNOWN')
C
      TREF = 298.15
      PREF = 1.00000
      MAXG=999
      NSGAS = 0
      NSMIN = 0
      DO 351  I=1,MAXG
      DO 350  II=1,21
        GG(I,II) = 0.
  350 CONTINUE
  351 CONTINUE
      DO 363 I = 1,40
      DO 362 II = 1,21
      DO 361 III = 1,10
        G(I,II,III) = 0.
        GM(I,II,III) = 0.
  361 CONTINUE
  362 CONTINUE
  363 CONTINUE
C
  501 FORMAT(1X,F11.5,2X,F5.0,1X,5(E12.6))
  502 FORMAT(16I5)
  503 FORMAT (20A4)
  505 FORMAT(1X,3A4//)
  506 FORMAT(F10.3,3G15.5)
  508 FORMAT(A8,2X,3F10.3,I5)
  509 FORMAT(2(F10.3,F15.3)/4F10.3)
  511 FORMAT(1X,/)
  514 FORMAT(1X,F11.5,5X,F5.0)
  515 FORMAT(4X,'P (BARS)',5X,'T (DEGREES C)',/)
  521 FORMAT(A8,A6,2F9.1,5F8.4)
  526 FORMAT(2A8,44X,A8)
  530 FORMAT(1X,7HMINERAL,3X,F10.2)
  531 FORMAT(I5,1X,7(I4,F6.0))
  532 FORMAT(8(I5,F5.0))
  533 FORMAT(1X,2A4,2X,F10.2)
  536 FORMAT(1X,5HH(2)0,5X,F10.2)
  537 FORMAT(2F5.2,2I5)
  551 FORMAT(1H1)
  560 FORMAT(//6X,'MINERAL',9X,'GREF',10X,'HREF',7X,'SREF',6X,'VREF',
     1 13X,'A',7X,'B*10**3',3X,'C*10**-5'/,6X,'-------',9X,'----',10X,
     2 '----',7X,'----',6X,'----',13X,'-',7X,'------',3X,'--------')
C
C
C
C     READS DATA FOR P AND T
C     ---------------------------------------
      READ(INP2,511)
      READ(INP2,502) JMAX,KMAX
      ITOTP=KMAX
      ITOTT=JMAX
      DO 13 K=1,KMAX
      INDP(K)=K
      DO 15 J=1,JMAX
   15 READ(INP2,501) P(J,K),TC(J)
   13 CONTINUE
      DO 18 J=1,JMAX
      INDT(J)=J
      T(J)=TC(J)+273.15
   18 CONTINUE
C
      WRITE(ITTY,710)
  710 FORMAT(//5X,'THE P-T DATA HAS BEEN READ....'/
     1 ,7X,'DO YOU WANT TO ECHO IT ON THE SCREEN ?'//)
      READ(ITTY,602) ANS
      IF(ANS.NE.POS) GO TO 19
      WRITE(ITTY,515)
      DO 17 K=1,KMAX
      DO 16 J=1,JMAX
   16 WRITE (ITTY,514) P(J,K),TC(J)
   17 CONTINUE
   19 CONTINUE
C
C
C     READS DATA FOR GASES
C     --------------------------------
      READ(INP3,511)
      DO 20 M=1,MAXG
   20 READ(INP3,521,END=24)(GAS(M,M1),M1=1,2),GGREF(M),HGREF(M),
     1 SGREF(M),AG(M),BG(M),CG(M),DG(M)
   24 IF(M.LT.MAXG) MAXG = M - 1
      WRITE(ITTY,715)
  715 FORMAT(//5X,'THE DATA FOR THE GASES HAS BEEN READ...'/)
C
   22 CONTINUE
C
C     THE FOLLOWING DISPLAYS THE T-P GRID AND GIVE THE USER THE CHOICE
C     OF P-T TO BE USED IN THE CALCULATION
C     ----------------------------------------------------------------
C
C       ITOTP: TOTAL NUMBER OF SELECTED PRESSURES
C       INDP(K) K=1,TOTP : ARRAY OF SELECTED PRESSURES
C       ITOTT: TOTAL NUMBER OF SELECTED TEMPERATURES
C       INDT(J) J=1,ITOTT :ARRAY OF SELECTED TEMPERATURES
C
      WRITE(ITTY,760)
  760 FORMAT(///5X,'DO YOU WANT TO CARRY CALCULATIONS ON THE WHOLE'
     1 ,' P-T GRID ?'/)
      READ(ITTY,602) ANS
      IF(ANS.EQ.POS) GO TO 50
C      IF SATURATION P-T CURVE
      IF(P(1,1).NE.1.) GO TO 51
      WRITE(ITTY,763)
  763 FORMAT(//5X,'WATER SATURATION P-T CURVE'/)
      GO TO 57
C
   51 WRITE(ITTY,761)
  761 FORMAT(4X,'INDEX   P BARS'/)
   70 DO 52 K=1,KMAX
      INDP(K)=0
   52 WRITE(ITTY,762) K, P(1,K)
  762 FORMAT(5X,3(I2,':',1X,F10.0,5X))
      WRITE(ITTY,764)
  764 FORMAT(/5X,'SELECT DESIRED PRESSURE INDEXES (END BY SLASH)'/)
      READ(ITTY,*) (INDP(K), K=1,KMAX)
C
   57 WRITE(ITTY,765)
  765 FORMAT(4X,3('INDEX  T DEG.C',5X),/)
      DO 54 I=1,11
      I2=I+11
      INDT(I)=0
      INDT(I2)=0
   54 WRITE(ITTY,762) I,TC(I),I2,TC(I2)
      WRITE(ITTY,766)
  766 FORMAT(/5X,'SELECT DESIRED TEMPERATURE INDEXES (END BY SLASH)'/)
      READ(ITTY,*) (INDT(I), I=1,JMAX)
C
      DO 56 K=1,KMAX
      IF(INDP(K).EQ.0) GO TO 58
   56 CONTINUE
   58 ITOTP=K-1
      DO 60 I=1,JMAX
      IF(INDT(I).EQ.0) GO TO 62
   60 CONTINUE
   62 ITOTT=I-1
      WRITE(ITTY,770)
  770 FORMAT(5X,'THE FOLLOWING PRESSURES AND TEMPERATURES HAVE BEEN',
     1 ' CHOSEN'/5X,'  IS THIS CORRECT ?'/)
      IF(P(1,1).NE.1.) GO TO 64
      WRITE(ITTY,763)
      GO TO 66
   64 WRITE(ITTY,772) (P(1,INDP(K)), K=1,ITOTP)
   66 WRITE(ITTY,772) (TC(INDT(I)), I=1,ITOTT)
  772 FORMAT(5X,6(F10.0,2X)/)
      READ(ITTY,602) ANS
      IF(ANS.NE.POS) GO TO 70
   50 CONTINUE
C
      WRITE(ITTY,774)
  774 FORMAT(//5X,'DO YOU WANT TO READ STOCHIOMETRIES FROM GASTHERM ?'/)
      READ(ITTY,602) ANS
      IF(ANS.NE.POS) GO TO 777
      CALL READER
      STOP
C
C     THE FOLLOWING READS STOCHIOMETRY FOR A GIVEN REACTION INVOLVING
C     MINERALS AND/OR GASES AND/OR IONS INTERACTIVELY
C-----------------------------------------------------------------------
  777 WRITE(ITTY,730)
  730 FORMAT(//5X,'INPUT REACTION NAME AND/OR TITLE (2 LINES)'/)
      READ(ITTY,600) (RENAM1(I), I=1,10)
      READ(ITTY,600) (RENAM2(I), I=1,10)
      WRITE(IOUT1,551)
      WRITE(IOUT1,601) (RENAM1(I),I=1,10)
      WRITE(IOUT1,601) (RENAM2(I),I=1,10)
      WRITE(IOUT1,603)
C
   95 WRITE(ITTY,732)
  732 FORMAT(/5X,'INPUT TOTAL NUMBER OF GASES AND MINERALS IN',
     1 ' THE REACTION'/)
      READ(ITTY,*) NSGAS,NSMIN
      WRITE(ITTY,731) NSGAS,NSMIN
  731 FORMAT(//5X,'THERE ARE',I2,' GASES AND ',I2,
     1 ' MINERALS'/5X,'INVOLVED IN THIS REACTION. IS THIS CORRECT ?'/)
      READ(ITTY,602) ANS
      IF(ANS.NE.POS) GO TO 95
C
C     READS STOCHIOMETRY FOR GASES, THEN CALC. FREE ENERGY
C    -----------------------------
  105 IF(NSGAS.EQ.0) GO TO 110
C*
      WRITE(ITTY,737)
  737 FORMAT(//5X,'DO YOU WANT THE GAS SPECIES PRINTED ON THE SCREEN?'/)
      READ(ITTY,602) ANS
      IF(ANS.NE.POS) GO TO 109
      DO 111 I=1,MAXG
  111 WRITE(ITTY,733) I,GAS(I,1)
  733 FORMAT(5X,4(I3,1X,A8,6X))
C
  109 WRITE(ITTY,738)
  738 FORMAT(//5X,'INPUT STOCH. COEFFICIENTS, INDEXES AND MOLE'
     1, ' FRACTION'/,5X,'FOR THE GASES (1 LINE)'/)
      READ(ITTY,*) (XRG(I),SPECG(I),XMOLG(I),I=1,NSGAS)
      DO 112 I=1,NSGAS
      IF(XMOLG(I).EQ.0.) XMOLG(I)=1.
  112 WRITE(ITTY,615) XRG(I),GAS(SPECG(I),1),XMOLG(I)
      WRITE(ITTY,736)
  736 FORMAT(/,5X,'IS THIS CORRECT?',/)
      READ(ITTY,602) ANS
      IF(ANS.NE.POS) GO TO 105
      CALL WRITGA
      CALL GASES
C
C     READS STOCHIOMETRY FOR MINERALS
C   ---------------------------------
  110 IF(NSMIN.EQ.0) GO TO 170
      DO 115 I=1,NSMIN
      WRITE(ITTY,739)
  739 FORMAT(//5X,'INPUT MINERAL NAME AS IT APPEARS IN DATA FILE'/)
      READ(ITTY,600) MINAM(I)
      WRITE(ITTY,741) MINAM(I)
  741 FORMAT(8X,'INPUT COEFFICIENT AND ACTIVITY OF ',A8/)
      READ(ITTY,*) XRM(I),ACTM(I)
  115 CONTINUE
      WRITE(ITTY,740)
  740 FORMAT(//5X,'DO YOU WANT TO CHANGE DATA FOR THE MINERALS ?'/)
      READ(ITTY,602) ANS
      IF(ANS.EQ.POS) GO TO 110
C
C     READS ONLY THE DATA FOR MINERALS THAT ARE  SPECIFIED
C     THEN CALCULATES FREE ENERGY
C
      IFLAG=0
      READ(INP4,511)
      DO 41 I=1,8
   41 SPECM(I)=BLANK
      WRITE(IOUT1,560)
   40 CONTINUE
      IF(IFLAG.GE.NSMIN) GO TO 125
      READ(INP4,526,END=100)(MIN(LN),LN=1,2),NMIN
      DO 42 I=1,NSMIN
      IF(MINAM(I).EQ.MIN(1)) GO TO 44
   42 CONTINUE
      READ(INP4,620,END=100)
      GO TO 40
   44 CONTINUE
      IND=I
      SPECM(I)=NMIN
      IFLAG=IFLAG+1
      WRITE(ITTY,744) MINAM(I),XRM(I),ACTM(I)
  744 FORMAT(5X,'FOUND : ',A8,' COEF= ',F6.3,'  ACT= ',F4.2)
      CALL README(MIN)
C
      CALL SOLIDS(IND)
C
      GO TO 40
C
  100 CONTINUE
C
      DO 120 I=1,NSMIN
      IF(SPECM(I).EQ.BLANK) WRITE(ITTY,748) MINAM(I)
  748 FORMAT(/5X,'I CANNOT FIND THE FOLLOWING MINERAL : ',A8)
  120 CONTINUE
      REWIND INP4
      WRITE(IOUT1,749)
  749 FORMAT(/10X,'***** CORRECTION,  NEW GIVEN MINERALS BELOW'/)
      GO TO 110
C
  125 CONTINUE
      WRITE(ITTY,746)
  746 FORMAT(/5X,11HTHAT'S IT !,'IS THIS CORRECT ?'/)
      READ(ITTY,602) ANS
      IF (ANS.EQ.POS) GO TO 170
      WRITE(IOUT1,749)
      GO TO 110
C
  170 CALL REAC
C
      WRITE(ITTY,756)
  756 FORMAT(//5X,'DO YOU WANT TO CALCULATE ANOTHER REACTION ?'/)
      READ(ITTY,602,END=416) ANS
  416 IF(ANS.NE.POS) STOP
      REWIND INP4
      WRITE(ITTY,768)
  768 FORMAT(5X,'DO YOU WANT TO KEEP THE SAME P-T GRID ?'/)
      READ(ITTY,602) ANS
      IF(ANS.EQ.POS) GO TO 50
      DO 460 K=1,KMAX
  460 INDP(K)=K
      DO 465 J=1,JMAX
  465 INDT(J)=J
      ITOTP=KMAX
      ITOTT=JMAX
      GO TO 22
C
  600 FORMAT(10A8)
  601 FORMAT(2X,10A8)
  602 FORMAT(A1)
  603 FORMAT(1X,120('+'))
  605 FORMAT(4I5)
  615 FORMAT(5X,F6.3,1X,':',A8,2X,F4.2)
  620 FORMAT(///)
C
      END
C
C
      SUBROUTINE GASES
      DOUBLE PRECISION GAS,MINAM,ION
      INTEGER SPECG,SPEC
      COMMON/GIBBS/G(40,21,10),GG(999,21),GM(8,21,10)
      COMMON/DATGAS/SGREF(999),GGREF(999),AG(999),BG(999),CG(999),
     1 DG(999),HGREF(999)
      COMMON/STOCH/ ACTI(10),ACTM(8),XMOLG(8),XRI(10),XRG(8),XRM(8),
     1 XRW,ACTW,NSION,NSGAS,NSMIN,SPEC(10),SPECG(8)
      COMMON/NOMS/GAS(999,2),ION(40),MINAM(8)
      COMMON/IO/INP1,INP2,INP3,INP4,INP5,IOUT1,IOUT2,IOUT3,ITTY
      COMMON/MAX/JMAX,KMAX,INDP(10),INDT(24),ITOTT,ITOTP,MAXG
      COMMON/DATA/TREF,PREF,T(24),P(24,10),TC(24)
C
C**********************************************************************C
C                                                                      C
C        GASES  CALCULATES STANDARD STATE THERMODYNAMIC PROPERTIES     C
C   FOR GASES AT GIVEN TEMPERATURE AND PRESSURE.                       C
C   THE STANDARD STATES WERE TAKEN AT T AND ONE BAR FOR THE IDEAL GAS. C
C                                                                      C
C       TREF  =  REFERENCE TEMPERATURE (298.15 K)                      C
C       PREF  =  REFERENCE PRESSURE (ONE ATMOSPHERE)                   C
C       T(J)  =  INDEX ON TEMPERATURE                                  C
C                                                                      C
C       AG  =                                                          C
C       BG  =     MAIER-KELLY HEAT CAPACITY COEFFICIENTS FOR GASES     C
C       CG  =        DG FROM BARIN AND KNACK (1973)                    C
C       DG  =                                                          C
C                                                                      C
C       GG  =  STANDARD STATE APPARENT FREE ENERGY FOR A GAS AT        C
C              T AND P                                                 C
C       SG  =  THIRD LAW ENTROPY FOR A GAS AT T AND P                  C
C       M   =  INDEX ON THE GAS                                        C
C                                                                      C
C       SGREF  =  STANDARD STATE THIRD LAW ENTROPY AT TREF AND PREF    C
C       GGREF  =  STANDARD STATE FREE ENERGY OF FORMATION AT TREF      C
C                 AND PREF                                             C
C                                                                      C
C**********************************************************************C
C     K=INDP(NK)
      WRITE(IOUT1,104)
  104 FORMAT(/5X,'FREE ENERGIES FOR GASES ARE CALCULATED AT ',
     1 'ONE BAR (STANDARD STATE)')
      WRITE(IOUT1,106) (GAS(SPECG(M),1), M=1,NSGAS)
  106 FORMAT(5X,'APPARENT FREE ENERGIES OF GASES AT GIVEN T:'/
     1 ,5X,42('-'),/3X,'TEMP DEG.C',5X,A8,7(6X,A8))
      DO 10 NJ=1,ITOTT
      J=INDT(NJ)
      DO 12 I = 1,NSGAS
      M=SPECG(I)
      BB=BG(M)*.001
      CC=CG(M)*100000.
      DD=DG(M)/1000000.
C
      G1=0.
      G0=GGREF(M)-SGREF(M)*(T(J)-TREF)
      IF(T(J).NE.TREF) G1=GTG(T(J),TREF,T(J),AG(M),BB,CC,DD)
      GG(M,J)= G0+G1
   12 CONTINUE
      WRITE(IOUT1,110) TC(J),(GG(SPECG(I),J), I=1,NSGAS)
  110 FORMAT(5X,F5.0,5X,8(F12.0,2X))
 10   CONTINUE
      RETURN
      END
      SUBROUTINE SOLIDS(IND)
      COMMON/DATA/TREF,PREF,T(24),P(24,10),TC(24)
      COMMON/MAX/JMAX,KMAX,INDP(10),INDT(24),ITOTT,ITOTP,MAXG
      COMMON/GIBBS/G(40,21,10),GG(999,21),GM(8,21,10)
      COMMON/DATMIN/A,B,C,VREF,HRF,SRF,GRF,A5,B5,C5,A6,B6,C6,TTRAN1,
     1TTRAN2,VRE1,VRE2,HTRAN1,HTRAN2,SLOPE1,SLOPE2,TTRAN3
      COMMON/IO/INP1,INP2,INP3,INP4,INP5,IOUT1,IOUT2,IOUT3,ITTY
      DIMENSION GI(5)
C
C**********************************************************************C
C                                                                      C
C        SOLIDS  CALCULATES THERMODYNAMIC DATA AT P AND T FOR MINERAL  C
C   PHASES USING MAIER-KELLY HEAT CAPACITY COEFFICIENTS AND CONSTANT   C
C   VOLUME (T,PREF).  THE STANDARD STATES WERE TAKEN AT T AND P FOR    C
C   THE PURE PHASE.                                                    C
C                                                                      C
C       TREF  =  REFERENCE TEMPERATURE (298.15 K)                      C
C       PREF  =  REFERENCE PRESSURE (ONE BAR)                          C
C       T(J)  =  TEMPERATURE (SUBSCRIPT)                               C
C       P(J,K)=  PRESSURE (SUBSCRIPT)                                  C
C                                                                      C
C       A  =                                                           C
C       B  =     MAIER-KELLY HEAT CAPACITY COEFFICIENTS                C
C       C  =     WHERE B IS X 10**3 AND C IS X 10**-5                  C
C                                                                      C
C       VREF  =  MINERAL MOLAR VOLUME AT 298.15 K AND ONE BAR          C
C                IN UNITS OF CAL/BAR/MOLE.                             C
C       GRF   =  STANDARD STATE FREE ENERGY OF FORMATION AT            C
C                298.15 K AND ONE BAR.                                 C
C       SRF   =  THIRD LAW ENTROPY AT 298.15 K AND ONE BAR.            C
C       MIN   =  PHASE NAME                                            C
C       GM    =  STANDARD STATE APPARENT FREE ENERGY OF THE PHASE      C
C                AT T AND P                                            C
C                                                                      C
C**********************************************************************C
C
C
C     G IS OBTAINED FROM DOUBLE CP INTEGRATION OVER TEMPRATURE
C     RANGE THEN INTGRATION OVER PRESSURE RANGE WITH VOLUME.
C     THE MAIER-KELLEY POWER FUNCTION COEFFICIENTS ARE FOR THE
C     FOLLOWING EQUATION  CP = A + BT + C/T**2  (NOTE POSITIVE C SIGN)
C     IF ANY PHASE TRANSITION OCCURS, INTEGRATION IS DONE FROM
C     TREF TO TTRAN WITH FISRT SET OF MAIER-KELLEY COEFFICIENTS,
C     THEN FROM TTRAN TO T(J) WITH SECOND SET OF CP COEFFICIENTS.
C     THEN THE TRANSITION PRESSURE AT T(J) IS CALCULATED FROM
C     THE DP/DT SLOPE (EITHER READ IN OR CALCULATED FROM THE
C     CLAPEYRON EQUATION). IF THE FINAL PRESSURE P(K) IS
C     HIGHER THAN THE TRANTITION PRESSURE AT T(J), PP, THEN
C     THE INTEGRATION IS MADE FROM PREF TO PP WITH VREF+VCH,
C     THEN FROM PP TO P(K) WITH VREF. IF THERE ARE TWO TRANSITIONS
C     TRAN1 AND TRAN2 AND THAT THE CALCULATED SECOND TRANSITION
C     PRESSURE PP2 IS HIGHER THAN THE FIRST TRANS. PRESSURE PP1, WHICH
C     MEANS THAN THE CLAPEYRON SLOPES HAVE CROSSED, INTEGRATION IS DONE
C     FROM PREF TO PPP (WHERE PPP IS A POINT IN BETWEEN PP1 AND PP2)
C     WITH VREF+VCH1+VCH2, THEN FROM PPP TO P(K) WITH VREF ONLY.
C
      DO 80 NK=1,ITOTP
      K=INDP(NK)
      DO 80 NJ=1,ITOTT
      J=INDT(NJ)
C
      DO 100 I=1,5
  100 GI(I)=0.0
      IF(T(J).EQ.TREF) GO TO 250
C
      IF(TTRAN2.EQ.0..OR.T(J).LT.TTRAN2) GO TO 200
C
C     SECOND TRANSITION
C     -----------------
C
      PP1=(T(J)-TTRAN1)*SLOPE1 + PREF
      PP2=(T(J)-TTRAN2)*SLOPE2 + PREF
C     INTEGRATES G AT T(J) OVER PRESSURE RANGE
C
      IF(PP2.LT.PP1) GO TO 10
      PPP=(PP1+PP2)/2.
      GI(5)=GP(VRE2,PREF,PPP) + GP(VREF,PPP,P(J,K))
      WRITE(ITTY,700) P(J,K),T(J)
      WRITE(IOUT1,700) P(J,K),T(J)
  700 FORMAT(5X,'CLAPEYRON SLOPES HAVE CROSSED  P=',F10.3,
     1 ' T=',F10.3)
      GO TO 50
C
   10 IF(P(J,K).LT.PP1) GO TO 15
      GI(5)=GP(VRE2,PREF,PP2) + GP(VRE1,PP2,PP1) + GP(VREF,PP1,P(J,K))
      GO TO 50
C
   15 IF(P(J,K).LT.PP2) GO TO 20
      GI(5)=GP(VRE2,PREF,PP2) + GP(VRE1,PP2,P(J,K))
      GO TO 50
C
   20 GI(5)=GP(VRE2,PREF,P(J,K))
C
   50 CONTINUE
C      INTEGRATES G AT P(J,K) FROM TTRAN1 TO TTRAN2 TO T(J)
C
      GI(4)=GT(T(J),TTRAN2,T(J),A6,B6,C6)
      GI(3)=GT(T(J),TTRAN1,TTRAN2,A5,B5,C5) + HTRAN2*(1.-T(J)/TTRAN2)
      GO TO 220
C
C     FIRST TRANSITION
C     ----------------
  200 CONTINUE
      IF(TTRAN1.EQ.0..OR.T(J).LT.TTRAN1) GO TO 240
C      INTEGRATES G AT T(J) OVER PRESSURE RANGE
C
      PP1=(T(J)-TTRAN1)*SLOPE1 + PREF
      IF(P(J,K).LT.PP1) GO TO 30
      GI(5)=GP(VRE1,PREF,PP1) + GP(VREF,PP1,P(J,K))
      GO TO 90
C
   30 GI(5)=GP(VRE1,PREF,P(J,K))
C
   90 CONTINUE
C      INTEGRATES G AT P(J,K) FROM TREF TO TTRAN1 TO T(J)
C
      GI(3)=GT(T(J),TTRAN1,T(J),A5,B5,C5)
  220 GI(2)=GT(T(J),TREF,TTRAN1,A,B,C) + HTRAN1*(1.-T(J)/TTRAN1)
      GO TO 300
C
C     NO TRANSITIONS
C     --------------
  240 GI(2)=GT(T(J),TREF,T(J),A,B,C)
  250 GI(5)=GP(VREF,PREF,P(J,K))
C
C
  300 GI(1)=GRF - SRF*(T(J)-TREF)
C
      GM(IND,J,K)=GI(1)+GI(2)+GI(3)+GI(4)+GI(5)
   80 CONTINUE
C
      RETURN
      END
      FUNCTION GP(V,P1,P2)
      GP=V*(P2-P1)
      RETURN
      END
      FUNCTION GT(T,T1,T2,A,B,C)
      ENTH=A*(T2-T1) + B*.5*(T2*T2-T1*T1) - C*(1./T2-1./T1)
      ENTR=A*ALOG(T2/T1) + B*(T2-T1) - C*.5*(1./(T2*T2)-1./(T1*T1))
      GT=ENTH-T*ENTR
      RETURN
      END
      FUNCTION GTG(T,T1,T2,A,B,C,D)
      ENTH=A*(T2-T1) + B*.5*(T2*T2-T1*T1) - C*(1./T2-1./T1)
     1 + D*(1./3.)*(T2*T2*T2-T1*T1*T1)
      ENTR=A*ALOG(T2/T1) + B*(T2-T1) - C*.5*(1./(T2*T2)-1./(T1*T1))
     1 + D*.5*(T2*T2-T1*T1)
      GTG=ENTH-T*ENTR
      RETURN
      END
      SUBROUTINE WRITGA
C     ----------------
C     WRITGA ECHOES THE GAS DATA (THAT WAS READ IN) ON THE OUTPUT
C
      DOUBLE PRECISION GAS,ION,MINAM
      INTEGER SPEC,SPECG
      COMMON/GIBBS/G(40,21,10),GG(999,21),GM(8,21,10)
      COMMON/DATGAS/SGREF(999),GGREF(999),AG(999),BG(999),CG(999),
     1 DG(999),HGREF(999)
      COMMON/STOCH/ ACTI(10),ACTM(8),XMOLG(8),XRI(10),XRG(8),XRM(8),
     1 XRW,ACTW,NSION,NSGAS,NSMIN,SPEC(10),SPECG(8)
      COMMON/NOMS/GAS(999,2),ION(40),MINAM(8)
      COMMON/IO/INP1,INP2,INP3,INP4,INP5,IOUT1,IOUT2,IOUT3,ITTY
C
      WRITE(IOUT1,519)
      DO 113 I=1,NSGAS
      N=SPECG(I)
  113 WRITE(IOUT1,522) GAS(N,1),GGREF(N),HGREF(N),SGREF(N),AG(N),BG(N)
     1  ,CG(N),DG(N)
      RETURN
C
  519 FORMAT(/6X,' GAS  ',10X,'GREF',10X,'HREF',7X,'SREF',11X,'A',7X,
     1 'B*10**3',3X,'C*10**-5',3X,'D*10**6',/7X,'---',12X,'----',10X,
     2 '----',7X,'----',11X,'-',7X,'-------',3X,'--------',3X'--------')
  522 FORMAT(5X,A8,5X,2(F12.0,2X),F7.3,5X,4F10.3)
      END
      SUBROUTINE README(MIN)
C
C     README READS THE MINERAL DATA FILE AND MAKES SOME NECESSARY
C     CONVERSIONS.
C
      DOUBLE PRECISION MIN
      COMMON/DATMIN/A,B,C,VREF,HRF,SRF,GRF,A5,B5,C5,A6,B6,C6,TTRAN1,
     1TTRAN2,VRE1,VRE2,HTRAN1,HTRAN2,SLOPE1,SLOPE2,TTRAN3
      COMMON/IO/INP1,INP2,INP3,INP4,INP5,IOUT1,IOUT2,IOUT3,ITTY
      DIMENSION MIN(2)
C
      READ(INP4,528)GRF,HRF,SRF,VREF,A,B,C,TTRAN1,HTRAN1,VCH1,SLOPE1,A5
     1 ,B5,C5,TTRAN2,HTRAN2,VCH2,SLOPE2,A6,B6,C6,TTRAN3
C
      WRITE(IOUT1,561) MIN(1),GRF,HRF,SRF,VREF,A,B,C
      IF(TTRAN1.EQ.0.) GO TO 46
      WRITE(IOUT1,562)
      WRITE(IOUT1,563)TTRAN1,HTRAN1,SLOPE1,VCH1,A5,B5,C5
      IF(TTRAN2.EQ.0.) GO TO 46
      WRITE(IOUT1,563)TTRAN2,HTRAN2,SLOPE2,VCH2,A6,B6,C6
C
   46 IF(SLOPE1.EQ.0.) SLOPE1 = 100000.
      IF(SLOPE2.EQ.0.) SLOPE2 = 100000.
      B = B*.001
      C = C*100000.
      B5 = B5*.001
      C5 = C5*100000.
      B6 = B6*.001
      C6 = C6*100000.
      VRE1 = VCH1 + VREF
      VRE2 = VCH2 + VRE1
      VREF = VREF*.0239
      VRE1 = VRE1*.0239
      VRE2 = VRE2*.0239
      IF(VRE1.NE.VREF.AND.HTRAN1.NE.0.) SLOPE1=HTRAN1/((VRE1-VREF)*
     1TTRAN1)
      IF(VRE2.NE.VRE1.AND.HTRAN2.NE.0.) SLOPE2 = HTRAN2/((VRE2-VRE1)*
     1TTRAN2)
      RETURN
C
  561 FORMAT(5X,A8,5X,2(F12.0,2X),F7.3,1X,F10.2,5X,3(F10.3,1X))
  528 FORMAT(2F15.0,2F15.5/7F10.3/7F10.3/7F10.3)
  562 FORMAT(18X,'TRANSITION T',4X,'H TRANS',3X,'EXP DP/DT',2X,
     1 'V CHANGE')
  563 FORMAT(13X,'----- ',F7.2,' K',4X,F12.0,2X,F6.2,2X,F10.3,5X,
     1 3(F10.3,1X),3X,E12.6)
      END
      SUBROUTINE REAC
C
C     THIS SUBROUTINE COMPUTES LOG K'S AND LOG Q'S FOR THE REACTIONS. IT
C     IS CALLED ONLY AFTER ALL THE NECESSARY FREE ENERGIES HAVE BEEN
C     CALCULATED
C
      DOUBLE PRECISION GAS,ION,MINAM
      INTEGER SPEC,SPECG
      COMMON/DATA/TREF,PREF,T(24),P(24,10),TC(24)
      COMMON/MAX/JMAX,KMAX,INDP(10),INDT(24),ITOTT,ITOTP,MAXG
      COMMON/GIBBS/G(40,21,10),GG(999,21),GM(8,21,10)
      COMMON/STOCH/ ACTI(10),ACTM(8),XMOLG(8),XRI(10),XRG(8),XRM(8),
     1 XRW,ACTW,NSION,NSGAS,NSMIN,SPEC(10),SPECG(8)
      COMMON/NOMS/GAS(999,2),ION(40),MINAM(8)
      COMMON/IO/INP1,INP2,INP3,INP4,INP5,IOUT1,IOUT2,IOUT3,ITTY
      COMMON/EQUIL/XLOGK(24)
C
C     PRINTS G MINERAL ON OUTPUT
C
      IF(NSMIN.EQ.0) GOTO 100
      WRITE(IOUT1,620)
      DO 150 NK=1,ITOTP
      K=INDP(NK)
      WRITE(IOUT1,620)
      WRITE(IOUT1,630) (MINAM(I), I=1,NSMIN)
      DO 155 NJ=1,ITOTT
      J=INDT(NJ)
      WRITE(IOUT1,635) TC(J),P(J,K),(GM(I,J,K), I=1,NSMIN)
  155 CONTINUE
  150 CONTINUE
C
C     CALCULATES DELTA G, LOGK AND LOG Q FOR THE GIVEN REACTION
C     --------------------------------------------------------
C
C       WRITES FIRST THE STOCHIOMETRY ON OUTPUT
  100 WRITE(IOUT1,750)
  750 FORMAT(///6X,'REACTION STOCHIOMETRY:',/6X,21('-'),
     1 /23X,'(ACT.)')
C
  205 IF(NSGAS.EQ.0) GO TO 215
      DO 210 I=1,NSGAS
  210 WRITE(IOUT1,615) XRG(I),GAS(SPECG(I),1),XMOLG(I)
C
  215 IF(NSMIN.EQ.0) GO TO 225
      DO 220 I=1,NSMIN
  220 WRITE(IOUT1,615) XRM(I),MINAM(I),ACTM(I)
C
C
  225 CONTINUE
      R=1.98719
      DO 400 NK=1,ITOTP
      K=INDP(NK)
      WRITE(IOUT1,752)
  752 FORMAT(///6X,'DEG.C',5X,'BARS',4X,'DELTA G REAC',4X,'LOG K',
     1 11X,'DELTA G REAC',4X,'LOG Q',/6X,'-----',5X,'----',7X,
     2'ACT = 1',6X,'-----',11X,'------------',4X,'-----',/27X,'-------')
      DO 410 NJ=1,ITOTT
      J=INDT(NJ)
      DELG=0.
      DELA=0.
  305 IF(NSGAS.EQ.0) GO TO 315
      DO 310 I=1,NSGAS
      NG=SPECG(I)
      DELG=DELG+XRG(I)*GG(NG,J)
  310 DELA=DELA+R*T(J)*ALOG(XMOLG(I))
C
  315 IF(NSMIN.EQ.0) GO TO 325
      DO 320 I=1,NSMIN
      DELG=DELG+XRM(I)*GM(I,J,K)
  320 DELA=DELA+R*T(J)*ALOG(ACTM(I))
C
  325 CONTINUE
      XLOGK(NJ)=-DELG/(2.302585*R*T(J))
      DELGQ=DELG+DELA
      XLOGQ=-DELGQ/(2.302585*R*T(J))
C
      WRITE(IOUT1,625) TC(J),P(J,K),DELG,XLOGK(NJ),DELGQ,XLOGQ
  410 CONTINUE
  400 CONTINUE
      RETURN
C
  617 FORMAT(5X,F7.3,' :H2O       ',F4.2)
  615 FORMAT(5X,F7.3,1X,':',A8,2X,F4.2)
  625 FORMAT(6X,F5.0,2X,F9.4,3X,2(E12.5,3X,F8.3,9X))
  620 FORMAT(///)
  630 FORMAT(5X,'APPARENT FREE ENERGIES OF MINERALS AT GIVEN P AND T:'/
     1 ,5X,50('-')/,1X,'DEG.C',4X,'BARS',5X,A8,7(6X,A8))
  635 FORMAT(1X,F5.0,1X,F9.4,1X,8(F12.0,2X))
      END
      SUBROUTINE READER
C
C     READER IS A ROUTINE TO READ MINERAL AND GAS STOCHIOMETRIES
C     IN THE FILE GASTHERM. THE DIFFERENT STOCHIOMETRIES ARE THEN
C     USED TO CALCULATE LOG K'S FROM THE THERMO DATA CONTAINED IN
C     GAS AND MINERAL DATA FILES. THE CALCULATED LOG K'S ARE COMPARED
C     TO THE LOG K'S FROM GASTHERM AND THE DIFFERENCE BETWEEN THE TWO
C     VALUES IS INCLUDED IN THE OUTPUT. THE ROUTINE ALSO CREATES
C     A COPY OF GASTHERM WITH THE CALCULATED NEW VALUES.
C
      DOUBLE PRECISION NAME,MIN,GAS,MINAM,ION,BLANK,DUMMY
     1 ,ETOIL,TEXT,REF,REF1
      INTEGER SPEC,SPECG,SAQ,BSTOT
      COMMON/MAX/JMAX,KMAX,INDP(10),INDT(24),ITOTT,ITOTP,MAXG
      COMMON/STOCH/ ACTI(10),ACTM(8),XMOLG(8),XRI(10),XRG(8),XRM(8),
     1 XRW,ACTW,NSION,NSGAS,NSMIN,SPEC(10),SPECG(8)
      COMMON/NOMS/GAS(999,2),ION(40),MINAM(8)
      COMMON/IO/INP1,INP2,INP3,INP4,INP5,IOUT1,IOUT2,IOUT3,ITTY
      COMMON/EQUIL/XLOGK(24)
      DIMENSION NAME(999),ITOT(999),
     1 IND(999,8),COEF(999,8),AKLOG(999,13),IFLAG(999),SAQ(50),
     2 MIN(2),WT(999),DELK(13),TEXT(10),REF(999),REF1(999)
C       SAQ ARE THE INDEXES FOR THE SAME SPECIES IN GAS DATA FILE
      DATA POS/'Y'/
      DATA BLANK/'        '/
      DATA ETOIL/'********'/
C      MAXG IS THE TOTAL OF GASES IN GAS DATA (TRANSFERED FROM MAIN)
C      BSTOT IS THE TOTAL NUMBER OF BASIC SPECIES IN GASTHERM
      BSTOT=42
      DO 13 J = 1, BSTOT
      SAQ(J) = J
   13 CONTINUE
C
C     THE FOLLOWING SKIPS THE FIRST PART OF GASTHERM DATA AND
C     POSITION READING AT THE BEGINNING OF THE GAS DATA
C      FIRST, WE SKIP TO THE FIRST DERIVED  GAS SPECIES
 1000 CONTINUE
   14 READ(INP5,803) (TEXT(I), I=1,10)
      WRITE(IOUT2,803) (TEXT(J), J=1,10)
      DUMMY = TEXT(1)
      IF(DUMMY.NE.ETOIL)GO TO 14
      DO 15 II=1,BSTOT
      READ(INP5,803) (TEXT(I), I=1,10)
   15 WRITE(IOUT2,803) (TEXT(J), J=1,10)
C      NOW WE SHOULD BE AT THE TOP OF THE GAS DATA
      WRITE(ITTY,501)
  501 FORMAT(///5X,'TOP PART OF FILE HAS BEEN READ IN GASTHERM...'/
     1 5X,'I AM NOW STARTING TO GAS DATA.'/)
C
      WRITE(ITTY,504)
  504 FORMAT(///5X,'DO YOU WANT ME TO WRITE A GASTHERM FORMATTED OUTPUT'
     1,' FILE',/5X,'WITH NEW LOG K VALUES ?'//)
      READ(ITTY,827) SOLANS
C
      WRITE(ITTY,502)
  502 FORMAT(///5X,'DO YOU WANT TO WORK ON GASES(1) OR MINERALS(2) ?'/)
      READ(ITTY,*) IMIN
C     SKIPS GASES IN GASTHERM IF IMIN EQUALS 2
      IF(IMIN.EQ.1) GO TO 96
   95 READ(INP5,803) (TEXT(I), I=1,10)
      WRITE(IOUT2,803) (TEXT(J), J=1,10)
      DUMMY = TEXT(1)
      IF(DUMMY.NE.BLANK)GO TO 95
      GO TO 407
C
   96 CONTINUE
C
C     GET ALL THE GASES OR MINERALS FROM GASTHERM
C     -------------------------------------------
  300 CONTINUE
C      READS ALL THE GASES
       NS = 0
   50 NS = NS + 1
      READ(INP5,820) NAME(NS),IT,(COEF(NS,I),IND(NS,I),I=1,IT)
      DUMMY = NAME(NS)
      IF(DUMMY.EQ.BLANK)GO TO 53
      ITOT(NS)=IT
      READ(INP5,825) (AKLOG(NS,I), I=1,8)
      READ(INP5,825) (AKLOG(NS,I), I=9,13)
      READ(INP5,825) DUM
      GO TO 50
   53 CONTINUE
      NNAM = NS - 1
      IF(IMIN.EQ.1) GO TO 410
C
  407 CONTINUE
C     READS ALL THE MINERALS
      DO 55 NS=1,999
      READ(INP5,822,END=409) NAME(NS),WT(NS),IT,(COEF(NS,I),IND(NS,I),
     1 I=1,IT)
      DUMMY = NAME(NS)
      IF(DUMMY.EQ.BLANK)GO TO 409
      ITOT(NS)=IT
      READ(INP5,825) (AKLOG(NS,I), I=1,8),REF1(NS)
      READ(INP5,825) (AKLOG(NS,I), I=9,13)
      READ(INP5,825) DUM
   55 CONTINUE
  409 NNAM=NS-1
C
C
C
  410 CONTINUE
C
C     NOW WE TRY TO MATCH THE MINERALS OR GASES WITH THE ONES
C     IN THE FILE GAS DATA AND THE FILE MIN DATA. MATCHING IS
C     DONE BY NAME, SO BE SURE THAT THE FIRST 8 CHARACTERS OF
C     NAMES ARE SIMILAR IN ALL THE FILES.
C
      DO 419 I=1,999
  419 IFLAG(I)=0
      IF(IMIN.GT.1) GO TO 450
C
C     LOOKS FOR THE GAS DATA AND CALCULATES K FOR THE REACTION
C     THE GAS THERMO DATA WAS READ IN MAIN
C     ------------------------------------
      IGAS=0
      DO 420 NS=1,NNAM
      DO 421 NG=1,MAXG
      IF(NAME(NS).EQ.GAS(NG,1)) GO TO 422
  421 CONTINUE
      GO TO 425
  422 CONTINUE
C
      IFLAG(NS)=NS
  425 WRITE(IOUT1,802)
      IGAS=IGAS+1
      NSGAS=ITOT(NS)
      XRG(1) = -1.
      SPECG(1) = NG
      XMOLG(1) = 1.
      DO 424 I=2,NSGAS
      XRG(I)=COEF(NS,I)
      KKK=IND(NS,I)
      SPECG(I)=SAQ(KKK)
      XMOLG(I) = 1.
  424 CONTINUE
C
      CALL WRITGA
      CALL GASES
      CALL REAC
      IF(SOLANS.NE.POS) GO TO 423
C      WRITES GASTHERM FORMATTED FILE WITH NEW LOG K'S AND DISPLAYS
C      THE DIFFERENCE BETWEEN OLD AND NEW LOG K'S VALUES
      IT=ITOT(NS)
      WRITE(IOUT2,821) NAME(NS),IT,(COEF(NS,I),IND(NS,I),I=2,IT)
      WRITE(IOUT2,826) NAME(NS),(XLOGK(I), I=1,8),GAS(NG,2)
      WRITE(IOUT2,826) NAME(NS),(XLOGK(I), I=9,13)
      WRITE(IOUT2,800) NAME(NS)
      DO 430 J=1,13
  430 DELK(J)= XLOGK(J)-AKLOG(NS,J)
      WRITE(IOUT3,826) NAME(NS),(DELK(J), J=1,8)
      WRITE(IOUT3,826) NAME(NS),(DELK(J), J=9,13)
C
C
  423 IF(IGAS.NE.NNAM) GO TO 420
      WRITE(ITTY,522)
  522 FORMAT(///5X,'OKAY.... I FINISHED THE CALCULATION !'/)
      GO TO 426
  420 CONTINUE
C      WHAT GAS WAS NOT FOUND ??
      WRITE(ITTY,522)
      DO 428 NS=1,NNAM
  428 IF(IFLAG(NS).EQ.0) WRITE(ITTY,520) NAME(NS)
  520 FORMAT(5X,'I COULD NOT FIND ',A8,' IN FILE GAS DATA')
C
  426 CONTINUE
      WRITE(IOUT2,803) BLANK
      WRITE(ITTY,525)
  525 FORMAT(///5X,'DO YOU WANT TO CALCULATE DATA FOR MINERALS NOW?',/)
      READ(ITTY,827) ANS
      IF(ANS.EQ.POS) GO TO 572
  570 READ(INP5,803,END=571) (TEXT(I), I=1,10)
      WRITE(IOUT2,803) (TEXT(J), J=1,10)
      GO TO 570
  571 RETURN
  572 IMIN = 2
      GO TO 407
C
C     LOOKS FOR THE MINERAL DATA AND CALCULATES K FOR THE REACTION
C     THE MINERAL THERMO DATA IS READ BY CALL TO README
C     -------------------------------------------------
C       THE DATA FILE IS SCANNED UNTIL WE REACH THE FIRST SELECTED
C       MINERAL. THEN WE SCAN FORWARD TO LOOK FOR THE NEXT CHOSEN
C       MINERAL. IF WE PASSED THE SECOND MINERAL BEFORE FINDING THE
C       FIRST ONE, THEN WE REWIND AND LOOK AGAIN. WE SAVE TIME IF THE
C       THE MINERALS ARE IN THE SAME ORDER OF APPARITION IN MINERAL DATA
C       AND IN GASTHERM DATA. THE USER CAN INPUT THE NAMES IN ANY ORDER
C       SINCE NAME(NS) WILL AUTOMATICALLY CONTAIN  THE NAMES IN THE
C       SAME ORDER OF APPARITION THEY ARE IN GASTHERM.
  450 CONTINUE
C
  451 NS = 1
      IREW=0
      READ(INP4,805)
  452 CONTINUE
      READ(INP4,806,END=455) MIN(1),REF(NS)
      IF(NAME(NS).EQ.MIN(1)) GO TO 453
      READ(INP4,807,END=455)
      GO TO 452
C
  453 WRITE(ITTY,530) MIN(1)
  530 FORMAT(5X,'FOUND ',A8,' THERMO DATA IN FILE MIN DATA')
      WRITE(IOUT1,560)
  560 FORMAT(//1H1,5X,'MINERAL',9X,'GREF',10X,'HREF',7X,'SREF',6X,'VREF'
     1 ,13X,'A',7X,'B*10**3',3X,'C*10**-5'/,6X,'-------',9X,'----',10X,
     2 '----',7X,'----',6X,'----',13X,'-',7X,'------',3X,'--------')
      CALL README(MIN)
C
      NSGAS=ITOT(NS)
      DO 456 I=1,NSGAS
      XRG(I)=COEF(NS,I)
      SPECG(I)=IND(NS,I)
  456 XMOLG(I) = 1.
      CALL WRITGA
      CALL GASES
  458 NSMIN=1
      XRM(1)=-1.
      ACTM(1)=1.
      MINAM(1)=NAME(NS)
      INDEX=1
      CALL SOLIDS(INDEX)
      CALL REAC
      IF(SOLANS.NE.POS) GO TO 457
C      WRITES GASTHERM FORMATTED FILE WITH NEW LOG K'S AND DISPLAYS
C      THE DIFFERENCE BETWEEN OLD AND NEW LOG K'S VALUES
      IT=ITOT(NS)
      WRITE(IOUT2,822) NAME(NS),WT(NS),IT,(COEF(NS,I),IND(NS,I),I=1,IT)
      WRITE(IOUT2,826) NAME(NS),(XLOGK(J), J=1,8),REF(NS)
      WRITE(IOUT2,826) NAME(NS),(XLOGK(J), J=9,13)
      WRITE(IOUT2,826) NAME(NS)
      DO 435 J=1,13
  435 DELK(J)= XLOGK(J)-AKLOG(NS,J)
      WRITE(IOUT3,826) NAME(NS),(DELK(J), J=1,8)
      WRITE(IOUT3,826) NAME(NS),(DELK(J), J=9,13)
C
C
  457 NS=NS+1
      IREW = 0
      IF(NS.GT.NNAM) GO TO 460
      GO TO 452
C      IF A MINERAL IS NOT FOUND AFTER THE FISRT SCAN OF THE FILE THEN
C      WE REWIND AND LOOK AGAIN BECAUSE WE MAY HAVE MISSED THE FIRST
C      PART OF THE FILE. IF WE STILL CANT FIND IT THEN WE REWIND AGAIN
C      AND LOOK FOR THE NEXT CHOSEN MINERAL.
  455 REWIND INP4
      IREW=IREW+1
      READ(INP4,805)
      IF(IREW.GT.1) GO TO 459
      GO TO 452
C
  459 WRITE(ITTY,532) NAME(NS)
  532 FORMAT(//5X,'I COULD NOT FIND ',A8,' IN FILE MINERAL DATA')
      IF(SOLANS.NE.POS) GO TO 534
C      SINCE THE MINERALS COULD NOT BE FOUND IN THE MINERAL DATA FILE
C      THE PROGRAM WRITES THE OLD LOG K'S INTO THE NEW FILE INSTEAD
      IT=ITOT(NS)
      WRITE(IOUT2,822) NAME(NS),WT(NS),IT,(COEF(NS,I),IND(NS,I),I=1,IT)
      WRITE(IOUT2,826) NAME(NS),(AKLOG(NS,J), J=1,8),REF1(NS)
      WRITE(IOUT2,826) NAME(NS),(AKLOG(NS,J), J=9,13)
      WRITE(IOUT2,826) NAME(NS)
  534 NS=NS+1
      IF(NS.GT.NNAM) GO TO 460
      IREW = 0
      GO TO 452
C
C
  460 WRITE(ITTY,522)
      RETURN
C
  800 FORMAT(A8)
  801 FORMAT(/A8)
  802 FORMAT(1H1)
  803 FORMAT(10(A8))
  805 FORMAT(/)
  806 FORMAT(A8,52X,A8)
  807 FORMAT(///)
  820 FORMAT(A8,8X,I1,7(F6.3,I2))
  821 FORMAT(A8,8X,I1,8X,7(F6.3,I2))
  822 FORMAT(A8,2X,F7.2,I2,1X,7(F6.3,I2)/)
  825 FORMAT(8X,8(F8.3),A8)
  826 FORMAT(A8,8(F8.3),A8)
  827 FORMAT(A1)
  830 FORMAT(10A8)
      END
