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.
* Copyright 2018, 2024 James Palandri
* 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
* Copyright 2023 James Palandri
* James Palandri hereby disclaims all copyright interest in the
* Subroutine "PT_INDEXER" (which defines a range of temperatures
* and pressures in which subsequent calculation are valid).
* written by James Palandri
* Signature of James Palandri, 31 January 2024
* James Palandri
* This file is part of subroutine PT_INDEXER.
* PT_INDEXER is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published
* by the Free Software Foundation, either version 3 of the License,
* or (at your option) any later version. PT_INDEXER is distributed in
* the hope that it will be useful, but WITHOUT ANY WARRANTY; without
* even the implied warranty of MERCHANTABILITY or FITNESS FOR A
* PARTICULAR PURPOSE. See the GNU General Public License for more
* details. You should have received a copy of the GNU
* General Public License along with SuprxnGT. If not, see
* .
*----------------------------------------------------------------------|
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)
* Copyright 2023 James Palandri
* James Palandri hereby disclaims all copyright interest in the
* Subroutine "INTRP3D" (which interpolates intermdiate logK values
* from known logK values that bracket the interpolated value).
* written by James Palandri
* Signature of James Palandri, 31 January 2024
* James Palandri
* This file is part of subroutine INTRP3D.
* INTRP3D is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published
* by the Free Software Foundation, either version 3 of the License,
* or (at your option) any later version. INTRP3D is distributed in
* the hope that it will be useful, but WITHOUT ANY WARRANTY; without
* even the implied warranty of MERCHANTABILITY or FITNESS FOR A
* PARTICULAR PURPOSE. See the GNU General Public License for more
* details. You should have received a copy of the GNU
* General Public License along with SuprxnGT. If not, see
* .
*----------------------------------------------------------------------|
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