PROGRAM XPT_LOGK * =====================================================================| * Highly modified version of SOLPRINT, to extract and interpolate | * log K's at T's and P's specifed by the user. | * J.Palandri, July.2018 | * INP1 = 1 : SOLTHERM datafile | * INP2 = 2 : Input datafile | * IOUT1 = 10 : Output file | * | * Input file format is T and P pairs separated by at least one blank | * space or comma, followed by the names of derived species, | * minerals, or gases, with no leading blanks, spelled exactly as | * in SOLTHERM. | * | * Current limit is 1000 T, P pairs. | * | * ---------------------------------------------------------------------| IMPLICIT NONE LOGICAL MINSWITCH, PT_FOUND, TWARN INTEGER INP1, INP2, IOUT1, PINDEX, TINDEX, TIMEVAL(8), ITEMPC INTEGER BSTOT, I, J, K, ISOL, ISPEC(9), ITOT, JJ, KK, L CHARACTER * 1 PT_ERROR, ISO_ERROR CHARACTER * 3 MONTH(12) CHARACTER * 4 PTFLAG CHARACTER * 8 DAT(10), REF, SPECNM(100), DATE CHARACTER * 10 TIME CHARACTER * 20 INP_NAME, NAME, NAME2 CHARACTER * 30 FORMULA CHARACTER * 80 DUM CHARACTER * 120 TITLE1, TITLE2 REAL*8 B, CHARGE, CHG(100), COEF(9), DENS, VOL, MAXT REAL*8 WEIGHT, WT(100), ZZILCH(37,36), ZILCH(15) REAL*8 TEMPC(1000), PRESS(1000), YILCH(8), VAL, LGK COMMON /STUFF/ CHG, COEF, WT, ISPEC, SPECNM DATA MONTH /'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug', + 'Sep','Oct','Nov','Dec'/ DATA PRESS, TEMPC /1000*0.d0,1000*0.d0 / * =====================================================================| PARAMETER (INP1 = 1, INP2 = 2, IOUT1 = 10) OPEN (INP1, FILE = 'SOLTHERM.XPT', STATUS = 'OLD') OPEN (INP2, FILE = 'XPT-logk.dat', STATUS = 'OLD') OPEN (IOUT1, FILE = 'XPT-LogK.out', STATUS = 'UNKNOWN') MINSWITCH = .FALSE. WEIGHT = 0.0 BSTOT = 0 FORMULA = ' ' INP_NAME = ' ' NAME = ' ' NAME2 = ' ' READ (INP1, *, ERR=277) BSTOT 277 IF (BSTOT .LE. 0) THEN WRITE (*,'(//2X, A)') + 'Input total number of component species in soltherm >>' READ (*,*) BSTOT END IF WRITE (*,*) * ---------------------------------------------------------------------| * Down to statement 40 we skip all the records before the component| * species, then read these species in | * ---------------------------------------------------------------------| 10 CONTINUE READ (INP1, '(A80)', END = 20) DUM IF (DUM(1:8) .NE. 'H+ ') GO TO 10 GO TO 30 20 CONTINUE WRITE (*,'(//2x,A,A)') 'ERROR: EOF encountered before ', + 'reaching H+ entry line in soltherm ***********' STOP 30 CONTINUE READ (DUM, 500,ERR=777) SPECNM(1), CHG(1), WT(1) 500 FORMAT (A8, 12X, F2.0, 17X, F8.4) DO 40 I = 2, BSTOT READ (INP1, 500,ERR=777) SPECNM(I), CHG(I), WT(I) 40 CONTINUE *** Read user name, write with date. CALL DATE_AND_TIME(VALUES=TIMEVAL) READ (INP2, '(A120)') TITLE1 READ (INP2, '(A120)') TITLE2 WRITE (IOUT1, '(A120)') TITLE1 WRITE (IOUT1, '(A120)') TITLE2 IF (TIMEVAL(3).GE.10) THEN WRITE (IOUT1,'(A,I4,A,A,A,I2,A,A)') + 'Log K''s extracted from Soltherm on ', + TIMEVAL(1),'-',MONTH(TIMEVAL(2)),'-',TIMEVAL(3) ELSE WRITE (IOUT1,'(A,I4,A,A,A,I1,A,A)') + 'Log K''s extracted from Soltherm on ', + TIMEVAL(1),'-',MONTH(TIMEVAL(2)),'-',TIMEVAL(3) ENDIF * ---------------------------------------------------------------------| * Interpolate logK's | * ---------------------------------------------------------------------| *** Read input file for T's and P's JJ = 0 110 JJ = JJ + 1 READ (INP2, *, ERR = 140, END = 3999) TEMPC(JJ), PRESS(JJ) GOTO 110 140 JJ = JJ - 1 BACKSPACE INP2 READ (INP2,'(A80)') DUM WRITE (*,*) WRITE (*,'(A)') ' Reading T, P pairs.' WRITE (*,'(A)') ' Non-numeric character found in line:' WRITE (*,'(4x,A)') DUM WRITE (*,*) WRITE (*,'(A)') ' Reading mineral, gas, and species names,' WRITE (*,'(A)') ' and calculating Log K''s' BACKSPACE INP2 *** If Non-numeric character, then current record is a species, mineral, *** or gas. Rewind soltherm, reset variables. 50 REWIND INP1 MINSWITCH = .FALSE. WEIGHT = 0.0 CHARGE = 0.0 DO 85 I = 1, 9 COEF(I) = 0.0 ISPEC(I) = 0 85 CONTINUE READ (INP2,'(A20)',END=3999) INP_NAME *** Find in SOLTHERM 155 READ (INP1,'(A20)',END=2999) NAME IF (NAME(1:7).EQ.'B12 COE') MINSWITCH = .TRUE. 157 IF (NAME.EQ.INP_NAME) THEN BACKSPACE INP1 WRITE (IOUT1, *) IF (MINSWITCH.EQV..FALSE.) THEN *** Then is aqs. species, read header READ (INP1, '(A20,A30,T91,A4)',ERR=777) + NAME, FORMULA, PTFLAG READ (INP1, '(8X,F3.0,T74,A8)',ERR=777) CHARGE, REF READ (INP1, '(I2,1X,9(2x,F8.3,I2))',ERR=777) + ITOT, (COEF(I), ISPEC(I),I = 1, ITOT) *** If H+ is not in the I = 1 position, flag it DO 60 I = 1, ITOT IF(I.GT.1.AND.ISPEC(I).EQ.1) THEN WRITE(IOUT1,663) NAME WRITE(*,663) NAME 663 FORMAT (' WARNING: Species ',A20, ' has H+ in a ', + 'position beyond I = 2.') ENDIF 60 CONTINUE ELSE *** Else is a gas or mineral, read header READ (INP1, '(A20,A30,T67,I2,T91,A4)',ERR=777) + NAME, FORMULA, ISOL, PTFLAG READ (INP1,'(T13,F8.3,T37,F8.3,T58,F7.3,T74,A8)',ERR=777) + WEIGHT, VOL, DENS, REF READ (INP1, '(I2,1X,9(2x,F8.3,I2))',ERR=777) + ITOT, (COEF(I), ISPEC(I),I = 1, ITOT) IF (ISOL .EQ. 6) READ (INP1, *) ENDIF IF (PTFLAG .EQ. 'All ') THEN *** Then read the full array DO 992, J = 1,37 READ (INP1,*,ERR=777) (ZZILCH(J,K), K = 1,36) 992 CONTINUE ELSE *** Else read the one line of L-V data READ (INP1,*,ERR=777) (YILCH(L), L = 1, 8) ENDIF ELSE GOTO 155 ENDIF *** Write header, stochiometry write (iout1,*) CALL PRINT (FORMULA, NAME, ITOT, CHARGE, WEIGHT, MINSWITCH) IF (MINSWITCH.EQV..TRUE.) THEN WRITE (IOUT1,'(2X, A10,F7.2,2X,A20,F8.2,2X,A13, F5.2)') + 'mol.wt. = ',WEIGHT,' mol.vol.(cc/mol) = ', + VOL,' rho(g/cc) = ',DENS ENDIF WRITE (IOUT1,'(2A8)') " Ref = ", REF WRITE (IOUT1,*) WRITE (IOUT1,'(A)') ' T(°C) P(bar) Log(K)' *** Interpolate logK's, write to output DO 1999, KK = 1, JJ IF (PTFLAG .EQ. 'All ') THEN *** Then full grid PT_ERROR = 'N' ISO_ERROR = 'N' CALL PT_INDEXER (PRESS(KK), TEMPC(KK), PINDEX, TINDEX, + PT_FOUND, PT_ERROR, ISO_ERROR) IF (PT_ERROR.EQ.'E'.OR.ISO_ERROR.EQ.'E') THEN WRITE (*,'(2X,A6,A20)') 'Entry ',NAME WRITE (*,994) TEMPC(KK), PRESS(KK) WRITE (IOUT1,994) TEMPC(KK), PRESS(KK) 994 FORMAT (1X,F7.2,2X,F8.3,1X + '99999.999 ERROR: Specified P and/or T are', + ' outside the P-T range of the data.') GOTO 1999 ENDIF CALL INTRP3D(PRESS(KK),TEMPC(KK),PINDEX,TINDEX,ZZILCH,VAL) LGK=VAL IF (PT_ERROR.EQ.'W'.OR.ISO_ERROR.EQ.'W') THEN WRITE (*,'(2X,A6,A20)') 'Entry ',NAME WRITE (IOUT1,999) TEMPC(KK), PRESS(KK), LGK WRITE (*,999) TEMPC(KK), PRESS(KK), LGK 999 FORMAT (1X,F7.2,2X,F8.3,1X,F9.3,2X'WARNING: P and/or T', + ' are in a grid block at the limit of the P-T', + ' range of data.') ELSE WRITE (IOUT1, '(1X,F7.2,2X,F8.3,1X,F9.3)') + TEMPC(KK), PRESS(KK), LGK ENDIF ELSE *** Else L-V saturation only CALL INTRP2D (TEMPC(KK), TWARN, MAXT, YILCH, VAL) LGK=VAL IF (YILCH(1).GT.99999.)THEN WRITE (*,'(2X,A6,A20)') 'Entry ',NAME WRITE (*,710) TEMPC(KK), PRESS(KK), LGK WRITE (IOUT1,710) TEMPC(KK), PRESS(KK),LGK 710 FORMAT (1X,F7.2,2X,F8.3,1x,F9.3 + ' ERROR: format error or data are ', + ' absent in soltherm') ELSEIF (MAXT.LT.(TEMPC(KK)-50.D0)) THEN ITEMPC = INT(MAXT) WRITE (*,'(2X,A6,A20)') 'Entry ',NAME WRITE (*,720) TEMPC(KK), PRESS(KK), LGK, ITEMPC WRITE (IOUT1,720) TEMPC(KK), PRESS(KK), LGK, ITEMPC 720 FORMAT (1X,F7.2,2X,F8.3,1X,F9.3, + ' ERROR: Data extend to ',I3,' °C only. Max', + ' extrapolation is 50 °C') ELSEIF (TEMPC(KK).LT.25.D0.OR.TEMPC(KK).GT.350.D0) THEN WRITE (*,'(2X,A6,A20)') 'Entry ',NAME WRITE (*,730) TEMPC(KK), PRESS(KK), LGK WRITE (IOUT1,730) TEMPC(KK), PRESS(KK), LGK 730 FORMAT (1X,F7.2,2X,F8.3,1X,F9.3, + ' WARNING: LogK extrapolated. Data limited to', + ' L-V saturation where 25 °C <= T <= 350 °C.') ELSE WRITE (IOUT1, '(1X,F7.2,2X,F8.3,1X,F9.3,2X,A,A)') + TEMPC(KK), PRESS(KK), LGK, 'Data for L-V sat', + ' only, pressure ignored.' ENDIF ENDIF 1999 CONTINUE *** Check for second soltherm entry in case of redox READ (INP1,'(/,A20)') NAME IF (NAME.EQ.INP_NAME) GOTO 157 GO TO 50 *** If input species or mineral not found in soltherm 2999 WRITE (IOUT1, *) WRITE (IOUT1, '(2X,A,A,A)') 'Entry ''',INP_NAME, + ''' not found in SOLTHERM.' WRITE (*, '(2X,A,A,A)') 'Entry ''',INP_NAME, + ''' not found in SOLTHERM.' GOTO 50 777 WRITE (*,'(2X,A25,A20)') 'Format error in record = ',NAME CLOSE (INP1) CLOSE (INP2) CLOSE (IOUT1) STOP 3999 CLOSE (INP1) CLOSE (INP2) CLOSE (IOUT1) WRITE (*,'(//,A,/)') ' DONE!' END SUBROUTINE PRINT (FORMULA, NAME, ITOT, CHARGE, WEIGHT, MINSWITCH) * =====================================================================| IMPLICIT NONE INTEGER INP1, IOUT1 INTEGER I, IPROD, IREAC, ISPEC(9), ITOT, PP, TRIGHT, TRIM CHARACTER * 3 CHAR1(9), CHAR2(9) CHARACTER * 8 PRODNM(9), REACNM(9), SPECNM(100) CHARACTER * 20 NAME CHARACTER * 30 FORMULA, REACNL, BLANK30 CHARACTER * 92 BLANK, BUFFER LOGICAL MINSWITCH REAL*8 CHARGE, CHBAL, CHG(100), COEF(9), PRODCF(9) REAL*8 REACF(9), WEIGHT, WT(100), WTM EXTERNAL TRIGHT COMMON /STUFF/ CHG, COEF, WT, ISPEC, SPECNM DATA CHAR1 /' ', 8*' + '/ DATA CHAR2 /' = ', 8*' + '/ DATA BLANK /' '/ DATA BLANK30 /' '/ * =====================================================================| INP1 = 1 IOUT1 = 10 IREAC = 1 IPROD = 0 CHBAL = 0.0 WTM = 0.0 REACF(1) = 1.0 REACNL = BLANK30 IF (FORMULA.NE.BLANK30) THEN IF (MINSWITCH.EQV..TRUE.) THEN WRITE (IOUT1,'(1X,A20)') NAME ENDIF REACNL = FORMULA ELSE REACNL(1:20) = NAME ENDIF DO 40 I = 1, ITOT * Add charges and weights of components species to check for balance, below. CHBAL = CHBAL + COEF(I) * CHG(ISPEC(I)) WTM = WTM + COEF(I) * WT(ISPEC(I)) * Parse positive and negative charge for components with branched IF IF (COEF(I)) 10, 10, 20 10 CONTINUE IREAC = IREAC + 1 REACF(IREAC) = -1.0 * COEF(I) REACNM(IREAC) = SPECNM(ISPEC(I)) GO TO 30 20 IPROD = IPROD + 1 PRODCF(IPROD) = COEF(I) PRODNM(IPROD) = SPECNM(ISPEC(I)) 30 CONTINUE 40 CONTINUE * ---------------------------------------------------------------------| * prints reaction in readable form | * ---------------------------------------------------------------------| PP = 1 BUFFER = ' ' DO 43 I = 1, IREAC IF (I .NE. 1) THEN WRITE (BUFFER(PP:), '(A3)') CHAR1(I) PP = PP + 3 END IF IF (TRIGHT(REACF(I)) .EQ. 1) THEN WRITE (BUFFER(PP:), '(F4.1)') REACF(I) PP = PP + 4 ELSE IF (TRIGHT(REACF(I)) .EQ. 2) THEN WRITE (BUFFER(PP:), '(F5.2)') REACF(I) PP = PP + 5 ELSE WRITE (BUFFER(PP:), '(F6.3)') REACF(I) PP = PP + 6 END IF IF (I.EQ.1) THEN WRITE (BUFFER(PP:), '(2X, A30)') REACNL TRIM = LEN_TRIM(REACNL) PP = PP + TRIM + 3 ELSE WRITE (BUFFER(PP:), '(1X, A8)') REACNM(I) TRIM = LEN_TRIM(REACNM(I)) PP = PP + TRIM + 2 ENDIF IF (PP + 20 .GE. 92) THEN WRITE (IOUT1, '(A92)') BUFFER PP = 2 BUFFER = ' ' END IF 43 CONTINUE IF (BUFFER .NE. BLANK) THEN WRITE (IOUT1, '(A92)') BUFFER END IF PP = 5 BUFFER = ' ' DO 45 I = 1, IPROD WRITE (BUFFER(PP:), '(A3)') CHAR2(I) PP = PP + 3 IF (TRIGHT(PRODCF(I)) .EQ. 1) THEN WRITE (BUFFER(PP:), '(F4.1)') PRODCF(I) PP = PP + 4 ELSE IF (TRIGHT(PRODCF(I)) .EQ. 2) THEN WRITE (BUFFER(PP:), '(F5.2)') PRODCF(I) PP = PP + 5 ELSE WRITE (BUFFER(PP:), '(F6.3)') PRODCF(I) PP = PP + 6 END IF WRITE (BUFFER(PP:), '(1X, A8)') PRODNM(I) TRIM = LEN_TRIM(PRODNM(I)) PP = PP + TRIM + 2 IF (PP + 20 .GE. 92) THEN WRITE (IOUT1, '(A92)') BUFFER PP = 5 BUFFER = ' ' END IF 45 CONTINUE IF (BUFFER .NE. BLANK) THEN WRITE (IOUT1, '(A92)') BUFFER END IF CHBAL = CHBAL - CHARGE WTM = WTM - WEIGHT IF (ABS(CHBAL) .LE. 0.000010) GO TO 50 WRITE (IOUT1, 412) CHBAL 412 FORMAT (/5X,' *** WARNING: charge balance is: ',F7.3/) WRITE (*,413) CHBAL, NAME 413 FORMAT(/5X,' *** WARNING: charge balance is: ',F7.3,' FOR ',A20) 50 CONTINUE IF(WEIGHT.EQ.0.) RETURN IF(ABS(WTM).LT.0.1) RETURN WRITE (IOUT1, 415) WTM 415 FORMAT(/5X,' *** WARNING: mol.weight balance is: ',F8.3/) WRITE (*,416) WTM, NAME 416 FORMAT(/5X,'*** WARNING: mol.weight balance is: ',F8.3, + ' FOR ',A20) RETURN END INTEGER FUNCTION TRIGHT (STUFF) * =====================================================================| * Returns 1 (for 0 or 1), 2 or 3, the number of significant figures for * reaction coefficients, to minimize the number of blanks. REAL * 8 CPUMIN, FRACTN, NUMBER, STUFF PARAMETER (CPUMIN = 1.0D-35) INTRINSIC DINT, DABS * =====================================================================| NUMBER = STUFF FRACTN = DINT(NUMBER) - NUMBER IF (DABS(FRACTN) .LE. CPUMIN) THEN TRIGHT = 1 RETURN END IF NUMBER = NUMBER * 10 FRACTN = DINT(NUMBER) - NUMBER IF (DABS(FRACTN) .LE. CPUMIN) THEN TRIGHT = 1 RETURN END IF NUMBER = NUMBER * 10 FRACTN = DINT(NUMBER) - NUMBER IF (DABS(FRACTN) .LE. CPUMIN) THEN TRIGHT = 2 RETURN END IF TRIGHT = 3 RETURN END SUBROUTINE PT_INDEXER (P,T,PIDX,TIDX,PT_FOUND,PT_ERROR,ISO_ERROR) *======================================================================| * * Finds array indices PIDX and TIDX using specified P, T * *----------------------------------------------------------------------| IMPLICIT NONE LOGICAL PT_FOUND INTEGER PIDX, TIDX CHARACTER PT_ERROR, ISO_ERROR REAL*8 P, T *======================================================================| **** Assign array index integer for pressure (PIDX). **** PT_ERROR = 'N' if no warning **** PT_ERROR = 'W' if P > 0.0 and P < 1.0 or **** if P > 5000.0 and P < 6000.0 or **** if T > 600.0 and T <= 700.0 **** PT_ERROR = 'E' if P <= 0.0 or **** if P >= 6000.0 or **** if T > 700.0 or **** if T < 0.01 PT_ERROR = 'N' IF (P.LE.0.0D0) THEN PIDX = 1 PT_ERROR = 'E' ELSE IF (P.GT.0.0D0.AND.P.LT.1.0D0) THEN PIDX = 1 PT_ERROR = 'W' ELSE IF (P.GE.1.0D0.AND.P.LT.1.0133D0) THEN PIDX = 1 ELSE IF (P.GE.1.0133D0.AND.P.LT.2.321D0) THEN PIDX = 2 ELSE IF (P.GE.2.321D0.AND.P.LT.4.758D0) THEN PIDX = 3 ELSE IF (P.GE.4.758D0.AND.P.LT.8.919D0) THEN PIDX = 4 ELSE IF (P.GE.8.919D0.AND.P.LT.15.537D0) THEN PIDX = 5 ELSE IF (P.GE.15.537D0.AND.P.LT.25.479D0) THEN PIDX = 6 ELSE IF (P.GE.25.479D0.AND.P.LT.39.737D0) THEN PIDX = 7 ELSE IF (P.GE.39.737D0.AND.P.LT.59.432D0) THEN PIDX = 8 ELSE IF (P.GE.59.432D0.AND.P.LT.85.839D0) THEN PIDX = 9 ELSE IF (P.GE.85.839D0.AND.P.LT.120.458D0) THEN PIDX = 10 ELSE IF (P.GE.120.458D0.AND.P.LT.165.212D0) THEN PIDX = 11 ELSE IF (P.GE.165.212D0.AND.P.LT.200.0D0) THEN PIDX = 12 ELSE IF (P.GE.200.0D0.AND.P.LT.250.0D0) THEN PIDX = 13 ELSE IF (P.GE.250.0D0.AND.P.LT.300.0D0) THEN PIDX = 14 ELSE IF (P.GE.300.0D0.AND.P.LT.350.0D0) THEN PIDX = 15 ELSE IF (P.GE.350.0D0.AND.P.LT.400.0D0) THEN PIDX = 16 ELSE IF (P.GE.400.0D0.AND.P.LT.450.0D0) THEN PIDX = 17 ELSE IF (P.GE.450.0D0.AND.P.LT.500.0D0) THEN PIDX = 18 ELSE IF (P.GE.500.0D0.AND.P.LT.550.0D0) THEN PIDX = 19 ELSE IF (P.GE.550.0D0.AND.P.LT.600.0D0) THEN PIDX = 20 ELSE IF (P.GE.600.0D0.AND.P.LT.650.0D0) THEN PIDX = 21 ELSE IF (P.GE.650.0D0.AND.P.LT.700.0D0) THEN PIDX = 22 ELSE IF (P.GE.700.0D0.AND.P.LT.750.0D0) THEN PIDX = 23 ELSE IF (P.GE.750.0D0.AND.P.LT.800.0D0) THEN PIDX = 24 ELSE IF (P.GE.800.0D0.AND.P.LT.850.0D0) THEN PIDX = 25 ELSE IF (P.GE.850.0D0.AND.P.LT.900.0D0) THEN PIDX = 26 ELSE IF (P.GE.900.0D0.AND.P.LT.950.0D0) THEN PIDX = 27 ELSE IF (P.GE.950.0D0.AND.P.LT.1000.0D0) THEN PIDX = 28 ELSE IF (P.GE.1000.0D0.AND.P.LT.1500.0D0) THEN PIDX = 29 ELSE IF (P.GE.1500.0D0.AND.P.LT.2000.0D0) THEN PIDX = 30 ELSE IF (P.GE.2000.0D0.AND.P.LT.2500.0D0) THEN PIDX = 31 ELSE IF (P.GE.2500.0D0.AND.P.LT.3000.0D0) THEN PIDX = 32 ELSE IF (P.GE.3000.0D0.AND.P.LT.3500.0D0) THEN PIDX = 33 ELSE IF (P.GE.3500.0D0.AND.P.LT.4000.0D0) THEN PIDX = 34 ELSE IF (P.GE.4000.0D0.AND.P.LE.5000.0D0) THEN PIDX = 35 ELSE IF (P.GT.5000.0D0.AND.P.LT.6000.0D0) THEN PIDX = 35 PT_ERROR = 'W' ELSE IF (P.GE.6000.0D0) THEN PIDX = 35 PT_ERROR = 'E' ENDIF **** Assign array index integer for temperature (TIDX) **** ISO_ERROR = 'N' if no warning **** ISO_ERROR = 'W' if P/T near L-V saturation or **** if P/T near critical isochore, rho H2O < 0.35g/cc **** ISO_ERROR = 'E' if P/T fully in the two-phase region or **** if P/T fully above critical isochore. IF (T.LT.0.01D0) THEN PT_ERROR = 'E' TIDX = 3 IF (PIDX.GT.34) PIDX = 34 ELSE IF (T.GE.0.01D0.AND.T.LE.50.0D0) THEN TIDX = 3 IF (PIDX.GT.34) PIDX = 34 ELSE IF (T.GT.50.0D0.AND.T.LE.75.0D0) THEN TIDX = 4 ELSE IF (T.GT.75.0D0.AND.T.LE.100.0D0) THEN TIDX = 5 IF (P.LT.1.0133D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF ELSE IF (T.GT.100.0D0.AND.T.LE.125.0D0) THEN TIDX = 6 IF (P.LT.2.321D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.1.013D0) ISO_ERROR = 'E' ELSE IF (T.GT.125.0D0.AND.T.LE.150.0D0) THEN TIDX = 7 IF (P.LT.4.758D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.2.321D0) ISO_ERROR = 'E' ELSE IF (T.GT.150.0D0.AND.T.LE.175.0D0) THEN TIDX = 8 IF (P.LT.8.919D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.4.758D0) ISO_ERROR = 'E' ELSE IF (T.GT.175.0D0.AND.T.LE.200.0D0) THEN TIDX = 9 IF (P.LT.15.537D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.8.919D0) ISO_ERROR = 'E' ELSE IF (T.GT.200.0D0.AND.T.LE.225.0D0) THEN TIDX = 10 IF (P.LT.25.479D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.15.537D0) ISO_ERROR = 'E' ELSE IF (T.GT.225.0D0.AND.T.LE.250.0D0) THEN TIDX = 11 IF (P.LT.39.737D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.25.479D0) ISO_ERROR = 'E' ELSE IF (T.GT.250.0D0.AND.T.LE.275.0D0) THEN TIDX = 12 IF (P.LT.59.432D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.39.737D0) ISO_ERROR = 'E' ELSE IF (T.GT.275.0D0.AND.T.LE.300.0D0) THEN TIDX = 13 IF (P.LT.85.839D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.59.432D0) ISO_ERROR = 'E' ELSE IF (T.GT.300.0D0.AND.T.LE.325.0D0) THEN TIDX = 14 IF (P.LT.120.458D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.85.839D0) ISO_ERROR = 'E' ELSE IF (T.GT.325.0D0.AND.T.LE.350.0D0) THEN TIDX = 15 IF (P.LT.165.212D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.120.458D0) ISO_ERROR = 'E' ELSE IF (T.GT.350.0D0.AND.T.LE.400.0D0) THEN TIDX = 16 IF (P.GE.165.212D0.AND.P.LT.200.0D0) THEN PIDX = PIDX + 3 ISO_ERROR = 'W' ENDIF IF (P.GE.200.0D0.AND.P.LT.250.0D0) THEN PIDX = PIDX + 2 ISO_ERROR = 'W' ENDIF IF (P.GE.250.0D0.AND.P.LT.300.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.165.212D0) ISO_ERROR = 'E' ELSE IF (T.GT.400.0D0.AND.T.LE.410.0D0) THEN TIDX = 17 IF (P.GE.200.0D0.AND.P.LT.350.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.200.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.410.0D0.AND.T.LE.420.0D0) THEN TIDX = 18 IF (P.GE.350.0D0.AND.P.LT.400.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.350.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.420.0D0.AND.T.LE.430.0D0) THEN TIDX = 19 IF (P.GE.350.0D0.AND.P.LT.400.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.350.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.430.0D0.AND.T.LE.440.0D0) THEN TIDX = 20 IF (P.GE.400.0D0.AND.P.LT.450.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.400.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.440.0D0.AND.T.LE.450.0D0) THEN TIDX = 21 IF (P.GE.450.0D0.AND.P.LT.500.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.450.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.450.0D0.AND.T.LE.460.0D0) THEN TIDX = 22 IF (P.GE.450.0D0.AND.P.LT.500.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.450.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.460.0D0.AND.T.LE.470.0D0) THEN TIDX = 23 IF (P.GE.500.0D0.AND.P.LT.550.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.500.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.470.0D0.AND.T.LE.480.0D0) THEN TIDX = 24 IF (P.GE.550.0D0.AND.P.LT.600.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.550.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.480.0D0.AND.T.LE.490.0D0) THEN TIDX = 25 IF (P.GE.550.0D0.AND.P.LT.600.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.550.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.490.0D0.AND.T.LE.500.0D0) THEN TIDX = 26 IF (P.GE.600.0D0.AND.P.LT.650.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.600.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.500.0D0.AND.T.LE.510.0D0) THEN TIDX = 27 IF (P.GE.600.0D0.AND.P.LT.650.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.600.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.510.0D0.AND.T.LE.520.0D0) THEN TIDX = 28 IF (P.GE.650.0D0.AND.P.LT.700.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.650.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.520.0D0.AND.T.LE.530.0D0) THEN TIDX = 29 IF (P.GE.700.0D0.AND.P.LT.750.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.700.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.530.0D0.AND.T.LE.540.0D0) THEN TIDX = 30 IF (P.GE.700.0D0.AND.P.LT.750.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.700.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.540.0D0.AND.T.LE.550.0D0) THEN TIDX = 31 IF (P.GE.750.0D0.AND.P.LT.800.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.700.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.550.0D0.AND.T.LE.560.0D0) THEN TIDX = 32 IF (P.GE.800.0D0.AND.P.LT.850.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.800.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.560.0D0.AND.T.LE.570.0D0) THEN TIDX = 33 IF (P.GE.800.0D0.AND.P.LT.850.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.800.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.570.0D0.AND.T.LE.580.0D0) THEN TIDX = 34 IF (P.GE.850.0D0.AND.P.LT.900.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.850.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.580.0D0.AND.T.LE.590.0D0) THEN TIDX = 35 IF (P.GE.900.0D0.AND.P.LT.950.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.900.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.590.0D0.AND.T.LE.600.0D0) THEN TIDX = 36 IF (P.GE.900.0D0.AND.P.LT.950.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.900.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.600.0D0.AND.T.LE.700.0D0) THEN TIDX = 36 IF (PT_ERROR.EQ.'N') PT_ERROR = 'W' IF (P.GE.900.0D0.AND.P.LT.950.0D0) THEN PIDX = PIDX + 1 ISO_ERROR = 'W' ENDIF IF (P.LT.800.0D0) ISO_ERROR = 'E' ELSE IF (T.GT.700.0D0) THEN TIDX = 36 PT_ERROR = 'E' ENDIF PT_FOUND = .TRUE. RETURN END SUBROUTINE INTRP3D (P, T, PIDX, TIDX, LOGKARRAY, LOGK) *======================================================================| * * Interpolates intermediate LogK(P,T) values from LogK(P,T) array * using specified P, T, and array indices PIDX and TIDX found using * subroutine PT_INDEXER * *----------------------------------------------------------------------| * * (Interpolates intermdiate f(x,y) values from f(x,y) array using * specified x, y, and array indices XIDX and YTIDX found using * subroutine XY_INDEXER) * *----------------------------------------------------------------------| IMPLICIT NONE INTEGER PDIM, TDIM, PIDX, TIDX, I, J, K, L PARAMETER (PDIM = 37) PARAMETER (TDIM = 36) REAL*8 P, T, PVAL(PDIM), TVAL(TDIM), LOGKARRAY (PDIM,TDIM) REAL*8 SUM, LK(3), LOGK DATA PVal / 1.0d0, 1.0133d0, 2.321d0, 4.758d0, 8.919d0, 15.537d0, + 25.479d0, 39.737d0, 59.432d0, 85.839d0, 120.458d0, + 165.212d0, 200.d0, 250.d0, 300.d0, 350.d0, 400.d0, + 450.d0, 500.d0, 550.d0, 600.d0, 650.d0, 700.d0, + 750.d0, 800.d0, 850.d0, 900.d0, 950.d0, + 1000.d0, 1500.d0, 2000.d0, 2500.d0, 3000.d0, 3500.d0, + 4000.d0, 4500.d0, 5000.d0 / DATA TVal / 0.01d0, 25.d0, 50.d0, 75.d0, 100.d0, 125.d0, 150.d0, + 175.d0, 200.d0, 225.d0, 250.d0, 275.d0, 300.d0, + 325.d0, 350.d0, 400.d0, 410.d0, 420.d0, 430.d0, + 440.d0, 450.d0, 460.d0, 470.d0, 480.d0, 490.d0, + 500.d0, 510.d0, 520.d0, 530.d0, 540.d0, 550.d0, + 560.d0, 570.d0, 580.d0, 590.d0, 600.d0 / *======================================================================| LK(1) = 0.0D0 LK(2) = 0.0D0 LK(3) = 0.0D0 LOGK = 0.0D0 * Interpolate three logk's at grid at the selected pressure, at * temperatures bracketing the selected temperature * Note that lowest TIDX value passed to this subroutine is 3 * and the lowest PIDX value passed is 1. L = 1 DO 500, I = TIDX - 2, TIDX DO 400, J = PIDX, PIDX + 2 SUM = LOGKARRAY(J,I) DO 300, K = PIDX, PIDX + 2 IF (J.NE.K) SUM = SUM*(P-PVAL(K))/(PVAL(J)-PVAL(K)) 300 CONTINUE LK(L) = LK(L) + SUM 400 CONTINUE L = L + 1 500 CONTINUE * Interpolate log k at the selected temperature K = 1 DO 700, I = TIDX - 2, TIDX SUM = LK(K) DO 600, J = TIDX - 2, TIDX IF (I.NE.J) SUM = SUM*(T-TVAL(J))/(TVAL(I)-TVAL(J)) 600 CONTINUE LOGK = LOGK + SUM K = K + 1 700 CONTINUE RETURN END SUBROUTINE INTRP2D (TEMPC, TWARN, MAXT, LGKDAT, LOGK) *======================================================================| * * Interpolates log K(s) at L-V saturation for non-Supcrt L-V data at * 25, 50, 100, 150, 200, 250, 300, 350 * * ---------------------------------------------------------------------| IMPLICIT NONE LOGICAL TWARN INTEGER I, IT, J, MAXTIX, IOUT1 REAL*8 TEMPC, LGKDAT(8), TDATA(8), LOGK, MAXT, SUM DATA TDATA + /25.D0, 50.D0, 100.D0, 150.D0, 200.D0, 250.D0, 300.D0, 350.D0/ *======================================================================| TWARN = .FALSE. **Max TEMPC value and array index for that value (TDATA(8) = 350 currently) IOUT1 = 11 MAXTIX = 8 MAXT = 350.D0 LOGK = 0.D0 IT = 0 **Find highest TEMPC for which data exist, and array index for that TEMPC ** log K = 99999.999 where data do not exist. DO 100, I = 8, 2, -1 IF (LGKDAT(I).EQ.99999.999D0) THEN MAXT = TDATA(I - 1) MAXTIX = I - 1 ENDIF 100 CONTINUE **Flag if Temperaute is outside the range for which log Ks exist in soltherm. IF (TEMPC.GT.TDATA(MAXTIX).OR.TEMPC.LT.TDATA(1)) TWARN = .TRUE. **If logk only at 25C, set logK = LGKDAT(1) IF (MAXTIX.EQ.1) THEN IF (TEMPC.LE.35.D0.AND.TEMPC.GE.10.D0) THEN LOGK = LGKDAT(1) ELSE LOGK = 99999.999D0 ENDIF **Else if logk only at 25 and 50C, inter/extrapolate w/ LGKDAT(1) and LGKDAT(2) ELSE IF (MAXTIX.EQ.2) THEN IF (TEMPC.GE.0.01D0.AND.TEMPC.LE.75.D0) THEN LOGK = LGKDAT(1) + (TEMPC - 25.D0)* + (LGKDAT(2)-LGKDAT(1))/(50.d0-25.D0) ELSE LOGK = 99999.999D0 ENDIF ELSE **Compute desired log K value for the given TEMPC-value (T)....... ** The following should work wherever a mineral has data to 100C or above. ** Find where the desired TEMPC value lies in the TEMPC data array. DO 150, I = 1,7 IF (TEMPC.GT.TDATA(I).AND.TEMPC.LE.TDATA(I+1)) IT = I 150 CONTINUE IF (TEMPC.GT.TDATA(8)) IT = 7 **If TEMPC .gt. the first TEMPC-array value that is 99999.999 or ** .gt. 400C, set logK to 99999.999. Max extrapolation is 50C. * IF (IT.GE.(MAXTIX+1).OR.TEMPC.GT.400.D0) THEN IF (TEMPC.GT.TDATA(MAXTIX)+50.D0.OR.TEMPC.GT.400.D0) THEN LOGK = 99999.999D0 ELSE IF (IT.LT.2) THEN *** If TEMPC .LT. the second TDATA value, use the 3 smallest TDATA values DO 300, I = 1, 3 SUM = LGKDAT(I) DO 200, J = 1,3 IF (I.NE.J) + SUM = SUM*(TEMPC-TDATA(J))/(TDATA(I)-TDATA(J)) 200 CONTINUE LOGK = LOGK + SUM 300 CONTINUE ELSE IF (IT.GT.MAXTIX-1) THEN **** If TEMPC .GT. next to last TDATA value, use 3 largest TDATA values DO 500, I = MAXTIX-2, MAXTIX SUM = LGKDAT(I) DO 400, J = MAXTIX-2, MAXTIX IF (I.NE.J) + SUM = SUM*(TEMPC-TDATA(J))/(TDATA(I)-TDATA(J)) 400 CONTINUE LOGK = LOGK + SUM 500 CONTINUE ELSE **** Else use TDATA values two below, one above TEMPC DO 700, I = IT-1, IT+1 SUM = LGKDAT(I) DO 600, J = IT-1, IT+1 IF (I.NE.J) + SUM = SUM*(TEMPC -TDATA(J))/(TDATA(I)-TDATA(J)) 600 CONTINUE LOGK = LOGK + SUM 700 CONTINUE ENDIF ENDIF RETURN END