

      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 file
*   INP2  =  2 : Input file
*   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                              |
* ---------------------------------------------------------------------|

      Dum(1:8) = '12345678'
      Do while (DUM(1:8) .NE. 'H+      ')
         READ (INP1, '(A80)', END = 333) DUM
      End do

      READ (DUM, 500,ERR=777) SPECNM(1), CHG(1), WT(1)
  500 FORMAT (A8, 12X, F2.0, 17X, F8.4)

      write (*, 501) SPECNM(1), CHG(1), WT(1)
  501 FORMAT (1x,A8, 12X, F2.0, 17X, F8.4)

      Bstot = 1
      Do while (specnm(bstot) .NE. '        ')
        bstot = bstot + 1
        READ (INP1,500, ERR=777) SPECNM(bstot), CHG(bstot), WT(bstot)
      End do
      bstot = bstot - 1

*** 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 it is aqs. the species, read it
            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

  333 WRITE (*,'(//2x,A,A,/2x,A)') 'ERROR: EOF encountered before ',
     +         'finding H+ component line in soltherm.',' Exiting...'
      STOP

  777 WRITE (*,'(2X,A,A,A)') 'Format error in record = ',NAME,
     +      ' Exiting...'
      STOP

 3999 CLOSE (INP1)
      CLOSE (INP2)
      CLOSE (IOUT1)

      WRITE (*,'(//,A,/)') '  Execution finished.'
      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
*    <https://www.gnu.org/licenses/>.

*----------------------------------------------------------------------|

      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
*    <https://www.gnu.org/licenses/>.

*----------------------------------------------------------------------|

      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
