C     MAIN PROGRAM NASA LEWIS EQUILIBRIUM. MODIFIED FOR MAC BY JRH
C       BEGINNING IN MAY, 1987
C     LAST MODIFICATION SEPT 22, 1987
C     Modified 09/22/03 by R.L. xxxxx
C     Commented out unused variables. Replaced TTY and Unit 9 with "*" 
C     to be compatable with Microsoft FORTRAN.  Changed read format in
C     subroutine SEARCH to read each record separately so no error on
C     last record (old version probably used double buffering). Changed
C     write format for namelist output in main program.  Commented out 
C     bypassed code in subroutine FCOEFF to prevent "go to" end of do 
C     loop.  Change test for FALSE to use ".NOT." instead of "FALSE" 
C     in subroutine OUT1. Also changed "go to 15" to "go to 1115" so
C     execution does not jump back into do loop in subroutine OUT1. 
C     Put four variables in local common block in subroutine EQLBRM
C     because they use the previous value they were set to in the
C     routine.
C
      DOUBLE PRECISION G,X
C
C#### F is not used
C      DOUBLE PRECISION COEF,S,EN,ENLN,HO,F
      DOUBLE PRECISION COEF,S,EN,ENLN,HO
C
C#### TTY replaced by "*" for Microsoft's FORTRAN
C      INTEGER DATA,ENSERT,REAC,STOP,OMIT,DATOUT,TTY
      INTEGER DATA,ENSERT,REAC,STOP,OMIT,DATOUT
      LOGICAL ALTER,SYSTM,PLUS,TEM,GAS,NEW,OS
        CHARACTER*20 INFILE,OUTFILE,THERMOFILE
C
      DIMENSION NCD(4),DAT(34),DATA(17)
C
      COMMON/POINTS/P(45),T(45),PPP(9),TTT(9),TOTN(9),RSIZE,TEM,GAS,NEW,
     1 GTOTN(9),FREE(9),OXFUG(9),FO2C(9),FO2H(9)
      COMMON/SPECES/COEF(2,7,100),S(100),EN(100,9),ENLN(100),HO(100)
     1,SUB(100,8),IUSE(100),TEMP(100,2),OMIT(8,100),ENSERT(8,100)
     2,JUSE(100),A(15,100)
      COMMON/MISC/ENN,SUMN,TT,LLMT(15),BO(15),TM,PP,R,ENNL,TRI,BASE(15)
     1,NAME(15,5),ANUM(15,5),TOP(15),ADD(15),VOL(100),GAMMA(100),NOCONV
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX
      COMMON /DOUBLE/ G(20,21), X(20)
C
      DATA MIT/4HOMIT/,REAC/4HREAC/,NMLT/4HNAME/,INSERT/4HINSE/
C#### END is not used
C     1 ,END/4HENDD/,STOP/4HSTOP/,LANK/1H /
     1 ,STOP/4HSTOP/,LANK/1H /
C#### ITY is now used as "*" and THERMOFILE changed to THERMO.DAT
C        DATA DATOUT/4/,TTY/9/, THERMOFILE/'THERMODAT'/
        DATA DATOUT/4/, THERMOFILE/'THERMO.DAT'/
C
    1 CONTINUE 
C#### TTY is not used, now used as "*"
C        WRITE(TTY,3200)
        WRITE(*,3200)
3200    FORMAT(' TYPE 1 TO CONTINUE, ZERO TO STOP')
C#### TTY is not used, now used as "*"
C        READ(TTY,*) NSTOP
        READ(*,*) NSTOP
        IF(NSTOP.NE.1) GOTO 111
C#### TTY is not used, now used as "*"
C        WRITE(TTY,3334)
        WRITE(*,3334)
3334    FORMAT(' Input name of INPUT file ',
     1' (the thermofile must be named THERMODAT)')
C#### TTY is not used, now used as "*"
C        READ(TTY,*) INFILE
        READ(*,8787) INFILE  
8787    FORMAT(A20)
C#### TTY is not used, now used as "*"
C        WRITE(TTY,3335)
        WRITE(*,3335)
3335    FORMAT(' Input name of OUTPUT file ')
C#### TTY is not used, now used as "*"
C        READ(TTY,*) OUTFILE
        READ(*,8787) OUTFILE
      INQUIRE(FILE=INFILE,OPENED=OS)
      IF(OS) CLOSE(5)
      INQUIRE(FILE=OUTFILE,OPENED=OS)
      IF(OS) CLOSE(6)
      INQUIRE(FILE=THERMOFILE,OPENED=OS)
      IF(OS) CLOSE(4)
      OPEN (UNIT=6,FILE=OUTFILE)
      OPEN (UNIT=5,FILE=INFILE)
      OPEN (UNIT=4,FILE=THERMOFILE)
      OPEN (UNIT=3,STATUS='SCRATCH')
      NEWMOL = 0
      NOCONV = 0
      TEM=.FALSE.
      ALTER = .FALSE.
      SYSTM = .FALSE.
      PLUS = .FALSE.
      GAS = .FALSE.
      RSIZE = 0.0
      NOMIT = 0
      R = 8314.3/4184.
C
C     READ AND WRITE CODE CARDS.
C
  203 READ (5,204) (DATA(I),I=1,17)
  204 FORMAT(A4,4X,8A4,6X,8A4)
      WRITE(6,2045)(DATA(I),I=1,17)
 2045 FORMAT(1X,A4,4X,2(8A4,6X))
      IF(DATA(1).EQ.REAC)  GO TO 11
      IF (DATA(1).EQ.MIT) GO TO 205
      IF (DATA(1).EQ.INSERT) GO TO 180
      IF(DATA(1).EQ.NMLT)  GO TO 210
      IF(DATA(1).EQ.STOP) GO TO 111
C#### TTY is not used, now used as "*"
C     WRITE(TTY,1024) (DATA(I),I=1,17)
      WRITE(*,1024) (DATA(I),I=1,17)
 1024 FORMAT(' ',17A4/' ERROR IN ABOVE CARD, CONTENTS IGNORED.')
      GO TO 203
C
C     CALL SUBROUTINE REACT TO READ REACTANT CARDS.
C
   11 NSERT = 0
      NMIN = 0
      CALL REACT
      GO TO 203
C
   18 READ(5,20)(DAT(I),I=1,5),NCD(2),(DAT(J),J=6,10),NCD(3),(DAT(K),
     1K=11,14),NCD(4)
   20 FORMAT(5E15.8,I5/5E15.8,I5/4E15.8,I20)
      WRITE (DATOUT,21)(DAT(I),I=1,14)
   21 FORMAT(5E15.8/5E15.8/4E15.8)
C     CHECK INSERT CARDS.
C
  180 DO 186 I=1,9,8
      IF(DATA(I+1).EQ.LANK) GO TO 186
      NSERT = NSERT+1
      DO 185 J=1,8
      ENSERT(J,NSERT)=DATA(I+J)
  185 CONTINUE
  186 CONTINUE
      GO TO 203
C
C     CHECK OMIT CARDS
C
  205 DO 209 I=1,9,8
      IF(DATA(I+1).EQ.LANK) GO TO 209
      NOMIT = NOMIT+1
      DO 208 J=1,8
      OMIT(J,NOMIT)=DATA(I+J)
  208 CONTINUE
  209 CONTINUE
      GO TO 203
C
C     VALUE,SYSTEM INCREMENT, TEM AND GAS OPTIONS.
C
  210 READ(5,77) TEM,GAS,KASE,RSIZE,TRI,(T(I),I=1,45),(P(J),J=1,45)
   77 FORMAT(2L1,I5,F11.0,F9.0/10(9F5.0/))
      GAS=.TRUE.
      WRITE(6,78) KASE,RSIZE,TRI,TEM,GAS,(T(I),I=1,45),(P(J),J=1,45)
C#### Error in format, wrong format for RSIZE
C   78 FORMAT(/13H CASE NUMBER ,I6,5X,6HTEM = ,L1,5X,6HTRI = ,F10.0//
   78 FORMAT(/13H CASE NUMBER ,I6,5X,8HRSIZE = ,F11.2,5X,6HTRI = ,F10.0//
     2 8H TEMP. =,2X,L1,2X,8H GAS   =,2X,L1//8H T     =,6F10.0/,6(8X,
     36F10.0/),8X,3F10.0//,' P     =',6F10.0/,6(8X,6F10.0/),8X,3F10.0/)
      DO 304 I=1,45
      IF(T(I).LT.0.0001) GO TO 304
      T(I) = T(I) + 273.15
  304 CONTINUE
      IF(NLM.EQ.0) GO TO 1
C
C     SET NP AND NT VALUES.
C
      DO 305 I=1,45
      IF(P(I ).EQ.0.) GO TO 322
      NP = I
  305 CONTINUE
  322 DO 307 IT = 1,45
      IF(T(IT).EQ.0.) GO TO 748
      NT = IT
  307 CONTINUE
  748 IF(PLUS.AND.TRI.GT.0.) GO TO 763
C  748 CONTINUE
C
C     CALL SUBROUTINE SEARCH TO STORE CURRENT THERMO DATA.
C
C#### TTY (Unit 9) is not used, now used as "*"
C        WRITE(9,3367)
      WRITE(*,3367)
3367  FORMAT(' CALLING SEARCH')
      CALL SEARCH
      IF(NS.EQ.0)  GO TO 1
C
C     INITIAL ESTIMATES.
C
      NPT = 1
      CALL INIT
C
C     CALL SUBROUTINE THERMO TO SET UP CURRENT TEMP AND PRESSURE.
C
  763 CALL THERMP
      GO TO 1
  111 CONTINUE
C#### TTY (Unit 9) is not used, now used as "*"
C        WRITE(TTY,3340)
        WRITE(*,3340)
3340    FORMAT(' CONGRATULATIONS! The program has ended normally ')
        STOP
      END
C
C This routine forces the compiler to expand the heap by 2k bytes so that
C the spool subroutine can be loaded by the print command.
C If this heap space is no allocated in this manner or with the linker
C the program will give an error 64 (out of memory) when the print option is
C selected from the menu.
C
        subroutine dum
        integer*2 ary(6000)
        common ary
        return
        end
      SUBROUTINE REACT
C
C     READ REACTANT CARDS AND COMPUTE MOLE FRACTIONS.
C
      DIMENSION ATOM(101),IATOM(101),DATA(15)
      LOGICAL ALTER,SYSTM,PLUS
C
      COMMON/MISC/ENN,SUMN,TT,LLMT(15),BO(15),TM,PP,R,ENNL,TRI,BASE(15)
     1,NAME(15,5),ANUM(15,5),TOP(15),ADD(15),VOL(100),GAMMA(100),NOCONV
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX
C
C     ATOMIC SYMBOLS.
C
      DATA  IATOM/
     A 2HH , 2HHE, 2HLI, 2HBE, 2HB , 2HC , 2HN , 2HO , 2HF , 2HNE, 2HNA,
     B 2HMG, 2HAL, 2HSI, 2HP , 2HS , 2HCL, 2HAR, 2HK , 2HCA, 2HSC, 2HTI,
     C 2HV , 2HCR, 2HMN, 2HFE, 2HCO, 2HNI, 2HCU, 2HZN, 2HGA, 2HGE, 2HAS,
     D 2HSE, 2HBR, 2HKR, 2HRB, 2HSR, 2HY , 2HZR, 2HNB, 2HMO, 2HTC, 2HRU,
     E 2HRH, 2HPD, 2HAG, 2HCD, 2HIN, 2HSN, 2HSB, 2HTE, 2HI , 2HXE, 2HCS,
     F 2HBA, 2HLA, 2HCE, 2HPR, 2HND, 2HPM, 2HSM, 2HEU, 2HGD, 2HTB, 2HDY,
     G 2HHO, 2HER, 2HTM, 2HYB, 2HLU, 2HHF, 2HTA, 2HW , 2HRE, 2HOS, 2HIR,
     I 2HPT, 2HAU, 2HHG, 2HTL, 2HPB, 2HBI, 2HPO, 2HAT, 2HRN, 2HFR, 2HRA,
     J 2HAC, 2HTH, 2HPA, 2HU , 2HNP, 2HPU, 2HAM, 2HCM, 2HBK, 2HCF, 2HES,
     K 2HFM, 2HD /
C
C     ATOMIC WEIGHTS.
C
      DATA ATOM/
     L  1.00797,  4.0026,   6.939,  9.0122,  10.811, 12.01115, 14.0067,
     M  15.9994, 18.9984,  20.183, 22.9898,  24.312,  26.9815,  28.086,
     N  30.9738,  32.064,  35.453,  39.948,  39.102,  40.0800,  44.956,
     O  47.90  ,  50.942,  51.996,  54.938,  55.847,  58.9332,  58.710,
     P  63.54  ,  65.370,  69.720,  72.590, 74.9216,  78.96  ,  79.909,
     Q  83.80  ,  85.470,  87.620,  88.905,  91.220,  92.906 ,  95.940,
     R  99.00  , 101.070, 102.905, 106.400, 107.870, 112.400 , 114.820,
     S 118.69  , 121.750, 127.600,126.9044, 131.300, 132.905 , 137.340,
     T 138.910 , 140.120, 140.907, 144.240, 145.000, 150.350 , 151.960,
     U 157.250 , 158.924, 162.500, 164.930, 167.260, 168.934 , 173.040,
     V 174.997 , 178.490, 180.948, 183.850, 186.200, 190.200 , 192.200,
     W 195.090 , 196.967, 200.590, 204.370, 207.190, 208.980 , 210.000,
     X 210.000 , 222.000, 223.000, 226.000, 227.000, 232.038 , 231.000,
     Y 238.030 , 237.000, 242.000, 243.000, 247.000, 249.000 , 251.000,
     Z 254.000 , 253.000, 2.014102/
C
      DATA LANK/1H /
C
      WP = 0.
      DO 8 J=1,15
      DATA(J) = 0.
      LLMT(J)=0
      BO(J)=0.
    8 CONTINUE
      N=1
      NLM=1
C
C     READ AND WRITE REACTANT CARDS
C
   20 IF(NEWMOL.NE.0) GO TO 80
      READ(5,21) (NAME(N,I),ANUM(N,I),I=1,5),BASE(N),TOP(N),ADD(N)
   21 FORMAT(5(A2,F7.0),3F9.0)
      IF(NAME(N,1).EQ.LANK) GO TO 200
      IF(NLM.EQ.0) GO TO 20
      WRITE(6,31) (NAME(N,I),ANUM(N,I),I=1,5),BASE(N),TOP(N),ADD(N)
   31 FORMAT(1X,5(A2,1X,F7.4),3F7.4)
      GO TO 85
   80 IF(NAME(N,1).EQ.LANK) GO TO 200
   85 RM = 0.
C
C     STORE ATOMIC SYMBOLS IN LLMT ARRAY.
C     CALCULATE MOLECULAR WEIGHT.
C
      DO 100 JJ=1,5
      IF(ANUM(N,JJ).EQ.0.)GO TO 101
      DO 41 J=1,15
      NJ = J
      IF(LLMT(J).EQ.0) GO TO 45
      IF(BASE(N).EQ.0.) GO TO 41
      IF(NAME(N,JJ).EQ.LLMT(J))GO TO 46
   41 CONTINUE
   45 IF(BASE(N).EQ.0.) NMIN = NMIN+1
      IF(BASE(N).EQ.0.) GO TO 100
      NLM = NJ 
      LLMT(NJ) = NAME(N,JJ)
   46 DO 48 KK=1,101
      KT = KK
      IF(IATOM(KK).EQ.NAME(N,JJ)) GO TO 50
   48 CONTINUE
C#### TTY (Unit 9) is not used, now used as "*"
C      WRITE(9,90) N
      WRITE(*,90) N
   90 FORMAT(1X,26HERROR IN REACTANT CARD NO.,I4)
      NLM = 0
      GO TO 20
   50 RM = RM+ANUM(N,JJ)*ATOM(KT)*BASE(N)
      DATA(NJ) = DATA(NJ) + (ANUM(N,JJ) * BASE(N))
  100 CONTINUE
C
C     ADD CONTRIBUTIONS TO WP.
C
  101 IF(BASE(N).EQ.0.) GO TO 300
      WP = WP+RM
  300 N=N+1
      IF(N.NE.16) GO TO 20
  200 NREAC =N-1
      IF(NEWMOL.NE.0) GO TO 600
C
C     DETERMINE TYPE OF SYSTEM BEING STUDIED.
C
      DO 555 I = 1,NREAC
      IF(BASE(I).EQ.0..AND.TOP(I).EQ.0.) GO TO 650
      IF(BASE(I).EQ.0.) GO TO 900
      GO TO 555
  900 PLUS = .TRUE.
      GO TO 600
  650 SYSTM = .TRUE.
      GO TO 600
  555 CONTINUE
      ALTER = .TRUE.
C
C     COMPUTE MOLE FRACTION OF EACH ELEMENT.
C
  600 DO 215 J=1,NLM
      BO(J) = DATA(J)/WP
  215 CONTINUE
C#### TTY (Unit 9) is not used, now used as "*"
C      WRITE(9,3333)
      WRITE(*,3333)
3333    FORMAT(' Normal exit from subroutine REACT')
        RETURN
      END
        SUBROUTINE SEARCH
C
C     SEARCH TAPE FOR THERMO DATA FOR SPECIES TO BE CONSIDERED.
C
C     THE FOLLOWING DOUBLE PRECISION TYPE STATEMENTS ARE REQUIRED FOR
C     IBM 360 MACHINE ONLY.
C
C#### F is not used
C      DOUBLE PRECISION COEF,S,EN,ENLN,HO,F
      DOUBLE PRECISION COEF,S,EN,ENLN,HO
C
C#### TTY is not used, now used as "*"
C      INTEGER SUB,END,TOOBIG,TSUB,ENSERT,DATE,TDATE,PHAZ,GAS,OMIT,TTY
      INTEGER SUB,END,TOOBIG,TSUB,ENSERT,DATE,TDATE,PHAZ,GAS,OMIT
      LOGICAL PLUS,TEST
C
      DIMENSION DATE(2,100),MT(10),B(10),NAM(8),TOOBIG(8,100),LMT(15)
      DIMENSION TDATE(2,100),TSUB(100,8)
C
      COMMON/SPECES/COEF(2,7,100),S(100),EN(100,9),ENLN(100),HO(100)
     1,SUB(100,8),IUSE(100),TEMP(100,2),OMIT(8,100),ENSERT(8,100)
     2,JUSE(100),A(15,100)
      COMMON/MISC/ENN,SUMN,TT,LLMT(15),BO(15),TM,PP,R,ENNL,TRI,BASE(15)
     1,NAME(15,5),ANUM(15,5),TOP(15),ADD(15),VOL(100),GAMMA(100),NOCONV
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX
C
C#### TTY is not used, now used as "*"
C      DATA GAS/1HG/,END/4HENDD/, TTY/9/
      DATA GAS/1HG/,END/4HENDD/
C
      I2B = 0
        REWIND 3
        REWIND 4
      NC= 0
      IX= 0
      JJ = 0
      DO 300 I=1,15
  300 LMT(I) = 0
      NS = 1
      IF(NMIN.LT.1) GO TO 500
C
C     DETERMINE WHICH ELEMENTS WILL BE PUT INTO STORAGE.
C
      DO 510 I=1,NREAC
      IF(BASE(I).NE.0.) GO TO 510
      DO 512 J=1,5
      IF(ANUM(I,J).EQ.0.) GO TO 512
      DO 505 K= 1,NLM
      IF(NAME(I,J).EQ.LLMT(K)) GO TO 512
  505 CONTINUE
      JJ = JJ+1
      LMT(JJ) = NAME(I,J)
  512 CONTINUE
  510 CONTINUE
  500 MS=1
      NM = 0
C
C     SET NUMBER OF SPECIES ALLOWED, CLEAR A(I,J).
C
      DO 3 I=1,100
      DO 3  J=1,NLM
      KOUNT = I
      A(J,I) = 0.
    3 CONTINUE
      MAXNS = KOUNT-1
      IF(NEWMOL.NE.0) GO TO 77
C
C     BEGIN LOOP FOR READING SPECIES DATA FROM TAPE.
C
    7 IF(NEWMOL.NE.0) GO TO 77
C#### Changed to read one record and check for ENDD record
C      READ(4,10) (NAM(I),I=1,8),DATE(1,NS),DATE(2,NS),PHAZ,T1,T2,
C     1(MT(J),B(J),J=1,10),VOL(NS)
   10 FORMAT(8A4,6X,2A3,2X,A1,2F10.3/10(A2,F3.0),E15.8)
      READ(4,1010)
     A    (NAM(I),I=1,8),DATE(1,NS),DATE(2,NS),PHAZ,T1,T2
 1010 FORMAT(8A4,6X,2A3,2X,A1,2F10.3)
      IF(NAM(1).EQ.END)  GO TO 171
      READ(4,1011) (MT(J),B(J),J=1,10),VOL(NS)
 1011 FORMAT(10(A2,F3.0),E15.8)
      READ     (4,20) ((COEF(I,J,NS),J=1,7),I=1,2)
   20 FORMAT (5E15.8)
      GO TO 21
   77 READ(3,10)(NAM(I),I=1,8),DATE(1,NS),DATE(2,NS),PHAZ,T1,T2,
     1(MT(J),B(J),J=1,10),VOL(NS)
      IF(NAM(1).EQ.END) GO TO 171
      READ(3,20)((COEF(I,J,NS),J=1,7),I=1,2)
   21 IF(NEWMOL.GE.1) GO TO 811
      IF(NOMIT.EQ.0) GO TO 810
C
C     IS THIS SPECIES TO BE OMITTED.
C
      DO 805 I=1,NOMIT
      DO 804 J=1,8
      IF(OMIT(J,I).NE.NAM(J)) GO TO 805
  804 CONTINUE
      GO TO 7
  805 CONTINUE
  810 IF(NSERT.EQ.0) GO TO 811
C
C     IS THIS SPECIES TO BE INCLUDED.
C
      DO 400 I=1,NSERT
      DO 410 J=1,8
      IF(ENSERT(J,I).NE.NAM(J)) GO TO 400
  410 CONTINUE
      GO TO 811
  400 CONTINUE
      GO TO 7
  811 TEST = .FALSE.
      DO 820 K=1,10
      IF(B(K).EQ.0.) GO TO 825
      DO 168 I=1,NLM
      NI = I
      IF(LLMT(I).EQ.MT(K)) GO TO 822
  168 CONTINUE
      IF(JJ.EQ.0) GO TO 167
      DO 169 I=1,JJ
      NI = I
      IF(LMT(I).EQ.MT(K)) GO TO 822
  169 CONTINUE
C
C     THIS SPECIES WILL NOT BE INCLUDED.
C
  167 DO 819 J=1,NLM
      A(J,NS) = 0.
  819 CONTINUE
      GO TO 7
  822 IF(NMIN.LT.1) GO TO 824
      DO 520 L=1,JJ
      IF(LMT(L).NE.MT(K)) GO TO 824
  520 CONTINUE
      TEST = .TRUE.
  824 IF((NS+NM).LE.MAXNS) A(NI,NS) = B(K)
  820 CONTINUE
  825 IF(.NOT.TEST) GO TO 855
      NM = NM+1
      DO 600 I=1,NLM
  600 A(I,NS)=0.0
  855 IF((NS+NM).LE.MAXNS) GO TO 828
      I2B = I2B+1
      DO 826 I=1,8
  826 TOOBIG(I,I2B) = NAM(I)
      GO TO 7
  828 DO 829 I=1,8
      IF(.NOT.TEST) GO TO 827
C
C     THIS SPECIES WILL BE PUT INTO TEMPORARY STORAGE.
C
      TSUB(MS,I)=NAM(I)
      GO TO 829
C
C     THIS SPECIES WILL BE PUT INTO THE CURRENT SYSTEM.
C
  827 SUB(NS,I)=NAM(I)
  829 CONTINUE
      IF(TEST) TDATE(1,MS)=DATE(1,NS)
      IF(TEST) TDATE(2,MS)=DATE(2,NS)
      IF(NEWMOL.LT.1) WRITE(3,10)(NAM(I),I=1,8),DATE(1,NS),DATE(2,NS),
     1 PHAZ,T1,T2,(MT(J),B(J),J=1,10),VOL(NS)
      IF(NEWMOL.LT.1) WRITE(3,20) ((COEF(I,J,NS),J=1,7),I=1,2)
      JUSE(NS) = 0
      IUSE(NS) = 0
      IF(PHAZ.EQ.GAS) GO TO 170
C
C     CONDENSED SPECIES.
C
      IF(TEST) GO TO 170
      NC = NC+1
      TEMP(NC,1)= T1
      TEMP(NC,2)= T2
      IX= IX+1
      IF(NS.EQ.1) GO TO 245
      IF(IUSE(NS-1).EQ.0) GO TO 245
      DO 830 I=1,NLM
      IF(A(I,NS).NE.A(I,NS-1)) GO TO 245
  830 CONTINUE
      IX= IX-1
  245 IUSE(NS)=-IX
      JUSE(NS) = -IX
  170 IF(TEST) MS=MS+1
      IF(.NOT.TEST) NS=NS+1
      GO TO 7
C
C     END CARD HAS BEEN READ.
C
  171 CONTINUE
C#### TTY is not used, now used as "*"
C        WRITE(TTY,3333)
        WRITE(*,3333)
3333    FORMAT(' End of data reached in subroutine SEARCH')
        IF(NEWMOL.LT.1.AND.NMIN.GT.0) WRITE(3,99) END
   99 FORMAT(A4/2H  /2H  /)
        REWIND 3
        REWIND 4
      NS= NS-1
      MS=MS-1
      IF(PLUS.AND.TRI.GT.0..AND.NEWMOL.EQ.1) GO TO 173
      IF(NEWMOL.NE.0) RETURN
  173 WRITE(6,172)
  172 FORMAT(' SPECIES BEING CONSIDERED IN THIS SYSTEM'   )
      DO 174 I=1,NS,2
      I5=I+1
      IF(NS.LT.I5) I5=NS
  174 WRITE(6,176) (DATE(1,J),DATE(2,J),(SUB(J,K),K=1,8),J=I,I5)
  176 FORMAT(2(5X,2A3,2X,8A4))
      IF(I2B.GT.0)  GO TO 870
      IF(MS.LT.1) RETURN
      WRITE(6,540)
  540 FORMAT(/'  SPECIES WHICH HAVE BEEN PUT INTO TEMPORARY STORAGE')
      DO 550 I=1,MS,2
      I5 = I+1
      IF(MS.LT.I5) I5=MS
  550 WRITE(6,176) (TDATE(1,J),TDATE(2,J),(TSUB(J,K),K=1,8),J=I,I5)
      RETURN
  870 WRITE(6,871)  I2B
  871 FORMAT(' INSUFFICIENT STORAGE FOR FOLLOWING',I3,' SPECIES')
      WRITE(6,880) ((TOOBIG(K,J),K=1,8),J=1,I2B)
  880 FORMAT(2(3X,8A4))
      NS = 0
      RETURN
      END
      SUBROUTINE EQLBRM
C
C     ROUTINE TO CALCULATE EQUILIBIUM COMPOSITION AND PROPERTIES.
C
C     THE FOLLOWING DOUBLE PRECISION TYPE STATEMENTS ARE REQUIRED FOR
C     MACHINES WITH LESS THEN 60 BIT WORDS.
C
      DOUBLE PRECISION G,X,SUM,DELN
C
C     THE FOLLOWING DOUBLE PRECISION TYPE STATEMENTS ARE REQUIRED FOR
C     IBM 360 MACHINE ONLY.
C
C#### F is not used
C      DOUBLE PRECISION COEF,S,EN,ENLN,HO,F,AA
      DOUBLE PRECISION COEF,S,EN,ENLN,HO,AA
        DOUBLE PRECISION K1,K2,K4,AK1,AK2,AK4,TI,TKK
C
      LOGICAL CONVG,SINGC,ISING,MISS,SINGC2,NEW,GAS,SING3,SING4
      INTEGER SUB
C
      COMMON/POINTS/P(45),T(45),PPP(9),TTT(9),TOTN(9),RSIZE,TEM,GAS,NEW,
     1 GTOTN(9),FREE(9),OXFUG(9),FO2C(9),FO2H(9)
      COMMON/SPECES/COEF(2,7,100),S(100),EN(100,9),ENLN(100),HO(100)
     1,SUB(100,8),IUSE(100),TEMP(100,2),OMIT(8,100),ENSERT(8,100)
     2,JUSE(100),A(15,100)
      COMMON/MISC/ENN,SUMN,TT,LLMT(15),BO(15),TM,PP,R,ENNL,TRI,BASE(15)
     1,NAME(15,5),ANUM(15,5),TOP(15),ADD(15),VOL(100),GAMMA(100),NOCONV
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX
      COMMON /DOUBLE/ G(20,21), X(20)
C
      DIMENSION DELN(100),DELG(100)
      DIMENSION AK1(8),AK2(8),AK4(8)
C####
C     PUT THESE FOUR VARIABLES IN LOCAL COMMON BLOCK TO SAVE
C     THEIR PREVIOUS VALUES WHICH ARE USED WHEN ENTERING THE 
C     SUBROUTINE AGAIN.
C
      REAL  FCO
      REAL  FCO2
      REAL  FH2
      REAL  FH2O
      COMMON/FUGACITY/
     A  FCO,
     B  FCO2,
     C  FH2,
     D  FH2O
C
      DATA SMALNO/1.E-6/,SMNOL/-13.815511/,ITN/35/,E/2.7182818/
      DATA IO2/'O2'/, NCO/'CO'/,NCO2/'CO2'/,NH2/'H2'/,NH2O/'H2O'/
C
C ***
C       COEFFICIENTS FOR EQUILIBRIUM CONSTANTS FROM DATA IN JANAF TABLES
C       K1: C+O2=CO2; K2:C+.5O2=CO;K4: H2+.5O2=H2O
C       EQUATION FORM IS LOG10(K)=A1+A2*T+A3*T**2 ... AND FIT FROM 500-2300K
C ***
        DATA AK1/.114463D3,-.243442D2,.248074D1,-.114863D0,.584013D-3,
     1.154925D-3,-.555708D-5,.592449D-7/
        DATA AK2/.36845D2,-.689197D1,.711475D0,-.337899D-1,.223693D-3,
     1.4365D-4,-.162723D-5,.178614D-7/
      DATA AK4/.687667D2,-.153556D2,.158928D1,-.759149D-1,.535197D-3,
     1.96958D-4,-.366479D-5,.407351D-7/
C
      SING4 = .FALSE.
      SINGC = .FALSE.
      SING3 = .FALSE.
      SINGC2=.FALSE.
      MISS=.FALSE.
      JDELG = 0
      SIZEG = 0.
      NINT = 0
      NOCON2 = 0
      ENL = ENNL
      SIZE = 18.420681
      IF(RSIZE.NE.0.) SIZE= RSIZE
      ISING = .FALSE.
C
C     CHECK FOR CONDENSED SPECIES THAT BELONG TO PRESENT TEMP. RANGE.
C
      INC=0
      DO 310 J=1,NS
      IF(IUSE(J).EQ.0) GO TO 310
      IF(JUSE(J).GT.0) JUSE(J)=-JUSE(J)
      INC=INC+1
      IF(TT.GE.TEMP(INC,1).AND.TT.LE.TEMP(INC,2)) GO TO 310
      IF(JUSE(J).LT.0) JUSE(J)=-JUSE(J)
  310 CONTINUE
      CONVG = .FALSE.
C
C     CALL SUBROUTINE FCOEFF TO CALCULATE FUGACITY OF GASES.
C
      CALL FCOEFF
      ITNUMB = ITN
C
C     CALL SUBROUTINE CPHS TO CALCULATE THERMO PROPERTIES.
C
      CALL CPHS 
      TM = ALOG(PP/ENN)
C
C     REMOVE CONDENSED SPECIES WHICH DONOT FIT PRESENT TEMPERATURE.
C
      INC= 0
      DO 270 J=1,NS
      IF(IUSE(J).EQ.0) GO TO 270
      IF(JUSE(J).LT.0) GO TO 270
      JJ=J-1
      IF(JJ.LT.1) GO TO 450
      IF(IABS(JUSE(J)).EQ.IABS(JUSE(JJ))) GO TO 270
  450 JJ=J+1
      IF(JJ.GT.NS) GO TO 270
      IF(IABS(JUSE(J)).EQ.IABS(JUSE(JJ))) GO TO 270
      IF(IUSE(J).LT.0) GO TO 270
      INC=INC+1
  270 CONTINUE
      IF(INC.GT.0) GO TO 150
C
C     BEGIN ITERATION
C     CALL SUBROUTINE MATRIX TO SET UP THE ELEMENTS OF THE MATRIX.
C
   62 CALL MATRIX
      ITST = IMAT
C#### TTY (Unit 9) is not used, now used as "*"
C        WRITE(9,2700) ITNUMB
C       WRITE(*,2700) ITNUMB
C2700    FORMAT(' IT NO = ',I4)
C
C     CALL SUBROUTINE GAUSS TO SOLVE THE MATRIX.
C
      CALL GAUSS
      IF(ITST.EQ.IMAT) GO TO 85
C
C     SINGULAR MATRIX.
C
      IF(.NOT.NEW) WRITE(6,75)
C#### TTY (Unit 9) is not used, now used as "*"
C      WRITE(9,74) TT
      WRITE(*,74) TT
   74 FORMAT(31H SINGULAR MATRIX AT TEMPERATURE,F10.0)
      NEW = .TRUE.
      IF(SING3.AND..NOT.SING4) GO TO 150
      IF(SINGC2) GO TO 873
      IF(SINGC) GO TO 875
      DO 970 JJ = 1, NS
      IF(IUSE(JJ).NE.0)  GO TO 970
      IF(EN(JJ,NPT).NE.0.)  GO TO 970
      EN(JJ,NPT) = SMALNO
      ENLN(JJ) = SMNOL
  970 CONTINUE
      IF(ISING)  GO TO 870
      ISING = .TRUE.
C#### TTY (Unit 9) is not used, now used as "*"
C     WRITE (9,776)
      WRITE (*,776)
  776 FORMAT (' RESTART')
      GO TO 62
C
C     TEST FOR SINGULARITY TO CONDENSED SPECIES.
C
  870 NCOND = IQ1-NLM-1
      IF(NCOND.LT.2.OR.SIZEG.EQ.0.) GO TO 874
      DO 872  J=1,NS
      IF(IUSE(J).LE.0)  GO TO 872
      IF(J.EQ.JDELG)  GO TO 872
      DO 671  I=1,NLM
      IF(A(I,J).EQ.A(I,JDELG))  GO TO 671
      IF(A(I,J).EQ.0..OR.A(I,JDELG).EQ.0.)  GO TO 872
  671 CONTINUE
      SINGC = .TRUE.
      IQ1 = IQ1-1
      EN(J,NPT) = 0.
      IUSE(J) = -IUSE(J)
      IF(JUSE(J).LT.0) JUSE(J) = -JUSE(J)
  872 CONTINUE
      IF(SINGC)  GO TO 40
  874 SINGC=.TRUE.
      DO 20 J = 1,NS
      IF(IUSE(J).LE.0) GO TO 20
      SUM = 0.
      DO 21 I = 1,NLM
   21 SUM = SUM + A(I,J) * X(I)
      DELG(J) = HO(J) - S(J) - SUM
      IF(PP.GT.1.) DELG(J)=HO(J)-S(J)+((VOL(J)*(PP-1.0))/(R*TT))-SUM
   20 CONTINUE
      DO 410 I=1,NS
      IF(IUSE(I).LE.0) GO TO 410
C#### TTY (Unit 9) is not used, now used as "*"
C      WRITE(9,420) (SUB(I,J),J=1,8),DELG(I),HO(I),S(I)
      WRITE(*,420) (SUB(I,J),J=1,8),DELG(I),HO(I),S(I)
  420 FORMAT(1X,8A4,E15.8,2F15.4)
  410 CONTINUE
C
C     REMOVE CONDENSED SPECIES WITH GREATEST DELG VALUE.
C
      SING3 = .TRUE.
  880 ASIZE = -1.E+36
      DO 30 J = 1,NS
      IF(IUSE(J).LE.0) GO TO 30
      IF(DELG(J).LT.0.) GO TO 31
      GO TO 30
   31 IF(DELG(J).LT.ASIZE) GO TO 30
      KOT = J
      ASIZE = DELG(J)
   30 CONTINUE
      IF(ASIZE.EQ.-1.E+36) GO TO 150
      IF(JUSE(KOT).LT.0) JUSE(KOT) = -JUSE(KOT)
      IUSE(KOT) = -IUSE(KOT)
      EN(KOT,NPT) = 0.
      IQ1 = IQ1 - 1
      GO TO 40
C
C     THIRD AND LAST ATTEMPT TO CORRECT SINGULAR MATRIX.
C
  875 SINGC2=.TRUE.
      MISS= .FALSE.
      GO TO 150
C
C     OBTAIN CORRECTIONS TO THE ESTIMATES
C
   85 ITNUMB= ITNUMB-1
      KK = NLM + 1
      SUM = X(IQ1)
      DO 101 J=1,NS
      IF (IUSE(J)) 101,98,100
   98 DELN(J)=-HO(J)+S(J)-ENLN(J)-TM-GAMMA(J)+SUM
      DO 99 K=1,NLM
      DELN(J)= DELN(J)+A(K,J)*X(K)
   99 CONTINUE
      GO TO 101
  100 DELN(J) = X(KK)
      KK = KK + 1
  101 CONTINUE
C
C     CALCULATE CONTROL FACTOR,AMBDA
C
      AMBDA= 1.
      AMBDA1= 1.
      IF(SUM.LT.0.)  SUM=-SUM
      DO 917 J=1,NS
      IF (IUSE(J).NE.0) GO TO 917
      IF((EN(J,NPT).GT.0.).AND.DELN(J).GT.SUM) SUM = DELN(J)
      IF((EN(J,NPT).NE.0.) .OR. DELN(J).LE.0.) GO TO 917
      IF(DABS(DELN(J)-X(IQ1)).LT..9D-99)  GO TO 150
      SUM1 = (-9.212-ENLN(J)+ ENL)/(DELN(J)-X(IQ1))
      IF(SUM1.LT.0.)  SUM1=-SUM1
      IF (SUM1.LT.AMBDA1) AMBDA1 = SUM1
  917 CONTINUE
      IF(SUM.GT.2.)AMBDA=2./SUM
      IF (AMBDA1.LT.AMBDA) AMBDA = AMBDA1
C
C     COMPUTE TOTAL FREE ENERGY.
C
C
C     APPLY CORRECTIONS TO ESTIMATES
      SUM = 0.
      DO 113 J=1,NS
      IF (IUSE(J)) 113,112,114
  112 ENLN(J)=ENLN(J)+AMBDA*DELN(J)
      EN(J,NPT) = 0.
      IF((ENLN(J)-ENL+SIZE).LE.0..AND.SUB(J,1).NE.IO2) GO TO 113
      EN(J,NPT) = E**ENLN(J)
      SUM = SUM+EN(J,NPT)
      GO TO 113
  114 EN(J,NPT) = EN(J,NPT) + AMBDA * DELN(J)
  113 CONTINUE
      SUMN = SUM
      ENL = ENL+AMBDA*X(IQ1)
      ENN = E**ENL
      TM = ALOG(PP/ENN)
C
C     TEST FOR CONVERGENCE
C
      IF(ITNUMB.EQ.0) GO TO 14
      IF(AMBDA.LT.1.) GO TO 62
      SUM = X(IQ1)
      IF (SUM.LT.0.) SUM = -SUM
      IF(SUM.GT.0.5E-5) GO TO 62
      DO 130 J=1,NS
      IF (IUSE(J).LT.0) GO TO 130
      AA= DELN(J)/SUMN
      IF(AA.LT.0.)  AA=-AA
      IF (IUSE(J).EQ.0) AA = AA*EN(J,NPT)
      IF(AA.GT.0.5E-5) GO TO 62
  130 CONTINUE
      CONVG= .TRUE.
      GO TO 160
C
C     NO CONVERGENCE, PRINT ERROR MESSAGE.
C
C#### TTY (Unit 9) is not used, now used as "*"
C   14 IF(.NOT.NEW) WRITE(9,75)
   14 IF(.NOT.NEW) WRITE(*,75)
   75 FORMAT(' ')
C#### TTY (Unit 9) is not used, now used as "*"
C      WRITE(9,973) ITN,NPT
      WRITE(*,973) ITN,NPT
  973 FORMAT(' ',I3,' ITERATIONS DID NOT SATISFY CONVERGENCE REQUIREMENT
     1S FOR THE POINT ',I5)
      NEW = .TRUE.
      NOCON2 = NOCON2 + 1
      IF(NOCON2.GT.2) GO TO 873
      IF(NOCON2.EQ.2) GO TO 880
C
C     RESET INITIAL ESTIMATES.
C
  150 IF(NINT.GT.6) GO TO 873
      SING4 = .TRUE.
      CALL INIT
      ENL = ENNL
      ITNUMB = ITN
      CONVG = .FALSE.
      NINT = NINT + 1
      GO TO 62
C
C     CONVERGENCE TESTS ARE SATISFIED, TEST CONDENSED SPECIES.
C
  160 IF(NC.EQ.0) GO TO 143
C
C     DETERMINE IF CONDENSED SPECIES SHOULD BE REMOVED.
C
      DO 146  J=1,NS
      IF(EN(J,NPT).GE.0.) GO TO 146
      IF(NINT.GT.15) GO TO 873
      IQ1= IQ1-1
      EN(J,NPT) = 0.
      NINT = NINT + 1
      IF(JUSE(J).LT.0) JUSE(J) = -JUSE(J)
      IF(IUSE(J).GT.0) IUSE(J) = -IUSE(J)
      GO TO 40
  146 CONTINUE
      SIZEG = 0.
      INC = 0
      DO 170 J = 1,NS
      IF(IUSE(J).EQ.0)  GO TO 170
      INC = INC + 1
      IF(IUSE(J).GT.0) GO TO 169
      DO 35 L=1,10
      IF(JUSE(J).GE.0) GO TO 169
      JJ=J-L
      IF(JJ.LT.1) GO TO 45
      IF(JUSE(J).EQ.-JUSE(JJ).AND.EN(JJ,NPT).GT.0.) GO TO 1165
   45 JJ=J+L
      IF(JJ.GT.NS) GO TO 35
      IF(JUSE(J).EQ.-JUSE(JJ).AND.EN(JJ,NPT).GT.0.) GO TO 1165
   35 CONTINUE
      GO TO 153
C
C      WRONG PHASE INCLUDED, SWITCH EN.
C
 1165 EN(J,NPT) = EN(JJ,NPT)
      IF(IUSE(J).LT.0) IUSE(J)=-IUSE(J)
      IF(IUSE(JJ).GT.0) IUSE(JJ)=-IUSE(JJ)
      DELG(J) = DELG(JJ)
      EN(JJ,NPT) = 0.
      GO TO 40
C
C     FIND THE SPECIE WITH THE MINIMUM NEGATIVE DELTA G.
C
  153 SUM = 0.
      DO 167 I=1,NLM
      SUM = SUM + A(I,J)*X(I)
  167 CONTINUE
      DELG(J) = HO(J) - S(J)-SUM
      IF(PP.GT.1.) DELG(J)=HO(J)-S(J)+((VOL(J)*(PP-1.0))/(R*TT))-SUM
      IF(.NOT.MISS) GO TO 169
      IF(DELG(J).GE.SIZEG.OR.DELG(J).GE.0.) GO TO 169
      SIZEG = DELG(J)
      JDELG = J
  169 IF(INC.EQ.NC) GO TO 1160
  170 CONTINUE
 1160 IF(MISS) GO TO 558
C
C     FIRST RUN, FIND SPECIE WITH NEXT LOWEST DELTA G VALUE.
C
      JKL=0
      LKJ=0
      KJ=0
  556 SIZEG=0.
      DO 555 J=1,NS
      IF(JUSE(J).GE.0) GO TO 555
      IF(IUSE(J).GE.0) GO TO 555
      IF(DELG(J).LT.0) LKJ=LKJ+1
      IF(DELG(J).GE.SIZEG.OR.DELG(J).GE.0.) GO TO 555
      KJ=J
      SIZEG=DELG(J)
  555 CONTINUE
      IF(JKL.GT.0) GO TO 557
      IF(LKJ.LT.2) GO TO 557
      DELG(KJ)=0.0
      JKL=JKL+1
      GO TO 556
  557 IF(KJ.EQ.0) GO TO 520
      SIZEG=DELG(KJ)
      JDELG=KJ
  520 MISS=.TRUE.
  558 IF(SIZEG.EQ.0.) GO TO 143
      JT = JDELG
      IQ1=IQ1+1
  166 IUSE(JT) = -IUSE(JT)
   40 CONVG = .FALSE.
  143 ITNUMB = ITN
      IF(.NOT.CONVG) GO TO 62
C
C     STORE PRESSURE, TEMP, AND N VALUES FOR THIS POINT.
C
      TTT(NPT)=TT
      ENNL = ENL
      PPP(NPT) = PP
      TOTN(NPT) = 0.
      GTOTN(NPT) = 0.
      DO 183 J=1,NS
      IF(IUSE(J).LT.0)  GO TO 183
      IF(IUSE(J).EQ.0) GTOTN(NPT) = GTOTN(NPT) + EN(J,NPT)
      TOTN(NPT) = TOTN(NPT) + EN(J,NPT)
  183 CONTINUE
C
C     COMPUTE TOTAL FREE ENERGY.
C
      FREE(NPT) = 0.
      DO 700 J = 1,NS
      IF(IUSE(J).LT.0) GO TO 700
      TEN = EN(J,NPT)/TOTN(NPT)
C       STORE VALUES OF CO, CO2, H2, AND H2O FUGACITIES FOR FO2 CALC
        IF(SUB(J,1).EQ.NCO)  FCO=TEN*EXP(GAMMA(J))
        IF(SUB(J,1).EQ.NCO2) FCO2=TEN*EXP(GAMMA(J))
        IF(SUB(J,1).EQ.NH2)  FH2= TEN*EXP(GAMMA(J))
        IF(SUB(J,1).EQ.NH2O) FH2O=TEN*EXP(GAMMA(J))
      IF(TEN.LT.5.E-6) GO TO 700
      IF(IUSE(J).EQ.0) GO TO 701
      TVOL = 0.
      IF(PP.GT.1.) TVOL=(VOL(J)*(PP-1.0))/(R*TT)
      FREE(NPT) = FREE(NPT)+((HO(J)-S(J)+TVOL)*TEN)
      GO TO 700
  701 FREE(NPT) = FREE(NPT) + ((HO(J)-S(J)+ENLN(J)+TM+GAMMA(J))*TEN)
  700 CONTINUE
      DO 705 J=1,NS
      IF(SUB(J,1).NE.IO2) GO TO 705
      OXFUG(NPT) = 0.434294*GAMMA(J)
  705 CONTINUE
C       CALC LOG10 FO2 FROM CO/CO2 AND H2/H2O
C *** CALC LOGS OF EQUILIBRIUM CONSTANTS
        TKK=0.01*TT
        K1=0.0
        K2=0.0
        K4=0.0
        DO 10 I=1,8
        TI=TKK**(I-1)
        K1=K1+AK1(I)*TI
        K2=K2+AK2(I)*TI
        K4=K4+AK4(I)*TI
10      CONTINUE
        XLKC=K1-K2
        XLKH=K4
C       FCO2,FCO,FH2O,AND FH2 ARE LOG10 VALUES
        IF(FCO.LT.1.E-20) FCO=1.E-20
        IF(FCO2.LT.1.E-20) FCO2=1.E-20
        FO2C(NPT)=2.0*(ALOG10(FCO2/FCO)-XLKC)
        IF(FH2.LT.1.E-20) FH2=1.E-20
        IF(FH2O.LT.1.E-20) FH2O=1.E-20
        FO2H(NPT)=2.0*(ALOG10(FH2O/FH2)-XLKH)
      RETURN
C
C     ERROR, SET TT=0
C
  873 TT=0.
      NPT = NPT-1
      NOCONV = NOCONV+1
      RETURN
      END
C
      SUBROUTINE CPHS
C
C     CALCULATES THERMODYNAMIC PROPERTIES FOR INDIVIDUAL SPECIES
C
C     THE FOLLOWING DOUBLE PRECISION TYPE STATEMENTS ARE REQUIRED FOR
C     IBM 360 MACHINE ONLY.
C
      DOUBLE PRECISION COEF,S,EN,ENLN,HO
C
      DIMENSION F(100)
C
      COMMON/SPECES/COEF(2,7,100),S(100),EN(100,9),ENLN(100),HO(100)
     1,SUB(100,8),IUSE(100),TEMP(100,2),OMIT(8,100),ENSERT(8,100)
     2,JUSE(100),A(15,100)
      COMMON/MISC/ENN,SUMN,TT,LLMT(15),BO(15),TM,PP,R,ENNL,TRI,BASE(15)
     1,NAME(15,5),ANUM(15,5),TOP(15),ADD(15),VOL(100),GAMMA(100),NOCONV
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX
C
      TLN = ALOG(TT)
      DO 90 J=1,NS
      K = 1
      IF(TT.LE.1000.) K=2
      KK = K
      IF(COEF(K,1,J).NE.0.) GO TO 97
      IF(IUSE(J).LT.0) GO TO 90
C
C     IF COEFFICIENTS ARE ZERO, USE OTHER TEMPERATURE INTERVAL.
C
      IF(KK.EQ.2) K=1
      IF (KK.EQ.1) K = 2
   97 S(J) =  ((((COEF(K,5,J)/4.)*TT+ COEF(K,4,J)/3.)*TT+ COEF(K,3,J)/2.
     1)* TT+COEF(K,2,J))*TT+ COEF(K,1,J)*TLN + COEF(K,7,J)
      HO(J) = ((((COEF(K,5,J)/5.)*TT+ COEF(K,4,J)/4.)*TT+ COEF(K,3,J)/3.
     1  ) *TT+ COEF(K,2,J)/2.)*TT+ COEF(K,1,J) + COEF(K,6,J)/TT
      IF(JUSE(J).GE.0) GO TO 90
      F(J) = HO(J) - S(J)
      IF(PP.GT.1.) F(J) =HO(J)-S(J)+((VOL(J)*(PP-1.))/(R*TT))
   90 CONTINUE
C
C     FIND PHASE OF CONDENSED SPECIE WITH MINIMUM DELTA F VALUE.
C
  200 DO 30 I=1,NS
      IF(JUSE(I).GE.0) GO TO 30
      DO 35 L=1,10
      IF((I-L).LE.0) GO TO 40
      IF(JUSE(I).NE.JUSE(I-L)) GO TO 40
      IF(F(I).GT.F(I-L)) GO TO 40
      JUSE(I-L) = -JUSE(I-L)
   40 IF((I+L).GT.NS) GO TO 30
      IF(JUSE(I).NE.JUSE(I+L)) GO TO 35
      IF(F(I).GT.F(I+L)) GO TO 35
      JUSE(I+L)=-JUSE(I+L)
   35 CONTINUE
   30 CONTINUE
      RETURN
      END
C
      SUBROUTINE MATRIX
C
C     SETS UP A MATRIX FOR LAGRANGIAN MULTIPLIER SOLUTION.
C
C     THE FOLLOWING DOUBLE PRECISION TYPE STATEMENTS ARE REQUIRED FOR
C     MACHINES WITH LESS THEN 60 BIT WORDS.
C
      DOUBLE PRECISION G,X
C
C     THE FOLLOWING DOUBLE PRECISION TYPE STATEMENTS ARE REQUIRED FOR
C     IBM 360 MACHINE ONLY.
C
      DOUBLE PRECISION COEF,S,EN,ENLN,HO,F
C
      COMMON/SPECES/COEF(2,7,100),S(100),EN(100,9),ENLN(100),HO(100)
     1,SUB(100,8),IUSE(100),TEMP(100,2),OMIT(8,100),ENSERT(8,100)
     2,JUSE(100),A(15,100)
      COMMON/MISC/ENN,SUMN,TT,LLMT(15),BO(15),TM,PP,R,ENNL,TRI,BASE(15)
     1,NAME(15,5),ANUM(15,5),TOP(15),ADD(15),VOL(100),GAMMA(100),NOCONV
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX
      COMMON /DOUBLE/ G(20,21), X(20)
C
      IQ2 = IQ1 + 1
      KMAT = IQ2
      IMAT = IQ1
C
C     CLEAR MATRIX STORAGES TO ZERO
C
      DO 211 I=1,IMAT
      DO 211 K=1,KMAT
      G(I,K) = 0.0
  211 CONTINUE
      G(IQ2,IQ1) = 0.
C
C     BEGIN SET UP OF ITERATION MATRIX
C
      KK = NLM
      DO 65  J=1,NS
      IF(IUSE(J).LT.0)  GO TO 65
      IF(IUSE(J).GT.0)  GO TO 70
      F = (HO(J)-S(J)+ENLN(J)+TM+GAMMA(J))*EN(J,NPT)
C
C     CALCULATE THE ELEMENTS G(I,K)
C
      DO 55 I=1,NLM
      IF (A(I,J) .EQ. 0.) GO TO 55
      TERM= A(I,J)*EN(J,NPT)
      DO 15 K = I,NLM
      G(I,K)= G(I,K) + A(K,J)*TERM
   15 CONTINUE
C
      G(I,IQ1)=G(I,IQ1)+TERM
      G(I,IQ2) = G(I,IQ2) + A(I,J) * F
   55 CONTINUE
      G(IQ1,IQ2) = G(IQ1,IQ2) + F
      GO TO 65
C
C     CONDENSED SPECIES.
C
   70 KK = KK + 1
      DO 75 I = 1,NLM
      G(I,KK) = A(I,J)
      G(I,KMAT) = G(I,KMAT) - A(I,J)*EN(J,NPT)
   75 CONTINUE
      G(KK,KMAT) = HO(J)-S(J)
      IF(PP.GT.1.) G(KK,KMAT)=HO(J)-S(J)+((VOL(J)*(PP-1.0))/(R*TT))
   65 CONTINUE
      G(IQ1,IQ1) = SUMN - ENN
C
C     REFLECT SYMMETRIC PORTIONS OF THE MATRIX
C
      DO 102 I=1,IQ1
      DO 102 J=I,IQ1
      G(J,I)=G(I,J)
  102 CONTINUE
C
C     COMPLETE THE RIGHT HAND SIDE
C
      DO 145 I=1,NLM
      G(I,KMAT) = G(I,KMAT)+(BO(I)-G(I,IQ1))
  145 CONTINUE
      G(IQ1,KMAT) = G(IQ1,KMAT)+ENN-SUMN
      RETURN
      END
C
C
      SUBROUTINE GAUSS
C
C     SOLVE ANY LINEAR SET OF UP TO 20 EQUATIONS
C     NUMBER OF EQUATIONS = IMAT
C
C     THE FOLLOWING DOUBLE PRECISION TYPE STATEMENTS ARE REQUIRED FOR
C     MACHINES WITH LESS THEN 60 BIT WORDS.
C
      DOUBLE PRECISION G,X,COEFX,SUM,Z
C
      COMMON/DOUBLE/G(20,21),X(20)
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX
C
      DIMENSION COEFX(20)
      DATA BIGNO/1.E+38/
C
C     BEGIN ELIMINATION OF NNTH VARIABLE
C
      IUSE1 = IMAT+1
      DO 45 NN=1,IMAT
      IF(NN.NE.IMAT) GO TO 8
      IF(G(NN,NN))31,23,31
C
C     SEARCH FOR MAXIMUM COEFFICIENT IN EACH ROW
C
    8 DO 18  I=NN,IMAT
      COEFX(I) = BIGNO
      IF(G(I,NN).EQ.0.)  GO TO 18
      COEFX(I) = 0.
      DO 10 J=NN,IUSE1
      SUM = G(I,J)
      IF(SUM.LT.0.)  SUM=-SUM
      IF(J.NE.NN)  GO TO 9
      Z = SUM
      GO TO 10
    9 IF(SUM.GT.COEFX(I))  COEFX(I)=SUM
   10 CONTINUE
      COEFX(I) = COEFX(I)/Z
   18 CONTINUE
C
C     LOCATE ROW WITH SMALLEST MAXIMUM COEFFICIENT
C
      TEMP = BIGNO
      I=0
      DO 22 J=NN,IMAT
      IF(COEFX(J).GE.TEMP) GO TO 22
      TEMP = COEFX(J)
      I=J
   22 CONTINUE
      IF(I.EQ.0) GO TO 23
C
C     INDEX I LOCATES EQUATION TO BE USED FOR ELIMINATING THE NTH
C     VARIABLE FROM THE REMAINING EQUATIONS
C
C      INTERCHANGE EQUATIONS I AND NN
C
      IF(NN.EQ.I) GO TO 31
      DO 30 J=NN,IUSE1
      Z=G(I,J)
      G(I,J)=G(NN,J)
      G(NN,J)=Z
   30 CONTINUE
C
C     DIVIDE NTH ROW BY NTH DIAGONAL ELEMENT AND ELIMINATE THE NTH
C     VARIABLE FROM THE REMAINING EQUATIONS
C
   31 K = NN + 1
      DO 36 J = K, IUSE1
      IF(G(NN,NN).EQ.0.)  GO TO 23
      G(NN,J) = G(NN,J)/G(NN,NN)
   36 CONTINUE
      IF(K.EQ.IUSE1) GO TO 45
      DO 44 I=K,IMAT
      DO 44 J=K,IUSE1
      G(I,J) = G(I,J) - G(I,NN)*G(NN,J)
   44 CONTINUE
   45 CONTINUE
C
C     BACKSOLVE FOR THE VARIABLES
C
      K = IMAT
   47 J = K + 1
      X(K)=0.0
      SUM = 0.0
      IF(IMAT.LT.J) GO TO 51
      DO 50 I=J,IMAT
      SUM = SUM+G(K,I)*X(I)
   50 CONTINUE
   51 X(K) = G(K,IUSE1) - SUM
      K = K - 1
      IF (K) 47,151,47
   23 IMAT = IMAT-1
  151 RETURN
      END
      SUBROUTINE THERMP
C
C     SETS CONDITIONS FOR SYSTEM BEING CONSIDERED.
C
C#### F is not used
C      DOUBLE PRECISION COEF,S,EN,ENLN,HO,F
      DOUBLE PRECISION COEF,S,EN,ENLN,HO
      LOGICAL ALTER,SYSTM,PLUS,TEM,NEW
C
      COMMON/POINTS/P(45),T(45),PPP(9),TTT(9),TOTN(9),RSIZE,TEM,GAS,NEW,
     1 GTOTN(9),FREE(9),OXFUG(9),FO2C(9),FO2H(9)
      COMMON/SPECES/COEF(2,7,100),S(100),EN(100,9),ENLN(100),HO(100)
     1,SUB(100,8),IUSE(100),TEMP(100,2),OMIT(8,100),ENSERT(8,100)
     2,JUSE(100),A(15,100)
      COMMON/MISC/ENN,SUMN,TT,LLMT(15),BO(15),TM,PP,R,ENNL,TRI,BASE(15)
     1,NAME(15,5),ANUM(15,5),TOP(15),ADD(15),VOL(100),GAMMA(100),NOCONV
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX
C
      NEW = .TRUE.
      NMIX = 0
C#### IRAT is not used (code commented out at end of this routine)
C      IF(PLUS.AND.TRI.NE.0.) IRAT=(BASE(1)-TOP(1))/TRI
C      IF(SYSTM) IRAT=(BASE(1)-TOP(1))/TRI
      NPT = 1
      NEWMOL = 1
      NRAT = 0
      IF(PLUS.AND.TRI.NE.0.) GO TO 95
      GO TO 96
C
C     CALL SUBROUTINE MOLEWT TO CALCULATE NEW MOLE FRACTIONS.
C
C#### Changed following code to get MOLEWT to work right
   95 CONTINUE
      ALTER = .FALSE.
      CALL MOLEWT
      NEWMOL = 2
C
C     SET ASSIGNED P
C
   96 IP = 0
      IF(TEM) GO TO 97
      IP = 0
  903 IP = IP + 1
      PP = P(IP)
C
C     SET ASSIGNED T
C
   97 DO 902 IT=1,NT
      IF(TEM) PP=P(IT)
      IF(TEM) IP=IP+1
      TT = T(IT)
C
C     CALL SUBROUTINE EQLBRM TO CALCULATE EQUILIBRIUM COMPOSITIONS.
C
      CALL EQLBRM
      IF(TEM.AND.NOCONV.EQ.NT) GO TO 900
      IF(NOCONV.GT.2) GO TO 860
      IF(.NOT.TEM.AND.NOCONV.EQ.(NT*NP)) GO TO 900
      IF(IP.EQ.NP.AND.IT.EQ.NT) GO TO 860
      ISV = NPT
      IF(NPT.NE.6) GO TO 870
  860 WRITE(6,5)
    5 FORMAT(' '//)
      NEW = .FALSE.
      IF(NPT.EQ.0) GO TO 900
C
C     CALL SUBROUTINE OUT1 TO PRINT FINAL RESULTS.
C
      CALL OUT1
      IF(TT.LT.1.E-6) GO TO 900
      NPT = 0
  870 NPT = NPT + 1
      IF(TT.GT.1.E-3) GO TO 800
      IF(NPT.EQ.0) NPT=1
      CALL INIT
      GO TO 902
  800 IF(.NOT.TEM.AND.IT.NE.NT) GO TO 302
      IF(IT.EQ.NT) GO TO 902
      IF(TEM.AND.ABS(P(IP)-P(IP+1)).GT.1.E-10) CALL INIT
      IF(TEM.AND.ABS(P(IP)-P(IP+1)).GT.1.E-10) GO TO 902
C
C     USE COMPOSITIONS FROM PREVIOUS POINT.
C
  302 DO 300 J=1,NS
      EN(J,NPT) = EN(J,ISV)
  300 CONTINUE
  902 CONTINUE
      IF(TEM) GO TO 900
      CALL INIT
      IF(IP.LT.NP) GO TO 903
C
C     CHECK TO DETERMINE IF MOLE FRACTIONS ARE TO BE COMPUTED.
C
  900 NOCONV = 0
C       SKIP THE FOLLOWING AT THIS TIME
        GO TO 30  
C#### Commented out following lines, because go to 30 went to end of do-loop
C#### which throws it back into the loop
C      IF(SYSTM) GO TO 766
C      IF(ALTER) GO TO 700
C      IF(PLUS.AND.TRI.EQ.0.) GO TO 700
C  766 IF(IRAT.GT.NMIX) GO TO 95
C      IF(IRAT.EQ.NMIX) RETURN
C  700 CONTINUE 
C      DO 30 I=1,NREAC
C      IF(BASE(I).LT.1.E-6) GO TO 95
C      IF(ADD(I).LT.0.) GO TO 755
C      IF((BASE(I)+ADD(I)*.5).LT.TOP(I)) GO TO 95
C      GO TO 30
C  755 IF((BASE(I)+ADD(I)*.5).GT.TOP(I)) GO TO 95
   30 CONTINUE
      RETURN
      END
C
      SUBROUTINE MOLEWT
C
C     COMPUTES MOLES OF EACH REACTANT.
C
      LOGICAL ALTER,SYSTM,PLUS
      DIMENSION BB(15),TO(15),AA(15)
C
      COMMON/MISC/ENN,SUMN,TT,LLMT(15),BO(15),TM,PP,R,ENNL,TRI,BASE(15)
     1,NAME(15,5),ANUM(15,5),TOP(15),ADD(15),VOL(100),GAMMA(100),NOCONV
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX 
      INTEGER NEWMOL_SAVE
C
C     DETERMINE WHICH SYSTEM IS BEING CONSIDERED.
C
      NMIN = 0
      IF(ALTER) GO TO 90
      IF(PLUS.AND.TRI.NE.0..AND.NEWMOL.LT.2) GO TO 66
      IF(PLUS.AND.TRI.NE.0..AND.NEWMOL.GE.2) GO TO 20
      IF(PLUS) GO TO 45
      IF(SYSTM.AND.NEWMOL.GE.2) GO TO 20
      NEW = 1
      LL = 1
      GO TO 66
C
C     SHOULD A NEW REACTANT BE ADDED OR SHOULD THE CONSTANT REACTANT
C     BE INCREMENTED.
C
   20 CONTINUE
      DO 35 I=1,NREAC
      IF(BASE(I).EQ.0.) GO TO 35
      IF(ADD(I).LT.0.) GO TO 750
      IF((BASE(I)+ADD(I)*.5).LT.TOP(I)) GO TO 90
      GO TO 35
  750 IF((BASE(I)+ADD(I)*.5).GT.TOP(I)) GO TO 90
   35 CONTINUE
      IF(SYSTM) GO TO 55
      IF(PLUS) GO TO 850
C
C     INCREMENT MOLES BETWEEN INCLUDED REACTANTS.
C
   90 CONTINUE
      DO 210 I=1,NREAC
        RTEST=ABS(ADD(I))
        ZEROT= 1.0E-4
        IF(RTEST.LT.ZEROT) GO TO 210
      IF(BASE(I).EQ.0.) GO TO 210
      IF(BASE(I).EQ.TOP(I)) GO TO 210
      BASE(I) = BASE(I)+ADD(I)
      IF(ADD(I).LT.0.) GO TO 220
      IF(BASE(I).GT.TOP(I)) BASE(I) = TOP(I)
      GO TO 210
  220 IF(BASE(I).LT.TOP(I)) BASE(I)=TOP(I)
  210 CONTINUE
      IF(NPT.EQ.0) NPT=1
      CALL REACT
      GO TO 25
C
C     ADD AND REMOVE ONE REACTANT.
C
   55 CONTINUE
      IF(NRAT.GE.1) GO TO 850
      IF(NEW.NE.2) GO TO 180
      LL = LL+1
      ADD(2) = ADD(1)
      BASE(1) = 0.
      TOP(1) = 0.
      ADD(3) = 0.
      NRAT = NRAT+1
  180 DO 70 I=LL,NREAC
      IF(LL.GE.2) GO TO 147
      IF(BASE(I).EQ.0.) GO TO 150
  147 IF(ADD(I).GT.0.) GO TO 160
      IF(ADD(I).EQ.0.) GO TO 150
      BASE(I) = BB(1)
      TOP(I) = TO(1)
      GO TO 70
  150 BASE(I) = BB(2)
      TOP(I) = TO(2)
      ADD(I) = AA(2)
      GO TO 70
  160 BASE(I) = 0.0
      TOP(I) = 0.0
      ADD(I) = 0.0
   70 CONTINUE
      NEW = NEW+1
      GO TO 999
C
C     STORE BASE(I), TOP(I), AND ADD(I) VALUES.
C
   66 CONTINUE
      DO 80 I=1,NREAC
      BB(I) = BASE(I)
      TO(I) = TOP(I)
      AA(I) = ADD(I)
   80 CONTINUE
      IF(SYSTM) GO TO 90
C
C     ADD ALL THE REACTANTS TO THE SYSTEM, OR INCREMENT EACH REACTANT.
C
  850 CONTINUE
      NMIX = NMIX+1
      TIM = NMIX
      IF(PLUS) NRAT = 1
      BASE(1) = BB(1)-(TRI*TIM)
      TOP(1) = TO(1)
      ADD(1) = -TRI
      BASE(2) = BB(2)
      TOP(2) = TO(2)-(TRI*TIM)
      ADD(2) = TRI
      BASE(3) = AA(3)+(TRI*TIM)
      TOP(3) = BASE(3)
      ADD(3) = 0.
      GO TO 999
C
C     ADD A NEW REACTANT TO THE SYSTEM.
C
   45 CONTINUE
      DO 450 I=1,NREAC
      IF(BASE(I).NE.0.) GO TO 450
      BASE(I) = TOP(I)
      ADD(I) = 0.0
      GO TO 999
  450 CONTINUE
C
C     CALL SUBROUTINE REACT TO COMPUTE NEW MOLE FRACTION.
C
  999 CONTINUE
      CALL REACT
      IF(NPT.EQ.0) NPT=1
      IF(PLUS.AND.TRI.NE.0..AND.NEWMOL.GE.2) GO TO 25
      IF(SYSTM.AND.NMIX.GT.1) GO TO 25
C
C     CALL SUBROUTINE SEARCH TO ADD OR SUBTRACT SPECIES.
C
C#### Following code changed for error in first call to SEARCH
      NEWMOL_SAVE = NEWMOL
      NEWMOL = 0
      CALL SEARCH
      NEWMOL = NEWMOL_SAVE
C
C     RESET THE INITIAL ESTIMATES.
C
   25 CONTINUE
      NPT=1
C#### TTY (Unit 9) is not used, now used as "*"
C        WRITE(9,2000)
      WRITE(*,2000) ALTER,SYSTM,PLUS
2000  FORMAT(' ALTER, SYSTEM, PLUS = ',3(L1,2X))
      CALL INIT
      RETURN
      END
C
      SUBROUTINE INIT
C
C#### F is not used
C      DOUBLE PRECISION COEF,S,EN,ENLN,HO,F
      DOUBLE PRECISION COEF,S,EN,ENLN,HO
C     SETS INITIAL ESTIMATES.
C
      COMMON/SPECES/COEF(2,7,100),S(100),EN(100,9),ENLN(100),HO(100)
     1,SUB(100,8),IUSE(100),TEMP(100,2),OMIT(8,100),ENSERT(8,100)
     2,JUSE(100),A(15,100)
      COMMON/MISC/ENN,SUMN,TT,LLMT(15),BO(15),TM,PP,R,ENNL,TRI,BASE(15)
     1,NAME(15,5),ANUM(15,5),TOP(15),ADD(15),VOL(100),GAMMA(100),NOCONV
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX
C
      ENN = .1
      ENNL = -2.3025851
      SUMN = ENN
      XI = NS-NC
      XI = ENN/XI
      XLN = ALOG(XI)
      DO 432 J=1,NS
      IF(IUSE(J).GT.0) IUSE(J)=-IUSE(J)
      EN(J,NPT) = 0.
      ENLN(J) = 0.
      IF(IUSE(J).NE.0) GO TO 432
      EN(J,NPT) = XI
      ENLN(J) = XLN
  432 CONTINUE
      IQ1 = NLM+1
      RETURN
      END
        SUBROUTINE FCOEFF
C
C       MODIFIED MAY 1987 TO INCLUDE SAXENA AND FEI CORRESPONDING STATES
C       METHOD ABOVE 4000 BARS.
C       MODIFIED AUGUST 1987 TO INCLUDE CHCL3, CH2CL2, CHCL3,CCL4, AND COCL2
C     RETURNS LN(FUGACITY COEFF.) OF GASEOUS SPECIES.
C
C     THE FOLLOWING DOUBLE PRECISION TYPE STATEMENTS ARE REQUIRED FOR
C     IBM 360 MACHINE ONLY.
C
C#### F is not used
C      DOUBLE PRECISION COEF,S,EN,ENLN,HO,F
      DOUBLE PRECISION COEF,S,EN,ENLN,HO
C
      INTEGER SUB
C#### AAA and BBB are not used
C      DIMENSION NAMES(26),AAA(26),BBB(26)
      DIMENSION NAMES(26)
C
      COMMON/SPECES/COEF(2,7,100),S(100),EN(100,9),ENLN(100),HO(100)
     1,SUB(100,8),IUSE(100),TEMP(100,2),OMIT(8,100),ENSERT(8,100)
     2,JUSE(100),A(15,100)
      COMMON/MISC/ENN,SUMN,TT,LLMT(15),BO(15),TM,PP,R,ENNL,TRI,BASE(15)
     1,NAME(15,5),ANUM(15,5),TOP(15),ADD(15),VOL(100),GAMMA(100),NOCONV
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX
C
      DATA NAMES/4HAR  ,4HH2  ,4HN2  ,4HO2  ,4HF2  ,4HCL2 ,4HNH3 ,4HN2H4
     A,4HHCN ,4HCO  ,4HCO2 ,4HN2O ,4HNO  ,4HN2O4,4HSO2 ,4HH2S ,4HHF  ,
     B4HSIF4,4HHCL ,4HCH4 ,4HH2O ,4HCHCL,4HCH2C,4HCH3C,4HCCL4,4HCOCL/
C
C
C     NNAMES IS THE NUMBER OF BUILT IN SPECIES
C
      NNAMES = 26
      DO 51 K=1,NS
      IF(IUSE(K).NE.0) GO TO 51
      IF(PP.LT.1.) GO TO 60
      DO 10 J=1,NNAMES
      IJ = J
      IF(SUB(K,1).EQ.NAMES(J)) GO TO 20
   10 CONTINUE
   60 GAMMA(K)=0.0
      GO TO 51
   20 CALL FPCALC(TT,PP,XLNFP,IJ)
        GAMMA(K)=XLNFP
51    CONTINUE
        RETURN
        END
        SUBROUTINE FPCALC(TK,P,XLNFP,IJ)
C       USES MRK AT P BELOW 4 KB AND SAXENA AT P ABOVE 4 KB
C       CALCULATES NATURAL LOG OF FUGACITY COEFFICIENTS.
C       TK IS IN KELVINS, P IS IN BARS. CORRECTED NOV. 15, 1990
C       FOR CO2 A PARAMETER IN THE RK PROGRAM.
        PO=4000.
        IF (P.LT.4000.) GO TO 100
C       CALC FP AT T AND 4 KB.
        CALL RKCALC(TK,PO,FO,IJ)
        CALL SAXENA(TK,P,XLNF,IJ)
C       ADD REFERENCE LN FUGACITY TO LN F2/F1
        PUREG=FO+XLNF
50      CONTINUE
        GO TO 200
100     CALL RKCALC(TK,P,PUREG,IJ)
200     CONTINUE
C       CONVERT FROM LN FUGACITY TO LN FUGACITY COEFFICIENT
250     XLNFP=PUREG-ALOG(P)
        RETURN
        END
        SUBROUTINE RKCALC(T,P,XLNF,IJ)
C       LAST MODIFYED JUN 13, 1987
C       CALCULATION OF PURE GAS MRK PROPERTIES FOLLOWING HOLLOWAY
C       1981.
C       P IS IN BARS, T IN KELVINS. RETURNS NATURAL LOG FUGACITY.
C
        DIMENSION A(26),B(26)
C
C       THE FOLLOWING ARRAYS GIVE THE REDLICH-KWONG A AND B PARAMETERS
C       FROM HOLLOWAY 1987 IN BAR UNITS.
C       REDLICH-KWONG A AND B PARAMETERS, FROM HOLLOWAY 1987
C       THE A PARAMETER IS DIVIDED BY ONE MILLION
C               AR    H2    N2    O2    F2    CL2   NH3  N2H4  HCN
      DATA A/16.9, 1.81, 15.6, 17.5, 58.3, 13.6, 86.8, 100., 244.,
C         C0   C02   N2O    NO   N2O2    S02  H2S  HF  SIF4
     A 17.2, 83., 68.7, 19.8, 112.6, 133.1, 89.,78., 85.6,
C         HCL  CH4  H2O CH3CL CH2CL2 CHCL3 CCL4 COCL2
     B 67.4, 32., 88.,156.0, 248.8,359.1,472.1,229.6/
C               AR    H2    N2    O2    F2    CL2   NH3  N2H4  HCN
      DATA B/22.3,15.34, 26.8, 22.1, 33.7,  39., 25.9,  60., 61.,
C         C0   C02   N2O    NO   N2O2    S02  H2S  HF  SIF4
     A 27.4,29.7, 30.8,   20., 30.7,  34.7, 20.,12.8,50.1,
C         HCL    CH4   H2O CH3CL CH2CL2 CHCL3 CCL4 COCL2
     B 28.12, 29.7, 14.6, 44.8, 60.4, 70.6, 87.8, 57.8/
C
C       GAS CONSTANT IN CC BAR/DEG MOLE
C#### RR is not used
C        DATA R/82.05/,RR/6732.2/
        DATA R/82.05/
        PBLN=ALOG(P)
        PA=P*0.98717
        TCEL=T-273.15
        RXT=R*T
        RT=R*T**1.5*1.E-6
C       CALCULATE T DEPENDENT A PARAMETERS FOR H2O AND CO2
        AH2OM=166.8-.19308*TCEL+.1864E-3*TCEL*TCEL-.71288E-7*TCEL**3
        IF(TCEL.GT.1200.) AH2OM=140.-0.05*TCEL
        IF(TCEL.LT.400.) GO TO 130
        IF(TCEL.LT.600.) AH2OM=4.221E3-3.1227E1*TCEL+8.7485E-2*TCEL
     1**2-1.07295E-4*TCEL**3+4.86111E-8*TCEL**4
C       SUB 400 DEG PARAMETERS BY JAMES CONNOLLY 7/30/81
  130   IF(TCEL.LT.400.) AH2OM=(7.4174E1+2.132E-1*TCEL-4.421E-4*
     1TCEL**2+2.8623E-7*TCEL**3)*0.987
        ACO2M=73.03-0.0714*TCEL+2.157E-5*TCEL*TCEL
        BSUM=B(IJ)
        ASUM=A(IJ)
        IF(IJ.EQ.21) ASUM=AH2OM
        IF(IJ.EQ.11) ASUM=ACO2M
        ASUM=ASUM/(BSUM*RT)
        BSUM=PA*BSUM/RXT
C
C       CALL REDKW TO FIND ln (FUGACITY COEFFICIENT OF PURE GAS).
        CALL REDKW(BSUM,ASUM,FP)
C       CONVERT TO LN FUGACITY
        XLNF=FP+PBLN
101     CONTINUE
        RETURN
        END
        SUBROUTINE REDKW(BP,A2B,XLNFP)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        REAL BP,A2B,XLNFP
C       A ROUTINE TO CALC COMPRESSIBILITY FACTOR AND FUGACITY
C       COEFFICIENT WITH THE REDLICH-KWONG EQUATION FOLLOWING
C       EDMISTER(1968).  THIS SOLUTION FOR SUPERCRITICAL FLUID.
C       WRITTEN BY JOHN R. HOLLOWAY, 1973,  REVISED APRIL, 1983
C       RETURNS A VALUE OF 1 FOR FP IF ARGUEMENTS OUT OF RANGE.
C       MODIFIED TO PREVENT HIGH P BLOWUP 6/12/87
C       RETURNS THE NATURAL LOG OF  FUGACITY COEFFICIENT
C#### DELTA is not used
C        DATA TH/0.333333D0/, DELTA/0.05D0/
        DATA TH/0.333333D0/
        IF(A2B.LT.1.D-100) A2B=1.D-6
        RR=-A2B*BP*BP
        QQ=BP*(A2B-BP-1.)
        XN=QQ*TH+RR-0.074074
        XM=QQ-TH
        XNN=XN*XN/4.
        XMM=XM*XM*XM/27.
        ARG=XNN+XMM
        IF(ARG.GT.0.0) GO TO 5
        IF(ARG.LT.0.0) GO TO 15
        FP=1.
        Z=1.
        RETURN
  15    COSPHI=DSQRT(-XNN/XMM)
        IF(XN.GT.0.0) COSPHI=-COSPHI
        TANPHI=DSQRT(1.-COSPHI*COSPHI)/COSPHI
        PHI=DATAN(TANPHI)*TH
        FAC=2.*DSQRT(-XM*TH)
C       SORT FOR LARGEST ROOT
        R1=DCOS(PHI)
        R2=DCOS(PHI+2.0944)
        R3=DCOS(PHI+4.18879)
C       SORT FOR LARGEST ROOT
        RH=R2
        IF(R1.GT.R2) RH=R1
        IF(R3.GT.RH) RH=R3
        Z=RH*FAC + TH
        GO TO 20
   5    X=DSQRT(ARG)
        F=1.
        XN2=-XN/2.
        XMM=XN2+X
        IF(XMM.LT.0.0) F=-1.
        XMM=F*((F*XMM)**TH)
        F=1.
        XNN=XN2-X
        IF(XNN.LT.0.0) F=-1.
        XNN=F*((F*XNN)**TH)
        Z=XMM+XNN+TH
  20    ZBP=Z-BP
        IF(ZBP.LT.1.D-10) ZBP=1.D-10
        BPZ=1.+BP/Z
        DFP=Z-1.-DLOG(ZBP)-A2B*DLOG(BPZ)
        IF(DFP.LT.-100..OR.DFP.GT.60) DFP=1.D-6
        XLNFP=DFP
        RETURN
        END
        SUBROUTINE SAXENA(TK,PBAR,XLNF,IJ)
C       HIGH PRESSURE CORRESPONDING STATES ROUTINES FROM SAXENA AND FEI
C       (1987) GCA VOL. 51, 783-791
C       CORRECTED FOR H2O TERM.
C       RETURNS NATURAL LOG OF THE RATIO F(P)/F(4000 BAR) AS XLNF ARRAY
        REAL LNF
        DIMENSION TCC(26), PCC(26)
        DATA PO/4000./
C       CRITICAL TEMPERATURES AND PRESSURES.
C
C                 AR    H2    N2    O2    F2    CL2  NH3  N2H4
      DATA TCC/150.8,43.6,126.2,154.58,144.3,417.,405.6,653.,
C          HCN   C0   C02   N2O    NO   N2O2   S02  H2S
     A 456.8,132.9,304.2,309.56,180., 431., 430.8,373.2,
C         HF  SIF4   HCL   CH4   H2O   CH3CL CH2CL2 CHCL3  CCL4 COCL2
     B  461.,259.,324.6,190.7,647.14,416.3, 510., 536.4,556.4, 455./
C                 AR    H2    N2    O2    F2  CL2   NH3  N2H4
      DATA PCC/48.1,20.46,33.5,49.77,51.47, 76., 111.3, 145.,
C         HCN   C0    C02   N2O    NO  N2O2    S02   H2S  HF
     A 53.2, 34.95,73.9, 71.5,  64.,  100.,134.3, 78.8, 64.,
C         SIF4   HCL  CH4   H2O  CH3CL CH2CL2 CHCL3  CCL4 COCL2
     B 36.7,  82.,  46.4, 218.,65.92, 60.0, 54.0, 45.0, 56.0/
C
        PB=PBAR
      TR=TK/TCC(IJ)
      PR=PB/PCC(IJ)
        PC=PCC(IJ)
C
        IF (IJ.EQ.16.OR.IJ.EQ.21) GO TO 3050
C       COEFFICIENTS FOR FLUIDS IN GENERAL
      A=2.0614 - 2.2351/TR**2 - .39411*LOG(TR)
      B=.055125/TR + .039344/TR**2
      C=-1.8935E-06/TR - 1.1092E-05/TR**2 - 2.1892E-05/TR**3
      D=5.0527E-11/TR - 6.3033E-21/TR**3
        GO TO 3100
C     SAXENA AND FEI EQUATIONS FOR H2O (AND H2S)
3050  A=1.4937-1.8626/TR**2 +.80003*TR**(-3)-.3941*LOG(TR)
      B=.0400241/TR +.024097/TR**2 - .0089634*TR**(-3)
      C=-9.016E-07/TR - 6.1345E-05/TR**2
        C=C + 2.238E-05/TR**3+5.2335E-07*LOG(TR)
      D=-7.6707E-09/TR + 4.1108E-08*TR**(-2)
        D=D - 1.4798E-08/TR**3-6.3033E-21*TR**3
C       CALC MOLAR VOLUME
3100  CONTINUE
C       INTEGRATE FROM PO (4000 BARS) TO P
      LNF = A*LOG(PB/PO)+(B/PC)*(PB-PO)+(C/(2.*PC**2))*(PB**2-PO**2)
      LNF=LNF+(D/(3.*PC**3))*(PB**3-PO**3)
        XLNF=LNF
        RETURN
        END
      SUBROUTINE OUT1
C
C     PRINTS OUTPUT DATA.  LAST MODIFIED JULY 23, 1989 FOR IW
C
C#### F is not used
C      DOUBLE PRECISION COEF,S,EN,ENLN,HO,F,V,TD,DTRA,DLTRA
      DOUBLE PRECISION COEF,S,EN,ENLN,HO,V,TD,DTRA,DLTRA
C
      INTEGER SUB
      DIMENSION TC(9), V(9), FUG(9), QFMREL(9)
      LOGICAL GAS,NEW
C
      COMMON/POINTS/P(45),T(45),PPP(9),TTT(9),TOTN(9),RSIZE,TEM,GAS,NEW,
     1 GTOTN(9),FREE(9),OXFUG(9),FO2C(9),FO2H(9)
      COMMON/SPECES/COEF(2,7,100),S(100),EN(100,9),ENLN(100),HO(100)
     1,SUB(100,8),IUSE(100),TEMP(100,2),OMIT(8,100),ENSERT(8,100)
     2,JUSE(100),A(15,100)
      COMMON/MISC/ENN,SUMN,TT,LLMT(15),BO(15),TM,PP,R,ENNL,TRI,BASE(15)
     1,NAME(15,5),ANUM(15,5),TOP(15),ADD(15),VOL(100),GAMMA(100),NOCONV
      COMMON/INDX/NP,NT,NPT,NLM,NS,IMAT,IQ1,NOMIT,NSERT,NC,KASE
     1,NREAC,NMIN,NEWMOL,NRAT,ALTER,SYSTM,PLUS,NMIX
C
      DATA IB/' '/,IO2/'O2  '/  
C#### Test for FALSE use's ".NOT." instead of "FALSE"
C#### Also changed go to 15 to go to 1115, cannot go to end of do loop
C#### It causes program to go back into the do loop
C       IF(NEW.EQ.FALSE) GO TO 15
        IF(.NOT.NEW) GO TO 1115
      WRITE(6,3) KASE
    3 FORMAT (9H CASE NO. ,I8)
        WRITE(6,1110)
1110    FORMAT(' INPUT COMPOSITION ')
      NCOND = 0 
      DO 15 N=1,NREAC
      DO 13 J=1,5
      JT = J
      IF(NAME(N,J).EQ.IB) GO TO 14
   13 CONTINUE
      JT = 6
   14 J = JT-1
1111    FORMAT(' ',5(A2,1X,F7.4,1X),5X,F9.3)
      WRITE(6,1111) (NAME(N,JJ),ANUM(N,JJ),JJ=1,5),BASE(N)
   15 CONTINUE 
C#### Added 1115 statement to prevent going back into do loop
1115  CONTINUE
      WRITE (6,863)
  863 FORMAT(/' THERMODYNAMIC PROPERTIES')
C
C     PRESSURE.
C
      WRITE(6,111)(PPP(J),J=1,NPT)
111     FORMAT(/' P,BARS ',T14,6F9.0)
C
C     TEMPERATURE.
C
      DO 65 I=1,NPT
      TC(I)= TTT(I) - 273.15
   65 CONTINUE
      WRITE(6,112)(TC(J),J=1,NPT)
112     FORMAT(' T, DEG. C ',T14,6F9.0)
C
C     FREE ENERGY.
C
      WRITE(6,113) (FREE(J),J=1,NPT)
113     FORMAT(' TOT. DELTA G ',T14,6F9.3)
C
C       DEFINE LIMIT FOR TRACE SPECIES
      DTRA = 9.D-7
        DLTRA=DLOG10(DTRA)
C
C     MOLE FRACTIONS - EQUILIBRIUM.
C
      WRITE(6,310)
  310 FORMAT(/' MOLE FRACTIONS(TOTAL)')
      DO 330 K=1,NS
      DO 315 I=1,NPT
        TD=TOTN(I)
      V(I) = EN(K,I)/TD
  315 CONTINUE
      DO 316 I=1,NPT
      IF(V(I).GE.DTRA) GO TO 320
  316 CONTINUE
      GO TO 330
  320 WRITE(6,114) (SUB(K,J),J=1,3),(V(I),I=1,NPT)
114     FORMAT(' ',3A4,6F9.5)
      IF(IUSE(K).GT.0) NCOND = NCOND + 1
  330 CONTINUE
      IF(.NOT.GAS.OR.NCOND.EQ.0) GO TO 20
      WRITE(6,80)
   80 FORMAT(/' MOLE FRACTIONS(GASES)')
      DO 430 K=1,NS
      IF(IUSE(K).NE.0) GO TO 430
      DO 415 I=1,NPT
        TD=GTOTN(I)
  415 V(I) = EN(K,I)/TD
      DO 416 I=1,NPT
      IF(V(I).GE.DTRA) GO TO 420
  416 CONTINUE
      GO TO 430
  420 WRITE(6,114) (SUB(K,J),J=1,3),(V(I),I=1,NPT)
  430 CONTINUE
   20 WRITE(6,335)
  335 FORMAT(/' LOG(10) (MOLE FRACTION)')
C
      DO 350 K=1,NS
      DO 340 I=1,NPT
      IF(GAS)GO TO 341
        TD=TOTN(I)
      V(I)=EN(K,I)/TD
      GO TO 340
  341 TD=GTOTN(I)
        V(I)=EN(K,I)/TD
  340 CONTINUE
      DO 370 I=1,NPT
      IF(V(I).EQ.0.) GO TO 370
      V(I)=DLOG10(V(I))
      IF(SUB(K,1).EQ.IO2) FUG(I) = V(I)+ALOG10(PPP(I))
  370 CONTINUE
      DO 365 J=1,NPT
      IF(V(J).LT.DLTRA) GO TO 360
  365 CONTINUE
      GO TO 350
  360 WRITE(6,115)(SUB(K,J),J=1,3),(V(I),I=1,NPT)
115     FORMAT(' ',3A4,6F9.3)
350     CONTINUE
C       OXYGEN FUGACITY CALCULATIONS
      DO 380 N = 1,NPT
C      QFM = 9.0-25738./TTT(N)+0.092*(PPP(N)-1.)/TTT(N)
C       THIS VERSION USES QFI (HUEBNER) AS REF VALUE
      XIW = 7.51-29382./TTT(N)+0.050*(PPP(N)-1.)/TTT(N)
        QFMREL(N) = XIW
      IF(PPP(N).LT.1.0.OR.V(N).LT.1.D-100) GO TO 380
        IF(OXFUG(N).GT.1.E4) OXFUG(N)=1.E4
      FUG(N) = FUG(N)+OXFUG(N)
  380 CONTINUE
        WRITE(6,121)
121     FORMAT(/' OXYGEN FUGACITY CALCULATIONS:')
      WRITE(6,116) (FUG(N),N=1,NPT)
116     FORMAT(' FO2=XO2*FGO2',T14,6F9.3)
      WRITE(6,119) (FO2C(N),N=1,NPT)
119     FORMAT(' FO2 = CO/CO2',T14,6F9.3)
      WRITE(6,120) (FO2H(N),N=1,NPT)
120     FORMAT(' FO2 = H2/H20',T14,6F9.3)
        WRITE(6,117) (QFMREL(N),N=1,NPT)
117     FORMAT(' QFI FO2',T14,6F9.3)
      RETURN
      END