PROGRAM KCAL C =====================================================================| C A1(i) = | C A2(i) = | C A3(i) = | C A4(i) = Coefficients in the | C C1(i) = Equation of State for aqueous ions | C C2(i) = | C GREF(i) = Standard state free energy of formation for the | C species at 298.15 K and 1 bar in Cal/mole. | C HREF(i) = Standard state enthalpy of formation for the species | C at 298.15 K and 1 bar in Cal/mole. | C ION(i) = The name of the aqueous ion | C INDP(k) = k = 1, ITOTP index array of selected pressures | C INDT(j) = j = 1, ITOTT index array of selected temperatures | C ITOTP = Total number of selected pressures | C ITOTT = Total number of selected temperatures | C NI = The maxium # of species. | C PMAX = The maximum number of pressure grids for water | C SREF(i) = The third law entropy for the species at 298.15 K | C and 1 bar in Cal/(mole K). | C T(j) = The temperature array in kelvin | C TC(j) = The temperature array in celcius | C TMAX = The maximum number of temperature grids for water | C ---------------------------------------------------------------------| INTEGER*2 NG, NI, PM, TM PARAMETER (NG = 10, NI = 40, PM = 5, TM = 24) INTEGER*2 I, I2, I3, I4, IFLAG, IMAX, IND, INDP(PM), INDT(TM) INTEGER*2 INDWAT, INP1, INP2, INP3, INP4, INP5, IOUT1, IOUT2 INTEGER*2 IOUT3, ITOTP, ITOTT, ITTY INTEGER*2 J, K, LN, M, M1, MAXG, NSGAS, NSION, NSMIN, PMAX INTEGER*2 SPEC(NI), SPECG(NG), SPECM(8), TMAX, Z(NI) REAL*8 A1(NI), A2(NI), A3(NI), A4(NI) REAL*8 ACTI(10), ACTM(8), ACTW REAL*8 AG(NG), BG(NG) REAL*8 C1(NI), C2(NI), CG(NG), CPREF(NI) REAL*8 DGDP(TM,PM), DGDT(TM,PM), DG2DT2(TM,PM), DIEL(TM,PM) REAL*8 G(NI,TM,PM), GG(10,21), GGREF(NG), GM(8,21,10) REAL*8 GREF(NI), GRF, GW(TM,PM), GWREF, HCT(0:2,7) REAL*8 HGREF(NG), HREF(NI), HRF, HTRAN1, HTRAN2, LITG(TM,PM) REAL*8 P(TM,PM), PREF, Q(TM,PM) REAL*8 SGREF(NG), SLOPE1, SLOPE2, SREF(NI), SRF REAL*8 T(TM), TC(TM), TREF, TTRAN1, TTRAN2, TTRAN3 REAL*8 VIREF(NI), VRE1, VRE2, VREF, WREF(NI) REAL*8 X(TM,PM), XLOGK(24), XMOLG(8), XRG(8), XRAD(NI) REAL*8 XRI(10), XRM(8), XRW, Y(TM,PM) CHARACTER* 1 ANS CHARACTER* 8 BLANK, GAS(NG,2), ION(NI), MIN(2), MINAM(8) CHARACTER* 8 RENAM1(10), RENAM2(10) CHARACTER*12 ENAME(8) CHARACTER*80 DUMMY, DUMMY2 COMMON /DATGAS/ SGREF, GGREF, AG, BG, CG, HGREF COMMON /DATION/ A1, A2, A3, A4, C1, C2, CPREF, GREF, HREF, & SREF, VIREF, WREF, XRAD, Z COMMON /DATMIN/ GRF, HRF, SRF, VREF, HCT, & TTRAN1, HTRAN1, VRE1, SLOPE1, & TTRAN2, HTRAN2, VRE2, SLOPE2, & TTRAN3 COMMON /GIBBS / G, GG, GM, GW, GWREF COMMON /IO / INP1, INP4, INP5, IOUT1, IOUT2, & IOUT3, ITTY COMMON /MAX / TMAX, PMAX, INDP, INDT, ITOTT, ITOTP, MAXG COMMON /NOMS / GAS, ION, MINAM COMMON /STOCH / ACTI, ACTM, XMOLG, XRI, XRG, & XRM, XRW, ACTW, NSION, NSGAS, NSMIN, INDWAT, & SPEC, SPECG COMMON /STUFF / P, PREF, T, TC, TREF COMMON /WATER / DIEL, LITG, Q, X, Y, DGDT, DGDP, DG2DT2 EXTERNAL AQUION, GASES, REAC, READER, README, WRITGA, WRITIO DATA BLANK /' '/ DATA ENAME /'NION.DAT ', 'NWATER.DAT ', 'NGAS.DAT ', & 'NMINERAL.DAT', 'SOLTHERM.DAT', 'KCALOUT.DAT ', & 'KLOG.DAT ', 'KDIFF.DAT '/ C =====================================================================| C Fortran I/O variables: | C INP1: Input file w/ ion data | C INP2: Input file w/ p-t data for water | C INP3: Input file w/ gases thermo data | C INP4: Input file w/ minerals thermo data | C INP5: SOLTHERM data file | C IOUT1: Main output | C IOUT2: SOLTHERM formatted output with new log k's | C IOUT3: Scratch file with differences between new and | C old log k's | C ITTY: TTY input or output | C ---------------------------------------------------------------------| INP1 = 10 INP2 = 11 INP3 = 12 INP4 = 13 INP5 = 14 ITTY = 5 IOUT1 = 6 IOUT2 = 15 IOUT3 = 16 OPEN (INP1, FILE = ENAME(1), STATUS = 'OLD') OPEN (INP2, FILE = ENAME(2), STATUS = 'OLD') OPEN (INP3, FILE = ENAME(3), STATUS = 'OLD') OPEN (INP4, FILE = ENAME(4), STATUS = 'OLD') OPEN (INP5, FILE = ENAME(5), STATUS = 'OLD') OPEN (IOUT1, FILE = ENAME(6), STATUS = 'UNKNOWN') OPEN (IOUT2, FILE = ENAME(7), STATUS = 'UNKNOWN') OPEN (IOUT3, FILE = ENAME(8), STATUS = 'UNKNOWN') TREF = 298.15D+00 PREF = 1.0D+00 GWREF = -56686.0D+00 IMAX = NI MAXG = NG DO 20 I = 1, TM T(I) = 0.0 TC(I) = 0.0 DO 10 J = 1, PM P(I,J) = 0.0 10 CONTINUE 20 CONTINUE C ---------------------------------------------------------------------| C INDWAT is the index of water | C ---------------------------------------------------------------------| INDWAT = 40 C ---------------------------------------------------------------------| C Read data for the ions | C ---------------------------------------------------------------------| READ (INP1, '(I5)') IMAX IF (IMAX .GT. NI) THEN PRINT '(A, A12, A, I2)', ' The number of ions selected in ', & ENAME(1), ' must be <= ', NI CLOSE (INP1) GOTO 700 END IF READ (INP1, '(A80)') DUMMY READ (INP1, '(A80)') DUMMY DO 30 I = 1, IMAX READ (INP1, 1000, END = 40) ION(I), Z(I), GREF(I), HREF(I), & SREF(I), CPREF(I), VIREF(I), A1(I) IF (ION(I) .EQ. BLANK) GO TO 40 READ (INP1, 1010) XRAD(I), A2(I), A3(I), A4(I), & C1(I), C2(I), WREF(I) 30 CONTINUE 40 CONTINUE PRINT *, ' The data for the ions has been read...' PRINT * C ---------------------------------------------------------------------| C Reads data for water at given P and T | C ---------------------------------------------------------------------| READ (INP2, '(A80)') DUMMY READ (INP2, '(A80)') DUMMY READ (INP2, '(2I5)') TMAX, PMAX IF (TMAX .GT. TM .OR. PMAX .GT. PM) THEN PRINT '(A, A, I2, /, A, A, I2, A, A12)', & ' The number of temperatures must ', & 'be less than ', TM, ' and the number of pressures must ', & 'be less than ', PM, ' in ', ENAME(2) CLOSE (INP2) GOTO 700 END IF ITOTT = TMAX ITOTP = PMAX DO 60 K = 1, PMAX DO 50 J = 1, TMAX READ (INP2, 1070) P(J,K), TC(J), GW(J,K), DIEL(J,K), Q(J,K), & X(J,K), Y(J,K), LITG(J,K), DGDT(J,K), & DGDT(J,K), DG2DT2(J,K) 50 CONTINUE 60 CONTINUE READ (INP2, '(A80)') DUMMY READ (INP2, '(A80)') DUMMY2 DO 70 J = 1, TMAX T(J) = TC(J) + 273.15 70 CONTINUE PRINT *, ' The data for water has been read....' PRINT * PRINT *, ' Do you want to echo it on the screen ? ' READ (ITTY, '(A1)') ANS IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') THEN DO 90 K = 1, PMAX WRITE (ITTY, 1020) DO 80 J = 1, TMAX WRITE (ITTY, 1030) P(J,K), TC(J), GW(J,K), DIEL(J,K), & Q(J,K), X(J,K), Y(J,K) 80 CONTINUE PRINT *, ' ' READ '(A1)', ANS 90 CONTINUE END IF C ---------------------------------------------------------------------| C Reads data for gases | C ---------------------------------------------------------------------| READ (INP3, '(I2)') MAXG READ (INP3, '(/)') DO 100 M = 1, MAXG READ (INP3, 1040, END = 110) (GAS(M,M1), M1 = 1, 2), GGREF(M), & HGREF(M), SGREF(M), AG(M), BG(M), & CG(M) 100 CONTINUE 110 CONTINUE C IF (M .LT. MAXG) MAXG = M - 1 PRINT *, ' The data for the gases has been read...' PRINT * C ---------------------------------------------------------------------| C Start of calculation loop | C ---------------------------------------------------------------------| 200 CONTINUE C ---------------------------------------------------------------------| C The following displays the T-P grid and give the user the | C choice of P-T to be used in the calculation | C ---------------------------------------------------------------------| PRINT *, 'The following pressures and temperatures have ', & 'been chosen for the default:' PRINT * READ (DUMMY, '(10I5)') (INDP(K), K = 1, PMAX) DO 203 K = 1, PMAX IF (INDP(K) .GT. 0) ITOTP = K 203 CONTINUE READ (DUMMY2, '(24I5)') (INDT(J), J = 1, TMAX) DO 208 J = 1, TMAX IF (INDT(J) .GT. 0) ITOTT = J 208 CONTINUE WRITE (ITTY, 3080) (P(1,INDP(K)), K = 1, ITOTP) WRITE (ITTY, 3090) (TC(INDT(I)), I = 1, ITOTT) PRINT * PRINT *, ' Do you want to carry calculations on the default', & ' P-T grid ? ' READ (ITTY, '(A1)') ANS IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GO TO 500 PRINT *, ' INDEX P bars' DO 210 K = 1, PMAX INDP(K) = 0 WRITE (ITTY, 1050) K, P(1,K) 210 CONTINUE PRINT *, ' (P = 1 bar means H2O saturation pressure)' PRINT *, ' Select desired pressure indexes (end by slash)' PRINT * READ (ITTY, *) (INDP(K), K = 1, PMAX) WRITE (ITTY, 1060) DO 220 I = 1, 8 I2 = I + 8 I3 = I + 16 INDT(I) = 0 INDT(I2) = 0 WRITE (ITTY, 1050) I, TC(I), I2, TC(I2), I3, TC(I3) 220 CONTINUE PRINT *, ' Select desired temperature indexes (end by slash)' PRINT * READ (ITTY, *) (INDT(I), I = 1, TMAX) DO 230 K = 1, PMAX IF (INDP(K) .EQ. 0) GO TO 240 ITOTP = K 230 CONTINUE 240 CONTINUE IF (ITOTP .NE. PMAX) ITOTP = K - 1 DO 250 I = 1, TMAX IF (INDT(I) .EQ. 0) GO TO 260 ITOTT = I 250 CONTINUE 260 CONTINUE IF (ITOTT .NE. TMAX) ITOTT = I - 1 PRINT *, 'The following pressures and temperatures have ', & 'been chosen:' WRITE (ITTY, 3080) (P(1,INDP(K)), K = 1, ITOTP) WRITE (ITTY, 3090) (TC(INDT(I)), I = 1, ITOTT) PRINT *, 'Is this correct ? ' READ (ITTY, '(A1)') ANS IF (ANS .EQ. 'N' .OR. ANS .EQ. 'n') GO TO 200 500 CONTINUE PRINT *, 'Do you want to read stochiometries from SOLTHERM ? ' READ (ITTY, '(A1)') ANS IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') THEN CALL READER GOTO 700 END IF C ---------------------------------------------------------------------| C The following reads stochiometry for a given reaction involving | C minerals and/or gases and/or ions interactively | C ---------------------------------------------------------------------| WRITE (ITTY, 2010) READ (ITTY, '(10A8)') (RENAM1(I), I = 1, 10) READ (ITTY, '(10A8)') (RENAM2(I), I = 1, 10) WRITE (IOUT1, '(1A1)') '1' WRITE (IOUT1, '(2X, 10A8)') (RENAM1(I), I = 1, 10) WRITE (IOUT1, '(2X, 10A8)') (RENAM2(I), I = 1, 10) WRITE (IOUT1, 1080) 510 CONTINUE WRITE (ITTY, 2030) READ (ITTY, *) NSION, NSGAS, NSMIN WRITE (ITTY, 2020) NSION, NSGAS, NSMIN READ (ITTY, '(A1)') ANS IF (ANS .NE. 'Y' .AND. ANS .NE. 'y') GO TO 510 IF (NSION .EQ. 0) GO TO 550 C ---------------------------------------------------------------------| C Reads stoichiometry for aqueous species, then calculate free | C energy | C ---------------------------------------------------------------------| 520 CONTINUE DO 530 I = 1, 10 I2 = I + 10 I3 = I + 20 I4 = I + 30 WRITE (ITTY, 2040) I, ION(I), I2, ION(I2), & I3, ION(I3), I4, ION(I4) 530 CONTINUE WRITE (ITTY, 2060) INDWAT WRITE (ITTY, 2050) READ (ITTY, *) (XRI(I), SPEC(I), ACTI(I), I = 1, NSION) DO 540 I = 1, NSION IF (ACTI(I) .EQ. 0.0) ACTI(I) = 1.0 IF (SPEC(I) .NE. INDWAT) THEN WRITE (ITTY, 1090) XRI(I), ION(SPEC(I)), ACTI(I) ELSE XRW = XRI(I) ACTW = ACTI(I) WRITE (ITTY, 2000) XRW, ACTW END IF 540 CONTINUE WRITE (ITTY, 2070) READ (ITTY, '(A1)') ANS IF (ANS .NE. 'Y' .AND. ANS .NE. 'y') GO TO 520 CALL WRITIO CALL AQUION C ---------------------------------------------------------------------| C Reads stoichiometry for gases, then calculate free energy | C ---------------------------------------------------------------------| 550 CONTINUE IF (NSGAS .EQ. 0) GO TO 580 DO 560 I = 1, MAXG WRITE (ITTY, 2040) I, GAS(I,1) 560 CONTINUE WRITE (ITTY, 2080) READ (ITTY, *) (XRG(I), SPECG(I), XMOLG(I), I = 1, NSGAS) DO 570 I = 1, NSGAS IF (XMOLG(I) .EQ. 0.0) XMOLG(I) = 1.0 WRITE (ITTY, 1090) XRG(I), GAS(SPECG(I),1), XMOLG(I) 570 CONTINUE WRITE (ITTY, 2070) READ (ITTY, '(A1)') ANS IF (ANS .NE. 'Y' .AND. ANS .NE. 'y') GO TO 550 CALL WRITGA CALL GASES C ---------------------------------------------------------------------| C Reads stoichiometry for minerals | C ---------------------------------------------------------------------| 580 CONTINUE IF (NSMIN .EQ. 0) GO TO 670 DO 590 I = 1, NSMIN WRITE (ITTY, 2090) READ (ITTY, '(A8)') MINAM(I) WRITE (ITTY, 3010) MINAM(I) READ (ITTY, *) XRM(I), ACTM(I) 590 CONTINUE WRITE (ITTY, 3000) READ (ITTY, '(A1)') ANS IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GO TO 580 C ---------------------------------------------------------------------| C Reads only the data for minerals that are specified then | C calculates free energy | C ---------------------------------------------------------------------| IFLAG = 0 595 CONTINUE READ (INP4, '(A8)') MIN(1) IF (MIN(1) .NE. '********') GO TO 595 DO 600 I = 1, 8 SPECM(I) = 0 600 CONTINUE 610 CONTINUE IF (IFLAG .GE. NSMIN) GO TO 660 READ (INP4, '(2A8)', END = 640) (MIN(LN), LN = 1, 2) DO 620 I = 1, NSMIN IF (MINAM(I) .EQ. MIN(1)) GO TO 630 620 CONTINUE READ (INP4, '(//////)', END = 640) GO TO 610 630 CONTINUE IND = I IFLAG = IFLAG + 1 SPECM(I) = IFLAG WRITE (ITTY, 3020) MINAM(I), XRM(I), ACTM(I) CALL README (MIN) CALL SOLIDS (IND) GO TO 610 640 CONTINUE DO 650 I = 1, NSMIN IF (SPECM(I) .EQ. 0) WRITE (ITTY, 3040) MINAM(I) 650 CONTINUE REWIND INP4 WRITE (IOUT1, 3050) GO TO 580 660 CONTINUE WRITE (ITTY, 3030) READ (ITTY, '(A1)') ANS IF (ANS .NE. 'Y' .AND. ANS .NE. 'y') THEN REWIND INP4 WRITE (IOUT1, 3050) GO TO 580 END IF 670 CONTINUE CALL REAC (XLOGK) C ---------------------------------------------------------------------| C Univariant curve not programmed in yet | C ---------------------------------------------------------------------| C WRITE (ITTY, 754) C 754 FORMAT (///5X, 'Do you want to calculate univariant curve ?'/) C READ (ITTY, '(A1)') ANS C IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') CALL UNIC WRITE (ITTY, 3060) READ (ITTY, '(A1)') ANS IF (ANS .NE. 'Y' .AND. ANS .NE. 'y') GOTO 700 REWIND INP4 WRITE (ITTY, 3070) READ (ITTY, '(A1)') ANS IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GO TO 500 DO 680 K = 1, PMAX INDP(K) = K 680 CONTINUE DO 690 J = 1, TMAX INDT(J) = J 690 CONTINUE ITOTP = PMAX ITOTT = TMAX GO TO 200 700 CONTINUE STOP 1000 FORMAT (A8, I4, 6E11.4) 1010 FORMAT (7X, F5.3, 6E11.4) 1020 FORMAT (6X, 'P', 8X, 'T', 8X, 'G', 9X, 'DIEL', & 7X, 'Q', 11X, 'X', 11X, 'Y',/) 1030 FORMAT (1X, F11.5, 1X, F5.0, 1X, 5(E11.5, 1X)) 1040 FORMAT (2A8, 4X, 6F10.2) 1050 FORMAT (5X, 3(I2, ':', 1X, F10.0, 5X)) 1060 FORMAT (4X, 3('INDEX T DEG.C', 5X), /) 1070 FORMAT (F8.3, F5.0, 1X, E11.4, 8(E11.3)) 1080 FORMAT (1X, 120('+')) 1090 FORMAT (5X, F6.3, 1X, ':', A8, 2X, F4.2) 2000 FORMAT (5X, F6.3, ' :H2O ', F4.2) 2010 FORMAT (//5X, 'Input reaction name and/or title (2 lines)'/) 2020 FORMAT (/5X, 'There are', I2, ' aqueous species, ', I2, & ' gases and ', I2, ' minerals involved in this ' & 'reaction. ', /5X, 'Is this correct ?'/) 2030 FORMAT (/5X, 'Input total # of aqueous species, gases and ' & 'minerals in the reaction.' & /5X, ' (Count water as an aqueous species) '/) 2040 FORMAT (5X, 4(I3, 1X, A8, 6X)) 2050 FORMAT (//5X, 'Input stoichiometric coefficients, indexes ', & 'and activities for ', & /5X, 'the ions (1 line). [Note: ', & 'Reactants should have negative coeffs.]',/) 2060 FORMAT (5X,' The current index for water is: ', I3) 2070 FORMAT (/5X, 'Is this correct ?'/) 2080 FORMAT (//5X, 'Input stoichiometric coefficients, indexes ' & 'and mole fraction', & /5X, 'for the gases (1 line)'/) 2090 FORMAT (//5X, 'Input mineral name as it appears in data file'/) 3000 FORMAT (//5X, 'Do you want to change data for the minerals ?') 3010 FORMAT (8X, 'Input coefficient and activity of ',A8) 3020 FORMAT (5X, 'Found : ', A8, ' Coef= ', F6.3, ' Act= ', F4.2) 3030 FORMAT (/5X, 'That''s it ! Is this correct ?'/) 3040 FORMAT (/5X, 'I cannot find the following mineral : ', A8) 3050 FORMAT (/10X, '***** Correction, New minerals given below'/) 3060 FORMAT (//5X, 'Do you want to calculate another reaction ?'/) 3070 FORMAT (5X, 'Do you want to keep the same P-T grid ?'/) 3080 FORMAT (5X, 'PRESSURES: ', 6(F7.2, 2X) / 19X, 4(F7.2, 2X)/) 3090 FORMAT (5X, 'TEMPERATURES: ', 6(F7.2, 2X) / 3(19X, 6(F7.2, 2X)/)/) END SUBROUTINE AQUION C =====================================================================| C *AQUION* calculates thermodynamic data at P and T for aqueous | C species using the revised HKF equation of state (Tanger & | C Helgeson, 1988; Shock and Helgeson, 1988; Shock, Helgeson & | C Sverjensky, 1989). Apparant standard state free energies are | C transferred back to *KCAL*. The standard states are at T and P | C for the ideal one molal solution. | C | C Originally written by R. STOESSEL | C Revised by J.D. DeBraal (4,92) | C | C GREF = Standard state free energy of formation for the | C species at 298.15 K and one bar. | C SREF = Third law entropy for the species at 298.15 K and | C one bar. | C TREF = Reference temperature (298.15 K) | C PREF = Reference pressure (one bar) | C | C T = Temperature = T(J) in kelvin | C P = Pressure = P(J,K) in bars | C | C C1 = " | C C2 = " | C THETA = Equation of state coefficients | C A1 = " | C A2 = " | C A3 = " | C A4 = " | C | C DIEL = The dielectric constant of water at T and P | C DIREF = The dielectric constant of water at TREF and PREF | C | C These properties are for the standard state: | C G = Apparent free energy for the aqueous specie at T and P | C S = Third law entropy of the aqueous species at T and P | C H = Apparent enthalpy of the aqueous species at T and P | C V = Partial molal volume at T and P in cal/mol/bar | C VCC = Partial molal volume at T and P in cc/mol | C CP = Partial molal heat capacity at T and P | C | C Y = Born Function defined by eqn. D-1, T & H | C X = Born Function defined by eqn. D-4, T & H | C Q = Born Function defined by eqn. D-2, T & H | C YREF = Y evaluated at Pref and Tref | C | C WREF = Born Coefficient at Tref and Pref | C WCALC = Calculated Born Coefficient at T and P defined by | C eqn. B-9, T & H | C X1 = d w / d g defined by eqn. B-10, T & H | C X2 = d2w / d2g defined by eqn. B-11, T & H | C | C | C The standard state apparent free energy of H+ is 0.0 at any | C T and P. | C ---------------------------------------------------------------------| INTEGER*2 NG, NI, PM, TM PARAMETER (NG = 10, NI = 40, PM = 5, TM = 24) INTEGER*2 I, INDP(PM), INDT(TM), INDWAT INTEGER*2 INP1, INP4, INP5, IOUT1, IOUT2, IOUT3 INTEGER*2 ITOTP, ITOTT, ITTY INTEGER*2 J, K, MAXG, NJ, NK, NS, NSGAS, NSION INTEGER*2 NSMIN, PMAX, SPEC(NI), SPECG(NG), TMAX, Z(NI) CHARACTER*8 GAS(NG,2), ION(NI), MINAM(8) REAL*8 A1(NI), A2(NI), A3(NI), A4(NI) REAL*8 ACTI(10), ACTM(8), ACTW REAL*8 C1(NI), C2(NI) REAL*8 CP, CPA, CPB, CPREF(NI) REAL*8 DIEL(TM,PM), DIEREF, DGDP(TM,PM), DGDT(TM,PM) REAL*8 DG2DT2(TM,PM) REAL*8 G(NI,TM,PM), G1, G2, G3, G4 REAL*8 GG(10,21), GM(8,21,10) REAL*8 GREF(NI), GW(TM,PM), GWREF REAL*8 H, HREF(NI), LITG(TM,PM) REAL*8 OMEGA, P(TM,PM), PREF, PSI, Q(TM,PM) REAL*8 S, S1, S2, S3 REAL*8 SREF(NI) REAL*8 T(24), TC(24), THETA, TREF REAL*8 V, VA, VB, VCC, VIREF(NI), WREF(NI), Y(TM,PM) REAL*8 YREF, X(TM,PM), X1, X2, XMOLG(8), XRAD(NI), XRG(8) REAL*8 XRI(10), XRM(8), XRW COMMON /DATION/ A1, A2, A3, A4, C1, C2, CPREF, GREF, HREF, & SREF, VIREF, WREF, XRAD, Z COMMON /GIBBS / G, GG, GM, GW, GWREF COMMON /IO / INP1, INP4, INP5, IOUT1, IOUT2, & IOUT3, ITTY COMMON /MAX / TMAX, PMAX, INDP, INDT, ITOTT, ITOTP, MAXG COMMON /NOMS / GAS, ION, MINAM COMMON /STOCH / ACTI, ACTM, XMOLG, XRI, XRG, & XRM, XRW, ACTW, NSION, NSGAS, NSMIN, INDWAT, & SPEC, SPECG COMMON /STUFF / P, PREF, T, TC, TREF COMMON /WATER / DIEL, LITG, Q, X, Y, DGDT, DGDP, DG2DT2 EXTERNAL WCALC C =====================================================================| DIEREF = 78.46D+00 PSI = 2600.00D+00 THETA = 228.00D+00 YREF = -5.81D-05 DO 100 NK = 1, ITOTP K = INDP(NK) WRITE (IOUT1, '(/)') WRITE (IOUT1, 510) WRITE (IOUT1, 515) DO 90 NJ = 1, ITOTT J = INDT(NJ) WRITE (IOUT1, 525) TC(J), P(J,K) DO 80 NS = 1, NSION I = SPEC(NS) IF (I .EQ. INDWAT) THEN WRITE (IOUT1, 500) GW(J,K) GO TO 80 END IF IF (Z(I) .EQ. 0) THEN OMEGA = WREF(I) X1 = 0.0D+00 X2 = 0.0D+00 ELSE CALL WCALC (OMEGA, X1, X2, LITG(J,K), XRAD(I), Z(I)) END IF C ---------------------------------------------------------------------| C Free Enegry (Eqn. 44 from Shock & Helegon) | C ---------------------------------------------------------------------| G1 = GREF(I) - SREF(I) * (T(J) - TREF) G1 = G1 - C1(I) * (T(J) * ALOG(T(J) / TREF) - T(J) + TREF) G1 = G1 + A1(I) * (P(J,K) - PREF) G1 = G1 + A2(I) * ALOG((PSI + P(J,K))/(PSI + PREF)) G2 = (1 / (T(J) - THETA)) - (1 / (TREF - THETA)) G2 = G2 * ((THETA - T(J)) / THETA) G3 = (TREF * (T(J) - THETA)) / (T(J) * (TREF - THETA)) G3 = T(J) / THETA**2 * ALOG(G3) G2 = G2 - G3 G1 = G1 - C2(I) * G2 G4 = (PSI + P(J,K)) / (PSI + PREF) G4 = A3(I) * (P(J,K) - PREF) + A4(I) * ALOG(G4) G1 = G1 + (1.0 / (T(J) - THETA)) * G4 IF (Z(I) .EQ. 0) THEN G1 = G1 + WREF(I) * YREF * (T(J) - TREF) G2 = WREF(I) * (1.0 / DIEL(J,K) - 1.0) G3 = WREF(I) * (1.0 / DIEREF - 1.0) G1 = G1 + G2 - G3 ELSE G1 = G1 + WREF(I) * YREF * (T(J) - TREF) G2 = OMEGA * (1.0 / DIEL(J,K) - 1.0) G3 = WREF(I) * (1.0 / DIEREF - 1.0) G1 = G1 + G2 - G3 END IF G(I,J,K) = G1 C ---------------------------------------------------------------------| C Entropy (Eqn. 38 from Shock & Helegson) | C ---------------------------------------------------------------------| S1 = SREF(I) + C1(I) * ALOG(T(J) / TREF) S3 = (TREF * (T(J) - THETA)) / (T(J) * (TREF - THETA)) S2 = 1.0 / (T(J) - THETA) - 1.0 / (TREF - THETA) S2 = S2 + 1.0 / THETA * ALOG(S3) S1 = S1 - C2(I) / THETA * S2 S3 = (PSI + P(J,K)) / (PSI + PREF) S3 = A3(I) * (P(J,K) - PREF) + A4(I) * ALOG(S3) S2 = (1.0 / (T(J) - THETA))**2 S1 = S1 + S2 * S3 IF (Z(I) .EQ. 0) THEN S1 = S1 + WREF(I) * (Y(J,K) - YREF) ELSE S1 = S1 + OMEGA * Y(J,K) S1 = S1 - (1.0 / DIEL(J,K) - 1) * (DGDT(J,K) * X1) S1 = S1 - WREF(I) * YREF END IF S = S1 C ---------------------------------------------------------------------| C Enthalpy | C ---------------------------------------------------------------------| H = HREF(I) + G(I,J,K) - GREF(I) + T(J)*S - TREF * SREF(I) C ---------------------------------------------------------------------| C Heat Capacity (Eqn. 36 from Shock & Helegson) | C ---------------------------------------------------------------------| CPA = C1(I) + C2(I) / (T(J) - THETA)**2 CPB = A4(I) * DLOG((PSI + P(J,K)) / (PSI + PREF)) CPB = A3(I) * (P(J,K) - PREF) + CPB CPB = CPB * (2.0 * T(J)) / (T(J) - THETA)**3 CPA = CPA - CPB IF (Z(I) .EQ. 0) THEN CP = CPA + WREF(I) * T(J) * X(J,K) ELSE CPA = CPA + OMEGA * T(J) * X(J,K) CPA = CPA + 2.0 * T(J) * Y(J,K) * (DGDT(J,K) * X1) CPB = DGDT(J,K)**2 * X2 + DG2DT2(J,K) * X1 CP = CPA - T(J) * (1 / DIEL(J,K) - 1) * CPB END IF C ---------------------------------------------------------------------| C Volume (Eqn. 11 from Shock & Helegson) | C ---------------------------------------------------------------------| VA = A1(I) + A2(I) * (1 / (PSI + P(J,K))) VB = A3(I) + A4(I) * (1 / (PSI + P(J,K))) VB = VB * (1 / (T(J) - THETA)) VA = VA + VB IF (Z(I) .EQ. 0) THEN V = VA + WREF(I) * Q(J,K) ELSE VA = VA - OMEGA * Q(J,K) VB = (1 / DIEL(J,K) - 1) * (DGDP(J,K) * X1) V = VA + VB END IF VCC = V * 41.8393 WRITE (IOUT1, 520) ION(I), G(I,J,K), H, S, CP, VCC 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN 500 FORMAT (1X, 'H2O ', 2X, F12.0, 13X, '--', 14X, '--', & 13X, '--', 16X, '--') 505 FORMAT (///) 510 FORMAT (1X, /, 1X, 'Species', 2X, 'App. Free Energy', & 2X, 'App. Enthalpy', 2X, 'Third Law Entropy', & 2X, 'Heat Capacity', 2X, 'Partial Molal Volume') 515 FORMAT (14X, 'Cal/Mole', 9X, 'Cal/Mole', & 7X, 'Cal/Deg/Mole', 5X, 'Cal/Deg/Mole', 9X, 'CC/Mole', & /1X, 96('-')) 520 FORMAT (1X, A8, 2X, F12.0, 4X, F12.0, 4X, F12.2, & 4X, F12.3, 6X, F12.3) 525 FORMAT (/,2X, F6.2,' Deg.C', 4X, F9.4,' Bars') END SUBROUTINE WCALC (WABS, X1, X2, LITG, XRAD, CHG) C =====================================================================| C This routine calculates the conventional Born coefficient, omega | C sub j, and its derivatives. References are to Shock, Helgeson, | C and Sverjensky, 1989. | C | C BIGGAM is big gamma defined by eqn. 19 | C CHG is the charge of the ion | C ETA is a combined constant in angstroms calories per mole | C LITG is a solvent function which varies with temperature and | C density. see Tanger and Helgeson, 1988 appendix C and H | C OMEGA is the conventional Born coefficient, defined by eqn. 16 | C WABS is the absolute Born coefficient, defined by eqn. 17 | C X1 is the first derivative of omega B-10 T & H | C X2 is the second derivative of omega B-11 T & H | C XRAD is the crystallographic radius of the ion | C ---------------------------------------------------------------------| INTEGER*2 CHG REAL*8 BIGGAM, ETA, LITG REAL*8 WABS, X1, X2, XRAD INTRINSIC FLOAT, IABS C =====================================================================| ETA = 1.66027D+05 IF (CHG .GT. 0) THEN BIGGAM = LITG + 0.94D+00 ELSE BIGGAM = LITG END IF WABS = CHG**2 / (XRAD + FLOAT(IABS(CHG)) * BIGGAM) WABS = WABS - CHG / (3.082D+00 + LITG) WABS = WABS * ETA X1 = IABS(CHG**3) / (XRAD + FLOAT(IABS(CHG)) * BIGGAM)**2 X1 = X1 - CHG / (3.082D+00 + LITG)**2 X1 = -ETA * X1 X2 = CHG**4 / (XRAD + FLOAT(IABS(CHG)) * BIGGAM)**3 X2 = X2 - CHG / (3.082D+00 + LITG)**3 X2 = 2 * ETA * X2 RETURN END SUBROUTINE GASES C =====================================================================| C | C GASES calculates standard state thermodynamic properties | C for gases at given temperature and pressure. | C The standard state is T and one bar for the ideal gas. | C | C TREF = Reference temperature (298.15 K) | C PREF = Reference pressure (one bar) | C T(J) = Index on temperature | C | C AG = | C BG = Maier-Kelly heat capacity coefficients for gases | C CG = | C | C GG = Standard state apparent free energy for a gas at | C T and P | C SG = Third law entropy for a gas at T and P | C M = Index on the gas | C | C SGREF = Standard state third law entropy at tref and pref | C GGREF = Standard state free energy of formation at tref | C and pref | C | C ---------------------------------------------------------------------| INTEGER*2 NG, NI, PM, TM PARAMETER (NG = 10, NI = 40, PM = 5, TM = 24) INTEGER*2 I, INDP(PM), INDT(TM) INTEGER*2 INDWAT, INP1, INP4, INP5, IOUT1, IOUT2, IOUT3 INTEGER*2 ITOTP, ITOTT, ITTY INTEGER*2 J, M, MAXG, NJ, NSGAS, NSION, NSMIN, PMAX INTEGER*2 SPEC(NI), SPECG(NG), TMAX REAL*8 ACTI(10), ACTM(8), ACTW, AG(NG) REAL*8 BB, BG(NG) REAL*8 CC, CG(NG), DD REAL*8 G(NI,TM,PM), G0, G1, GG(10,21), GGREF(NG), GM(8,21,10) REAL*8 GT, GW(TM,PM), GWREF REAL*8 HGREF(NG) REAL*8 P(TM,PM), PREF REAL*8 SGREF(NG) REAL*8 T(TM), TC(TM), TREF REAL*8 XMOLG(8), XRG(8), XRI(10), XRM(8), XRW CHARACTER * 8 GAS(NG,2), ION(NI), MINAM(8) EXTERNAL GT COMMON /DATGAS/ SGREF, GGREF, AG, BG, CG, HGREF COMMON /GIBBS / G, GG, GM, GW, GWREF COMMON /IO / INP1, INP4, INP5, IOUT1, IOUT2, & IOUT3, ITTY COMMON /MAX / TMAX, PMAX, INDP, INDT, ITOTT, ITOTP, MAXG COMMON /NOMS / GAS, ION, MINAM COMMON /STOCH / ACTI, ACTM, XMOLG, XRI, XRG, & XRM, XRW, ACTW, NSION, NSGAS, NSMIN, INDWAT, & SPEC, SPECG COMMON /STUFF / P, PREF, T, TC, TREF C =====================================================================| C K=INDP(NK) DD = 0.0D+0 WRITE (IOUT1, 104) WRITE (IOUT1, 106) (GAS(SPECG(M),1), M = 1, NSGAS) DO 10 NJ = 1, ITOTT J = INDT(NJ) DO 12 I = 1, NSGAS M = SPECG(I) BB = BG(M) * 0.001 CC = CG(M) * 100000.0 G1 = 0.0 G0 = GGREF(M) - SGREF(M) * (T(J) - TREF) IF (T(J) .NE. TREF) & G1 = GT (T(J), TREF, T(J), AG(M), BB, CC, DD, DD, DD, DD) GG(M,J) = G0 + G1 12 CONTINUE WRITE (IOUT1, 110) TC(J), (GG(SPECG(I),J), I = 1, NSGAS) 10 CONTINUE 104 FORMAT (///5X, 'Free energies for gases are calculated at ', & ' one bar (standard state).') 106 FORMAT (5X, 'Apparent free energies of gases at given T:' & /5X, 42('-'), /3X, 'TEMP DEG.C', 5X, A8, 7(6X, A8)) 110 FORMAT (5X, F6.2, 5X, 8(F12.0, 2X)) RETURN END SUBROUTINE REAC (XLOGK) C =====================================================================| C This subroutine computes log k's and log q's for the reactions. | C It is called only after all the necessary free energies have | C been calculated. | C ---------------------------------------------------------------------| INTEGER*2 NG, NI, PM, TM PARAMETER (NG = 10, NI = 40, PM = 5, TM = 24) INTEGER*2 I, INDP(PM), INDT(TM) INTEGER*2 INDWAT, INP1, INP4, INP5, IOUT1, IOUT2, IOUT3 INTEGER*2 ITOTP, ITOTT, ITTY INTEGER*2 J, K, MAXG, NGAS, NJ, NK, NS, NSGAS, NSION, NSMIN, PMAX INTEGER*2 SPEC(NI), SPECG(NG), TMAX REAL*8 ACTI(10), ACTM(8), ACTW REAL*8 DELA, DELG, DELGQ REAL*8 G(NI,TM,PM), GG(10,21), GM(8,21,10) REAL*8 GW(TM,PM), GWREF REAL*8 P(TM,PM), PREF, R REAL*8 T(TM), TC(TM), TREF REAL*8 XLOGK(24), XLOGQ, XMOLG(8), XRG(8), XRI(10), XRM(8) REAL*8 XRW CHARACTER* 8 GAS(NG,2), ION(NI), MINAM(8) COMMON /GIBBS / G, GG, GM, GW, GWREF COMMON /IO / INP1, INP4, INP5, IOUT1, IOUT2, & IOUT3, ITTY COMMON /MAX / TMAX, PMAX, INDP, INDT, ITOTT, ITOTP, MAXG COMMON /NOMS / GAS, ION, MINAM COMMON /STOCH / ACTI, ACTM, XMOLG, XRI, XRG, & XRM, XRW, ACTW, NSION, NSGAS, NSMIN, INDWAT, & SPEC, SPECG COMMON /STUFF / P, PREF, T, TC, TREF 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, Log K and Log Q for the given reaction | C ---------------------------------------------------------------------| C Writes first the stochiometry on output | C ---------------------------------------------------------------------| 100 WRITE (IOUT1, 750) 750 FORMAT (/1X, 'Reaction stochiometry:', & /1X, 21('-'), /20X, '(ACT.)') IF (NSION .EQ. 0) GO TO 205 DO 200 I = 1, NSION IF (SPEC(I) .NE. INDWAT) GO TO 202 WRITE (IOUT1, 617) XRW, ACTW GO TO 200 202 WRITE (IOUT1, 615) XRI(I), ION(SPEC(I)), ACTI(I) 200 CONTINUE 205 CONTINUE IF (NSGAS .EQ. 0) GO TO 215 DO 210 I = 1, NSGAS WRITE (IOUT1, 615) XRG(I), GAS(SPECG(I),1), XMOLG(I) 210 CONTINUE 215 CONTINUE IF (NSMIN .EQ. 0) GO TO 225 DO 220 I = 1, NSMIN WRITE (IOUT1, 615) XRM(I), MINAM(I), ACTM(I) 220 CONTINUE 225 CONTINUE R = 1.98726 DO 400 NK = 1, ITOTP K = INDP(NK) WRITE (IOUT1, 752) DO 410 NJ = 1, ITOTT J = INDT(NJ) DELG = 0.0 DELA = 0.0 IF (NSION .EQ. 0) GO TO 305 DO 300 I = 1, NSION IF (ACTI(I) .EQ. 0.0D0) ACTI(I) = 1.0D0 NS = SPEC(I) IF (SPEC(I) .NE. INDWAT) GO TO 302 DELG = DELG + XRW * GW(J,K) DELA = DELA + R * T(J) * ALOG(ACTW) * XRW GO TO 300 302 CONTINUE DELG = DELG + XRI(I) * G(NS,J,K) DELA = DELA + R * T(J) * ALOG(ACTI(I)) * XRI(I) 300 CONTINUE 305 CONTINUE IF (NSGAS .EQ. 0) GO TO 315 DO 310 I = 1, NSGAS NGAS = SPECG(I) DELA = DELA + R * T(J) * ALOG(XMOLG(I)) * XRG(I) DELG = DELG + XRG(I) * GG(NGAS,J) 310 CONTINUE 315 CONTINUE IF (NSMIN .EQ. 0) GO TO 325 DO 320 I = 1, NSMIN DELG = DELG + XRM(I) * GM(I,J,K) DELA = DELA + R * T(J) * ALOG(ACTM(I)) * XRM(I) 320 CONTINUE 325 CONTINUE XLOGK(NJ) = -DELG / (2.302585 * R * T(J)) DELGQ = DELG + DELA XLOGQ = -DELGQ / (2.302585 * R * T(J)) WRITE (IOUT1, 625) TC(J), P(J,K), DELG, XLOGK(NJ), & DELGQ, XLOGQ 410 CONTINUE 400 CONTINUE RETURN 615 FORMAT (1X, F7.3, 1X, ': ', A8, 2X, F4.2) 617 FORMAT (1X, F7.3, ' : H2O ', F4.2) 620 FORMAT (/) 625 FORMAT (1X, F5.0, 2X, F9.4, 2X, 2(E12.5, 2X, F8.3, 2X)) 630 FORMAT (1X, 'Apparent free energies of minerals at given P and T:' & /1X, 50('-'),/1X, 'Deg.C', 4X, 'Bars', 5X, A8, 7(6X, A8)) 635 FORMAT (1X, F5.0, 1X, F9.4, 1X, 8(F12.0, 2X)) 752 FORMAT (/1X, 'Deg.C', 4X, 'Bars', 5X, 'Delta G Reac', & 3X, 'Log K', 4X, 'Delta G Reac', 3X, 'Log Q', & /1X, '-----', 2X, '---------', 2X, '---ACT = 1--', & 2X, '--------', 2X, '------------', 2X, '--------') END SUBROUTINE README (MIN) C =====================================================================| C README reads the mineral data file and makes some necessary | C conversions. | C ---------------------------------------------------------------------| INTEGER*2 I, INP1, INP4, INP5 INTEGER*2 IOUT1, IOUT2, IOUT3, ITTY, J CHARACTER*8 JOU CHARACTER*8 MIN(2) REAL*8 GRF, HCT(0:2,7), HRF, HTRAN1, HTRAN2, SRF REAL*8 SLOPE1, SLOPE2, TTRAN1, TTRAN2, TTRAN3 REAL*8 VCH1, VCH2, VRE1, VRE2, VREF COMMON /DATMIN/ GRF, HRF, SRF, VREF, HCT, & TTRAN1, HTRAN1, VRE1, SLOPE1, & TTRAN2, HTRAN2, VRE2, SLOPE2, & TTRAN3 COMMON /IO / INP1, INP4, INP5, IOUT1, IOUT2, & IOUT3, ITTY C =====================================================================| READ (INP4, 500) JOU, GRF, HRF, SRF, VREF, & HCT(0,1), HCT(0,2), HCT(0,3), & HCT(0,4), HCT(0,5), HCT(0,6), HCT(0,7), & TTRAN1, HTRAN1, VCH1, SLOPE1, & HCT(1,1), HCT(1,2), HCT(1,3), & HCT(1,4), HCT(1,5), HCT(1,6), HCT(1,7), & TTRAN2, HTRAN2, VCH2, SLOPE2, & HCT(2,1), HCT(2,2), HCT(2,3), & HCT(2,4), HCT(2,5), HCT(2,6), HCT(2,7), TTRAN3 WRITE (IOUT1, 501) MIN(1), GRF, HRF, SRF, VREF IF (TTRAN1 .EQ. 0.0) GO TO 10 WRITE (IOUT1, 502) WRITE (IOUT1, 503) TTRAN1, HTRAN1, SLOPE1, VCH1, & HCT(1,1), HCT(1,2), HCT(1,3), & HCT(1,4), HCT(1,5), HCT(1,6), HCT(1,7) IF (TTRAN2 .EQ. 0.0) GO TO 10 WRITE (IOUT1, 503) TTRAN2, HTRAN2, SLOPE2, VCH2, & HCT(2,1), HCT(2,2), HCT(2,3), & HCT(2,4), HCT(2,5), HCT(2,6), HCT(2,7) 10 CONTINUE IF (JOU .EQ. 'J' .OR. JOU .EQ. 'j') THEN GRF = GRF / 4.184 HRF = HRF / 4.184 SRF = SRF / 4.184 HTRAN1 = HTRAN1 / 4.184 HTRAN2 = HTRAN2 / 4.184 DO 30 I = 1,7 DO 20 J = 0,2 HCT(J,I) = HCT(J,I) / 4.184 20 CONTINUE 30 CONTINUE END IF IF (SLOPE1 .EQ. 0.0) SLOPE1 = 100000.0D+00 IF (SLOPE2 .EQ. 0.0) SLOPE2 = 100000.0D+00 C CT1 = TTRAN1 - 273.15 C CT2 = TTRAN2 - 273.15 C CT3 = TTRAN3 - 273.15 DO 40 I = 0, 2 HCT(I,2) = HCT(I,2) * 0.001 HCT(I,3) = HCT(I,3) * 100000.0 HCT(I,4) = HCT(I,4) * 100.0 HCT(I,5) = HCT(I,5) * 1.0E-5 HCT(I,6) = HCT(I,6) * 1.0E+07 HCT(I,7) = HCT(I,7) * 1.0E+04 40 CONTINUE VRE1 = VCH1 + VREF VRE2 = VCH2 + VRE1 VREF = VREF * 0.0239 VRE1 = VRE1 * 0.0239 VRE2 = VRE2 * 0.0239 IF (VRE1 .NE. VREF .AND. HTRAN1 .NE. 0.0) SLOPE1 = HTRAN1 / & ((VRE1 - VREF) * TTRAN1) IF (VRE2 .NE. VRE1 .AND. HTRAN2 .NE. 0.0) SLOPE2 = HTRAN2 / & ((VRE2 - VRE1) * TTRAN2) RETURN 500 FORMAT (A1, 2F15.0, 2F15.5 / 3(6F12.4 / 4F12.4 /)) 501 FORMAT (/1X, A8, 2X, F12.0, 4X, F12.0, 8X, F6.2, 26X, F10.3) 502 FORMAT (11X, 'Transition T', 4X, 'H trans', 3X, 'Exp dP/dT', & 2X, 'V change') 503 FORMAT (6X, '-----', 3X, F7.2, ' K', 2X, F12.0, 2X, F6.2, & 2X, F10.3, 2X, 6(F10.3, 2X), 2X, E12.6) END SUBROUTINE SOLIDS (IND) C =====================================================================| C SOLIDS calculates thermodynamic data at P and T for mineral | C phases using Maier-Kelly heat capacity coefficients and constant | C volume (t,pref). | C The standard state is T and P for the pure solid. | C | C TREF = Reference temperature (298.15 K) | C PREF = Reference pressure (one bar) | C T(J) = Temperature (subscript) | C P(J,K)= Pressure (subscript) | C | C VREF = Mineral molar volume at 298.15 K and one bar | C in units of cal/bar/mole. | C GRF = Standard state free energy of formation at | C 298.15 K and one bar. | C SRF = Third law entropy at 298.15 K and one bar. | C MIN = Phase name | C GM = Standard state apparent free energy of the phase | C at T and P | C | C G is obtained from double CP integration over temprature | C range then integration over pressure range with volume. | C The Maier-Kelley power function coefficients are for the | C following equation | C 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 | C done from pref to ppp (where ppp is a point in between pp1 and | C pp2) with vref+vch1+vch2, then from ppp to P(k) with vref only. | C | C ---------------------------------------------------------------------| INTEGER*2 NG, NI, PM, TM PARAMETER (NG = 10, NI = 40, PM = 5, TM = 24) INTEGER*2 IND, INDP(PM), INDT(TM), INP1, INP4, INP5 INTEGER*2 IOUT1, IOUT2, IOUT3, ITOTP, ITOTT, ITTY INTEGER*2 I, J, K, MAXG, NJ, NK, PMAX, TMAX REAL*8 G(NI,TM,PM), GG(10,21), GI(5), GM(8,21,10) REAL*8 GP, GRF, GT, GW(TM,PM), GWREF REAL*8 HCT(0:2,7), HRF, HTRAN1, HTRAN2 REAL*8 P(TM,PM), PP1, PP2, PPP, PREF REAL*8 SRF, SLOPE1, SLOPE2 REAL*8 T(TM), TC(TM), TREF, TTRAN1, TTRAN2, TTRAN3 REAL*8 VRE1, VRE2, VREF EXTERNAL GP, GT COMMON /DATMIN/ GRF, HRF, SRF, VREF, HCT, & TTRAN1, HTRAN1, VRE1, SLOPE1, & TTRAN2, HTRAN2, VRE2, SLOPE2, & TTRAN3 COMMON /GIBBS / G, GG, GM, GW, GWREF COMMON /IO / INP1, INP4, INP5, IOUT1, IOUT2, & IOUT3, ITTY COMMON /MAX / TMAX, PMAX, INDP, INDT, ITOTT, ITOTP, MAXG COMMON /STUFF / P, PREF, T, TC, TREF C =====================================================================| DO 230 NK = 1, ITOTP K = INDP(NK) DO 220 NJ = 1, ITOTT J = INDT(NJ) DO 100 I = 1, 5 GI(I) = 0.0D+00 100 CONTINUE IF (T(J) .EQ. TREF) GO TO 200 IF (TTRAN2 .EQ. 0.0 .OR. T(J) .LT. TTRAN2) GO TO 150 C ---------------------------------------------------------------------| C Second Transition | C ---------------------------------------------------------------------| PP1 = (T(J) - TTRAN1) * SLOPE1 + PREF PP2 = (T(J) - TTRAN2) * SLOPE2 + PREF C ---------------------------------------------------------------------| C Integrates G at T(J) over pressure range | C ---------------------------------------------------------------------| IF (PP2 .GT. PP1) THEN PPP = (PP1 + PP2) / 2.0D+0 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) ELSE IF (P(J,K) .GT. PP1) THEN GI(5) = GP (VRE2, PREF, PP2) + GP (VRE1, PP2, PP1) + & GP (VREF, PP1, P(J,K)) ELSE IF (P(J,K) .GT. PP2) THEN GI(5) = GP (VRE2, PREF, PP2) + GP (VRE1, PP2, P(J,K)) ELSE GI(5) = GP (VRE2, PREF, P(J,K)) END IF C ---------------------------------------------------------------------| C Integrates G at P(J,K) from TTRAN1 to TTRAN2 to T(J) | C ---------------------------------------------------------------------| GI(4) = GT (T(J), TTRAN2, T(J), HCT(2,1), HCT(2,2), & HCT(2,3), HCT(2,4), HCT(2,5), HCT(2,6), HCT(2,7)) GI(3) = GT (T(J), TTRAN1, TTRAN2, HCT(1,1), HCT(1,2), & HCT(1,3), HCT(1,4), HCT(1,5), HCT(1,6), HCT(1,7)) & + HTRAN2 * (1.0 - T(J) / TTRAN2) GO TO 180 C ---------------------------------------------------------------------| C First Transition | C ---------------------------------------------------------------------| 150 CONTINUE IF (TTRAN1 .EQ. 0.0 .OR. T(J) .LT. TTRAN1) GO TO 190 C ---------------------------------------------------------------------| C Integrates G at T(J) over pressure range | C ---------------------------------------------------------------------| PP1 = (T(J) - TTRAN1) * SLOPE1 + PREF IF (P(J,K) .GT. PP1) THEN GI(5) = GP (VRE1, PREF, PP1) + GP (VREF, PP1, P(J,K)) ELSE GI(5) = GP (VRE1, PREF, P(J,K)) END IF C ---------------------------------------------------------------------| C Integrates G at P(J,K) from TREF to TTRAN1 to T(J) | C ---------------------------------------------------------------------| 170 CONTINUE GI(3) = GT (T(J), TTRAN1, T(J), HCT(1,1), HCT(1,2), & HCT(1,3), HCT(1,4), HCT(1,5), HCT(1,6), HCT(1,7)) 180 CONTINUE GI(2) = GT (T(J), TREF, TTRAN1, HCT(0,1), HCT(0,2), & HCT(0,3), HCT(0,4), HCT(0,5), HCT(0,6), HCT(0,7)) & + HTRAN1 * (1.0 - T(J) / TTRAN1) GO TO 210 C ---------------------------------------------------------------------| C No Transitions | C ---------------------------------------------------------------------| 190 CONTINUE GI(2) = GT (T(J), TREF, T(J), HCT(0,1), HCT(0,2), & HCT(0,3), HCT(0,4), HCT(0,5), HCT(0,6), HCT(0,7)) 200 CONTINUE GI(5) = GP (VREF, PREF, P(J,K)) 210 CONTINUE GI(1) = GRF - SRF * (T(J) - TREF) GM(IND,J,K) = GI(1) + GI(2) + GI(3) + GI(4) + GI(5) 220 CONTINUE 230 CONTINUE RETURN 700 FORMAT (5X, 'CLAPEYRON SLOPES HAVE CROSSED P= ', F10.3, & ' T= ', F10.3) END REAL * 8 FUNCTION GP (V, P1, P2) REAL * 8 P1, P2, V GP = V * (P2 - P1) RETURN END REAL*8 FUNCTION GT (T, T1, T2, A, B, C, D, E, F, G) REAL*8 A, B, C, D, E, F, G, ENTH, ENTR, T, T1, T2 ENTH = A * (T2 - T1) & + B * 0.5 * (T2**2 - T1**2) & - C * (1.0 / T2 - 1.0 / T1) & + D * 2.0 * (DSQRT(T2) - DSQRT(T1)) & + E * (1.0 / 3.0) * (T2**3 - T1**3) & - F * 0.5 * (1.0 / T2**2 - 1.0 / T1**2) & + G * ALOG(T2/T1) ENTR = A * ALOG(T2/T1) & + B * (T2 - T1) & - C * 0.5 * (1.0 / T2**2 - 1.0 / T1**2) & - D * 2.0 * (1.0 / DSQRT(T2) - 1.0 / DSQRT(T1)) & + E * 0.5 * (T2**2 - T1**2) & - F * (1.0 / 3.0) * (1.0 / T2**3 - 1.0 / T1**3) & - G * (1.0 / T2 - 1.0 / T1) GT = ENTH - T * ENTR RETURN END SUBROUTINE WRITIO C =====================================================================| C WRITIO is used to echo the ion data (that was read in) on the | C output. | C ---------------------------------------------------------------------| INTEGER*2 NG, NI, PM, TM PARAMETER (NG = 10, NI = 40, PM = 5, TM = 24) INTEGER*2 I, INDWAT, INP1, INP4, INP5 INTEGER*2 IOUT1, IOUT2, IOUT3, ITTY INTEGER*2 N, NSGAS, NSION, NSMIN INTEGER*2 SPEC(NI), SPECG(NG), Z(NI) CHARACTER* 8 GAS(NG,2), ION(NI), MINAM(8) REAL*8 A1(NI), A2(NI), A3(NI), A4(NI) REAL*8 ACTI(10), ACTM(8), ACTW, AG(NG), BG(NG) REAL*8 C1(NI), C2(NI), CG(NG), CPREF(NI) REAL*8 G(NI,TM,PM), GG(10,21), GGREF(NG), GM(8,21,10) REAL*8 GREF(NI), GW(TM,PM), GWREF REAL*8 HGREF(NG), HREF(NI) REAL*8 SGREF(NG), SREF(NI) REAL*8 VIREF(NI), WREF(NI) REAL*8 XMOLG(8), XRG(8), XRAD(NI), XRI(10), XRM(8), XRW COMMON /DATGAS/ SGREF, GGREF, AG, BG, CG, HGREF COMMON /DATION/ A1, A2, A3, A4, C1, C2, CPREF, GREF, HREF, & SREF, VIREF, WREF, XRAD, Z COMMON /GIBBS / G, GG, GM, GW, GWREF COMMON /IO / INP1, INP4, INP5, IOUT1, IOUT2, & IOUT3, ITTY COMMON /NOMS / GAS, ION, MINAM COMMON /STOCH / ACTI, ACTM, XMOLG, XRI, XRG, & XRM, XRW, ACTW, NSION, NSGAS, NSMIN, INDWAT, & SPEC, SPECG C =====================================================================| WRITE (IOUT1, 507) DO 104 I = 1, NSION N = SPEC(I) IF (N .EQ. INDWAT) GO TO 106 WRITE (IOUT1, 510) ION(N), GREF(N), HREF(N), & SREF(N), A1(N), A2(N), & A3(N), A4(N), C1(N), & C2(N), WREF(N), XRAD(N) GO TO 104 106 WRITE (IOUT1, 500) GWREF 104 CONTINUE RETURN ENTRY WRITGA C ---------------------------------------------------------------------| C WRITGA echoes the gas data (that was read in) on the output | C ---------------------------------------------------------------------| WRITE (IOUT1, 519) DO 113 I = 1, NSGAS N = SPECG(I) WRITE (IOUT1, 522) GAS(N,1), GGREF(N), HGREF(N), SGREF(N), & AG(N), BG(N), CG(N) 113 CONTINUE RETURN 500 FORMAT (1X, 'H2O ', 2X, F12.0/) 507 FORMAT (//1X, 'Species ', 6X, 'Gref', 10X, 'Href', 7X, 'Sref', & 8X, 'A1', 12X, 'A2', 12X, 'A3', 12X, 'A4', & /1X, '--------', 2X, '------------', 2X, '------------', & 2X, '------', 2X, '-----C1-----', 2X, '-----C2-----', & 2X, '---W ref----', 2X, '---X rad----') 510 FORMAT (1X, A8, 2X, 2(F12.0, 2X), F6.2, 2X, 4(E12.4, 2X)/, & 47X, 4(E12.4, 2X)) 519 FORMAT (//6X, ' Gas ', 10X, 'Gref', 10X, 'Href', 7X, 'Sref', & 11X, 'A', 7X, 'B*10**3', 3X, 'C*10**-5', & /7X, '---', 12X, '----', 10X, '----', 7X, '----', & 11X, '-', 7X, '-------', 3X, '--------') 522 FORMAT (5X, A8, 5X, 2(F12.0, 2X), F6.2, 5X, 3(F10.3, 1X)) END SUBROUTINE READER C =====================================================================| C Reader is a routine to read mineral and gas stochiometries in | C the file soltherm. These given stochiometries are used to | C calculate log k's from the thermo data contained in the ion, min | C and gas data files. The user can choose to calculate new | C log k's for a set of either given mineral names or given mineral | C compositions, or to recalculate log k's for the whole set of | C gases and minerals in SOLTHERM. In each case, the calculated | C log k's are compared with the original log k's from SOLTHERM and | C the differences between the two sets of values are computed. | C The routine also creates SOLTHERM formatted output with the | C calculated new values. | C | C ** Revised May 87 for compatibility w/ new SOLTHERM N.S. | C ** Revised Jul 87 to inlude "fake" complexes as minerals | C ---------------------------------------------------------------------| INTEGER*2 NG, NI, NM, PM, TM PARAMETER (NG = 10, NI = 40, NM = 400, PM = 5, TM = 24) INTEGER*2 COMTOT, GASTOT INTEGER*2 I, I2, I3, IEQ(NI), IFL, IFLAG(NM), IGAS INTEGER*2 IMATCH, IMIN, IND(NM,11), INDEX, INDP(PM) INTEGER*2 INDT(TM), INDWAT, INP1, INP4, INP5 INTEGER*2 IOUT1, IOUT2, IOUT3, IREW INTEGER*2 ISEL, IT, IT1, ITOT(NM), ITOTP, ITOTT, ITTY INTEGER*2 J, K, KKK, KSPEC(8), KTOT, MAXG, MINTOT INTEGER*2 N, NGAS, NNAM, NS, NSGAS, NSION, NSMIN, NSPEC INTEGER*2 PLXTOT, PMAX, SAQ(NI), SPEC(NI), SPECG(NG) INTEGER*2 SS2(NM), TMAX, Z(NI) REAL*8 A1(NI), A2(NI), A3(NI), A4(NI) REAL*8 ACTI(10), ACTM(8), ACTW, AKLOG(NM,9) REAL*8 C1(NI), C2(NI), COEF(NM,11), CPREF(NI) REAL*8 DELK(8), DENS(NM), DGDP(TM,PM), DGDT(TM,PM) REAL*8 DG2DT2(TM,PM), DIEL(TM,PM) REAL*8 G(NI,TM,PM), GG(10,21), GM(8,21,10) REAL*8 GREF(NI), GW(TM,PM), GWREF, HREF(NI), LITG(TM,PM) REAL*8 P(TM,PM), PREF, Q(TM,PM), SREF(NI), SS1(NM) REAL*8 T(TM), TC(TM), TREF REAL*8 VIREF(NI), WREF(NI), WT(NM) REAL*8 X(TM,PM), XLOGK(24), XMOLG(8), XRAD(NI) REAL*8 XRG(8), XRI(10), XRM(8), XRW, Y(TM,PM) CHARACTER*1 ANS, SOLANS CHARACTER*8 BLANK, C123, CHOICE(2), DUM(25,10), DUMMY CHARACTER*8 GAS(NG,2), HION, ION(NI), IONAM(NI) CHARACTER*8 MIN(2), MINAM(8), NAME(NM), NOM(10), REF(NM) EXTERNAL AQUION, GASES, REAC, README, SOLIDS, WRITGA, WRITIO COMMON /DATION/ A1, A2, A3, A4, C1, C2, CPREF, GREF, HREF, & SREF, VIREF, WREF, XRAD, Z COMMON /GIBBS / G, GG, GM, GW, GWREF COMMON /IO / INP1, INP4, INP5, IOUT1, IOUT2, & IOUT3, ITTY COMMON /MAX / TMAX, PMAX, INDP, INDT, ITOTT, ITOTP, MAXG COMMON /NOMS / GAS, ION, MINAM COMMON /STOCH / ACTI, ACTM, XMOLG, XRI, XRG, & XRM, XRW, ACTW, NSION, NSGAS, NSMIN, INDWAT, & SPEC, SPECG COMMON /STUFF / P, PREF, T, TC, TREF COMMON /WATER / DIEL, LITG, Q, X, Y, DGDT, DGDP, DG2DT2 C ---------------------------------------------------------------------| C SAQ are the indexes of aq. species in SOLTHERM (as in ionam) | C IEQ are the indexes for the same species in ion data file | C IEQ is the conversion array between the numbers of the ions in | C SOLTHERM and in the ion data file. Updated SEP 88, in mostly | C same order. Watch that Au and Hg are actualy chlorides in | C SOLTHERM, therefore, SOLTHERM Au and Hg soichiometries should | C be reexpressed w/o chlorides. | C ---------------------------------------------------------------------| DATA IEQ / 1, 40, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, & 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, & 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, & 40/ DATA SAQ / 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, & 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, & 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, & 40/ DATA BLANK /' '/ DATA C123 /'C123 COE'/ DATA CHOICE /'Gases ', 'Minerals'/ DATA HION /'H+ '/ C =====================================================================| C The following reads the component species in SOLTHERM data, | C then skips until the beginning of the gas data | C first, we skip until H+ is found | C ---------------------------------------------------------------------| ISEL = 0 IMIN = 0 READ (INP5, '(4I5)') COMTOT, PLXTOT, GASTOT, MINTOT IF (COMTOT .EQ. 0) THEN WRITE (ITTY, 5000) READ (ITTY, '(A8)') DUMMY IF (DUMMY .EQ. BLANK) THEN COMTOT = 37 ELSE READ (DUMMY, '(I5)') COMTOT END IF END IF IF (MINTOT .EQ. 0 .OR. GASTOT .EQ. 0) THEN KTOT = NM ELSE IF (MINTOT+GASTOT .LE. NM) THEN KTOT = MINTOT + GASTOT ELSE PRINT *, 'The number of minerals + gases in SOLTHERM exceed ' PRINT *, ' the max. array size of ', NM STOP END IF 100 CONTINUE READ (INP5, '(A8)', END = 580) IONAM(1) IF (IONAM(1) .NE. HION) GO TO 100 DO 110 I = 2, COMTOT READ (INP5, '(A8)') IONAM(I) 110 CONTINUE C ---------------------------------------------------------------------| C Now we look for the C123 entry after the complexes | C ---------------------------------------------------------------------| 120 CONTINUE READ (INP5, '(A8)', END = 580) DUMMY IF (DUMMY .NE. C123) GO TO 120 C ---------------------------------------------------------------------| C Now we should be at the top of the gas data. | C GASTOT is the total number of gases in SOLTHERM | C MAXG is the total of gases in gas data (transfered from main) | C COMTOT is the total number of component species in SOLTHERM | C KTOT is the maximum number of gases + minerals in SOLTHERM | C ---------------------------------------------------------------------| WRITE (ITTY, 5020) KTOT IF (GASTOT .EQ. 0) THEN WRITE (ITTY, 5030) READ (ITTY, '(A1)') DUMMY IF (DUMMY .EQ. BLANK) THEN GASTOT = 20 ELSE READ (DUMMY, '(I5)') GASTOT END IF END IF IF (GASTOT .GT. KTOT) GASTOT = KTOT WRITE (ITTY, 5040) READ (ITTY, '(A1)') SOLANS C ---------------------------------------------------------------------| C Now, there is a trick to calculate data for complexes. The | C derived species must be included in the mineral data file, | C and it will be treated as a mineral (no use of subroutine | C AQUION). In SOLTHERM, all the derived aq. species of interest | C must be placed after the gases, and only the 'all minerals' | C option will be used. Suppress all minerals in SOLTHERM. | C ---------------------------------------------------------------------| WRITE (ITTY, 5050) READ (ITTY, *) IMIN IF (IMIN .EQ. 3) GO TO 300 125 CONTINUE WRITE (ITTY, 5060) CHOICE(IMIN), CHOICE(IMIN) READ (ITTY, *) ISEL IF (ISEL .LT.1 .OR. ISEL .GT. 3) GO TO 125 IF (ISEL .EQ. 2) GO TO 200 IF (ISEL .EQ. 3) GO TO 300 C ---------------------------------------------------------------------| C Selection by name. | C Names in both SOLTHERM and Mineral.dat should match | C ---------------------------------------------------------------------| WRITE (ITTY, 5070) CHOICE(IMIN) DO 130 I = 1, 10 READ (ITTY, '(A8)') NOM(I) IF (NOM(I) .EQ. BLANK) GO TO 140 130 CONTINUE 140 CONTINUE NNAM = I - 1 C ---------------------------------------------------------------------| C Match selected names with SOLTHERM and get STOCHIO | C ---------------------------------------------------------------------| NS = 1 DO 150 I = 1, NM IFLAG(I) = 0 150 CONTINUE IF (IMIN .EQ. 1) KTOT = GASTOT DO 170 K = 1, KTOT READ (INP5, 5080, END = 180) NAME(NS), WT(NS), IT, & (COEF(NS,I), IND(NS,I), I = 1, IT), & SS1(NS), SS2(NS) ITOT(NS) = IT DO 160 I = 1, NNAM IF (NOM(I) .NE. NAME(NS)) GO TO 160 WRITE (ITTY, 5090) NAME(NS) READ (INP5, 5100, END = 180) (AKLOG(NS,J), J = 1, 8), . REF(NS) READ (INP5, 5110) DENS(NS) IF (SS2(NS) .EQ. 6) READ (INP5, 5120) (DUM(NS,J), J = 1, 10) IF (NS .EQ. NNAM) GO TO 380 C ---------------------------------------------------------------------| C NS incremented only when names match | C ---------------------------------------------------------------------| NS = NS + 1 IFLAG(I) = K GO TO 170 160 CONTINUE C ---------------------------------------------------------------------| C Skips log k's and regression coefficients (2 records) | C and gas virial coefficients if SS2(NS) = 6 | C ---------------------------------------------------------------------| READ (INP5, '(/)', END = 180) IF (SS2(NS) .EQ. 6) READ (INP5, '(A)') NAME(NS) 170 CONTINUE C ---------------------------------------------------------------------| C What mineral(s) was (were) not found in SOLTHERM | C ---------------------------------------------------------------------| 180 CONTINUE NS = 0 DO 190 I = 1, NNAM IF (IFLAG(I) .EQ. 0) THEN WRITE (ITTY, 5130) NOM(I) ELSE NS = NS + 1 END IF 190 CONTINUE REWIND INP5 GO TO 120 C ---------------------------------------------------------------------| C Selection by composition. The compostion is given by the | C species indexes from SOLTHERM (SAQ) | C ---------------------------------------------------------------------| 200 CONTINUE WRITE (ITTY, 5140) DO 210 I = 1, 10 I2 = I + 10 I3 = I + 20 WRITE (ITTY, 5150) I, IONAM(I), I2, IONAM(I2), I3, IONAM(I3) 210 CONTINUE DO 220 I = 1, 8 KSPEC(I) = 0 220 CONTINUE READ (ITTY, *) (KSPEC(I), I = 1, 8) DO 230 I = 1, 8 IF (KSPEC(I) .EQ. 0) GO TO 240 230 CONTINUE 240 CONTINUE NSPEC = I - 1 C ---------------------------------------------------------------------| C Read SOLTHERM and get names of gases or minerals | C which contains given species | C ---------------------------------------------------------------------| IF (IMIN .EQ. 1) KTOT = GASTOT NS = 1 DO 280 K = 1, KTOT READ (INP5, 5080, END = 290) NAME(NS), WT(NS), IT, (COEF(NS,I), & IND(NS,I), I = 1, IT) ITOT(NS) = IT IMATCH = 0 DO 260 N = 1, NSPEC DO 250 I = 1, IT IF (IND(NS,I) .EQ. KSPEC(N)) IMATCH = IMATCH + 1 IF (IMATCH .EQ. NSPEC) GO TO 270 250 CONTINUE 260 CONTINUE C ---------------------------------------------------------------------| C Skip log k data and regression coefficients (2 records) | C and gas virial coefficients if SS2(NS) = 6 | C ---------------------------------------------------------------------| READ (INP5, '(/)', END = 290) IF (SS2(NS) .EQ. 6) READ (INP5, '(/)') GO TO 280 270 CONTINUE WRITE (ITTY, 5170) NAME(NS) READ (INP5, 5100, END = 290) (AKLOG(NS,J), J = 1, 8), REF(NS) READ (INP5, 5110) DENS(NS) IF (SS2(NS) .EQ. 6) READ (INP5, 5120) (DUM(NS,J), J = 1, 10) C ---------------------------------------------------------------------| C NS incremented only when stochiometries match | C ---------------------------------------------------------------------| NS = NS + 1 280 CONTINUE 290 CONTINUE NNAM = NS - 1 IF (NNAM .GT. 0) GO TO 380 WRITE (ITTY, 5160) CHOICE(IMIN) GO TO 570 C ---------------------------------------------------------------------| C Get all the minerals or gases from SOLTHERM | C ---------------------------------------------------------------------| 300 CONTINUE IFL= 0 C ---------------------------------------------------------------------| C Reads all the gases | C ---------------------------------------------------------------------| K = GASTOT 310 CONTINUE DO 320 NS = 1, K READ (INP5, 5080, END = 330) NAME(NS), WT(NS), IT, & (COEF(NS,I), IND(NS,I), I = 1, IT), SS1(NS), SS2(NS) ITOT(NS) = IT READ (INP5, 5100, END = 330) (AKLOG(NS,J), J = 1, 8), REF(NS) READ (INP5, 5110) DENS(NS) IF (SS2(NS) .EQ. 6) READ (INP5, 5120) (DUM(NS,J), J = 1, 10) 320 CONTINUE 330 CONTINUE NNAM = GASTOT IF (IMIN .EQ. 1) GO TO 380 IF (IMIN .EQ. 3) GO TO 350 IF (IFL .EQ. 0) GO TO 340 NNAM = NS - 1 GO TO 380 C ---------------------------------------------------------------------| C Reads all the minerals | C ---------------------------------------------------------------------| 340 CONTINUE IFL = 1 K = MINTOT GO TO 310 C ---------------------------------------------------------------------| C Reads list of complexes placed after the gases in SOLTHERM | C that will be treated like minerals | C ---------------------------------------------------------------------| 350 CONTINUE DO 360 NS = 1, KTOT READ (INP5, 815, END = 370) NAME(NS), WT(NS), IT, & (COEF(NS,I), IND(NS,I), I = 1, IT) IF (IT .EQ. 0) GO TO 370 ITOT(NS) = IT - 1 IF (IT .EQ. 0) GO TO 370 READ (INP5, 5100) (AKLOG(NS,J), J = 1, 8), REF(NS) READ (INP5, 5110) DENS(NS) 360 CONTINUE 370 CONTINUE NNAM = NS - 1 380 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 390 I = 1, NM IFLAG(I) = 0 390 CONTINUE IF (IMIN .GT. 1) GO TO 460 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 430 NGAS = 1, MAXG DO 420 NS = 1, NNAM IF (NAME(NS) .NE. GAS(NGAS,1)) GO TO 420 WRITE (IOUT1, '(A1)') '1' IFLAG(NS) = NS IGAS = IGAS + 1 NSION = ITOT(NS) DO 400 I = 1, NSION XRI(I) = COEF(NS,I) KKK = IND(NS,I) SPEC(I) = IEQ(KKK) IF (SPEC(I) .NE. INDWAT) GO TO 400 XRW = XRI(I) ACTW = 1.0 ACTI(I) = 1.0 400 CONTINUE CALL WRITIO CALL AQUION NSMIN = 0 NSGAS = 1 XMOLG(1) = 1 SPECG(1) = NGAS XRG(1) = -1.0 CALL WRITGA CALL GASES CALL REAC (XLOGK) C ---------------------------------------------------------------------| C Writes SOLTHERM formatted file with new log k's and displays | C the difference between old and new log k's values | C ---------------------------------------------------------------------| IF (SOLANS .EQ. 'Y' .OR. SOLANS .EQ. 'y') THEN IT = ITOT(NS) WRITE (IOUT2, 5080) NAME(NS), WT(NS), IT, & (COEF(NS,I), IND(NS,I), I = 1, IT), SS2(NS) WRITE (IOUT2, 826) NAME(NS), (XLOGK(J), J = 1, 8), REF(NS) WRITE (IOUT2, 823) NAME(NS), DENS(NS) IF (SS2(NS) .EQ. 6) WRITE (IOUT2, 5120) & (DUM(NS,J), J = 1, 10) DO 410 J = 1, 8 DELK(J) = XLOGK(J) - AKLOG(NS,J) 410 CONTINUE WRITE (IOUT3, 826) NAME(NS), (DELK(J), J = 1, 8) END IF IF (IGAS .NE. NNAM) GO TO 430 GO TO 450 420 CONTINUE 430 CONTINUE C ---------------------------------------------------------------------| C What gas was not found ?? | C ---------------------------------------------------------------------| DO 440 NS = 1, NNAM IF (IFLAG(NS) .EQ. 0) WRITE (ITTY, 5180) NAME(NS) 440 CONTINUE 450 CONTINUE WRITE (ITTY, 5190) WRITE (ITTY, 5200) READ (ITTY, '(A1)') ANS IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') RETURN REWIND INP5 GO TO 120 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 SOLTHERM 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 SOLTHERM. | C ---------------------------------------------------------------------| 460 CONTINUE C ---------------------------------------------------------------------| C If all SOLTHERM minerals have to be computed, or if more | C than 5 minerals are given, all the ions are calculated first | C ---------------------------------------------------------------------| IF (NNAM .LE. 5) GO TO 480 NSION = COMTOT DO 470 I = 1, NSION SPEC(I) = IEQ(SAQ(I)) 470 CONTINUE CALL WRITIO CALL AQUION 480 CONTINUE NS = 1 IREW = 0 485 CONTINUE READ (INP4, '(A8)') MIN(1) IF (MIN(1) .NE. '********') GO TO 485 490 CONTINUE READ (INP4, '(A8)', END = 560) MIN(1) IF (NAME(NS) .EQ. MIN(1)) GO TO 500 READ (INP4, '(//////)', END = 560) GO TO 490 500 CONTINUE WRITE (ITTY, 5210) MIN(1) WRITE (IOUT1, 5230) CALL README(MIN) IFLAG(NS) = NS NSION = ITOT(NS) DO 510 I = 1, NSION XRI(I) = COEF(NS,I) KKK = IND(NS,I) SPEC(I) = IEQ(KKK) IF (SPEC(I) .NE. INDWAT) GO TO 510 XRW = XRI(I) ACTW = 1.0 ACTI(I) = 1.0 510 CONTINUE IF (NNAM .GT. 5) GO TO 520 CALL WRITIO CALL AQUION 520 CONTINUE NSMIN = 1 NSGAS = 0 XRM(1) = -1.0 ACTM(1) = 1.0 MINAM(1) = NAME(NS) INDEX = 1 CALL SOLIDS (INDEX) CALL REAC (XLOGK) C ---------------------------------------------------------------------| C Writes SOLTHERM formatted file with new log k's and displays | C the difference between old and new log k's values. (Kind of | C sloppy to repeat everything here, but c'est la vie !) | C ---------------------------------------------------------------------| IF (SOLANS .EQ. 'Y' .OR. SOLANS .EQ. 'y') THEN IT = ITOT(NS) IF (IMIN .EQ. 3) GO TO 530 WRITE (IOUT2, 5080) NAME(NS), WT(NS), IT, (COEF(NS,I), & IND(NS,I), I = 1, IT) GO TO 540 C ---------------------------------------------------------------------| C Writes complexe style stoichiometry | C ---------------------------------------------------------------------| 530 CONTINUE C IT1 = IT+1 WRITE (IOUT2, 815) NAME(NS), WT(NS), IT, (COEF(NS,I), & IND(NS,I), I = 1, IT) 540 CONTINUE WRITE (IOUT2, 826) NAME(NS), (XLOGK(J), J = 1, 8), REF(NS) WRITE (IOUT2, 823) NAME(NS), DENS(NS) IF (SS2(NS). EQ. 6) WRITE (IOUT2, 5120) (DUM(NS,J), J = 1, 10) DO 550 J = 1, 8 DELK(J) = XLOGK(J) - AKLOG(NS,J) 550 CONTINUE WRITE (IOUT3, 826) NAME(NS), (DELK(J), J = 1, 8) END IF NS = NS + 1 IREW = 0 IF (NS .GT. NNAM) GO TO 570 GO TO 490 C ---------------------------------------------------------------------| C If a mineral is not found after the first 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. | C ---------------------------------------------------------------------| 560 CONTINUE REWIND INP4 IREW = IREW + 1 IF (IREW .LT. 2) GO TO 485 IREW = 0 WRITE (ITTY, 5220) NAME(NS) NS = NS + 1 IF (NS .LE. NNAM) GO TO 490 570 CONTINUE WRITE (ITTY, 5190) WRITE (ITTY, 5200) READ (ITTY, '(A1)') ANS IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') THEN REWIND INP5 GO TO 120 END IF RETURN 580 CONTINUE WRITE (ITTY, 5010) STOP 5000 FORMAT (/5X,'Enter the number of components in SOLTHERM [37]') 5010 FORMAT (/5X,'Cannot find the H+ or the C123 entry in SOLTHERM'//) 5020 FORMAT (///5X, 'Data for complexes has been skipped in ', & 'SOLTHERM...', & /5X, 'Now starting to read data for gases and ', & 'minerals.', & /5X, '(Will read only up to ', I3, & ' minerals and gases)'/) 5030 FORMAT (/5X, 'Enter total number of gases in SOLTHERM [20] '/) 5040 FORMAT (/5X, 'Do you want to write a SOLTHERM formatted ', & 'output file', & /5X, 'with new Log K values ? [N] '//) 5050 FORMAT (///5X, 'Do you want to work on (1) Gases, (2) Minerals ', & /5X, 'or on (3) Complexes (placed after the gases in ', & 'SOLTHERM) ?'/) 5060 FORMAT (///5X, 'Do you want to select ', A8, ' by (1) Name, ', & '(2) Compostion, ' & /5X, 'or do you want (3) All of the ',A8,' ?'/) 5070 FORMAT (///5X, 'Input ', A8, ', one per line, blank to end,' & /5X, 'no more than 10 names'/) 5080 FORMAT (A8, 2X, F7.2, I2, 1X, 11(F6.3,I2), F6.3, I2) 5090 FORMAT (5X, 'Found ', A8, ' stochiometry in file SOLTHERM !') 5100 FORMAT (8X, 8F8.3, A8) 5110 FORMAT (68X, F6.3) 5120 FORMAT (10A8) 5130 FORMAT (/5X, 'Cannot find ',A8,' in SOLTHERM !') 5140 FORMAT (///5X, 'Input the indexes of desired species, ', & 'slash to end') 5150 FORMAT (5X, 3(I2, ': ', A8, 5X)) 5160 FORMAT (/5X, 'No ', A8, ' contains the selected species!'/) 5170 FORMAT (5X, A8, ' contains given species') 5180 FORMAT (5X, 'Could not find ',A8,' in the gas data file !') 5190 FORMAT (///5X, 'Okay.... the calculation is finished !'/) 5200 FORMAT (/5X, 'Do you want to go back to the selection process ?'/) 5210 FORMAT (5X, 'Found ',A8,' thermo data in Mineral file') 5220 FORMAT(//5X,'Could not find ',A8,' in Mineral file') 5230 FORMAT ('1', //, 1X, 'Mineral ', 6X, 'Gref', 10X, 'Href', & 7X, 'Sref', 6X, 'Vref', 9X, 'A', & 8X, 'B*10**3', 5X, 'C*10**-5'/, & 1X, '--------', 2X, '------------', & 2X, '------------', 2X, '------', & 2X, '----------', 2X, '----------', & 2X, '----------', 2X, '----------') 815 FORMAT (A8, A8, I1, A8, 6(F6.3, I2), F5.3, I2) 823 FORMAT (A8, 60X, F6.3) 826 FORMAT (A8, 8F8.3, A8) END