*=======================================================================
*
*     Program Alltst
*
*=======================================================================

      LOGICAL    PTfind

      CHARACTER  PTerr, isoerr, repeat

      INTEGER    Pidx, Tidx, i, j
      INTEGER    INP1, INP2, IOUT1

      REAL*8     P, T, LKarra(20,21), LogK

      INP1  = 11
      INP2  = 12
      IOUT1 = 21
      OPEN (INP1,  FILE = 'quartz.dat',  STATUS = 'OLD')
      OPEN (INP2,  FILE = 'PT.dat',      STATUS = 'OLD')
      OPEN (IOUT1, FILE = 'logk.out',    STATUS = 'Unknown')

      Do 10, i = 1,20
         Read (inp1,*) (LKarra (i,j), j = 1,21)
10    Continue

      PTfind = .false.

100   Write (*,*)

      PTerr  = 'N'
      isoerr = 'N'
      Pidx = -100
      Tidx = -100
      P = -9999.9d0
      T = -9999.9d0

      read  (INP2,*) P,T

      Call FindPT (P, T, Pidx, Tidx, PTfind, PTerr, isoerr)
      Call FindLK (P, T, Pidx, Tidx, LKarra, LogK)

      write (IOUT1,'(2x,A8,2x,A11)') 'Pressure','Temperature'
      write (IOUT1,'(2x,F8.3,2x,F8.3)') P,T
      Write (IOUT1,*)
      Write (IOUT1,'(2x,A6,4x,A6)') 'Pidx','Tidx'
      write (IOUT1,'(2x,I4,6x,I4)') Pidx, Tidx
      Write (IOUT1,*)
      Write (IOUT1,'(A8,F8.3)') 'Log K = ', logK

      If (isoerr.eq.'E') then
         Write (IOUT1,'(A25,A38)') ' Specified P-T conditions',
     +         ' are too far into the two-phase region'
         Write (IOUT1,'(A24,A41)') '  or too far outside the',
     +         ' critical isochore, rho(H2O) << 0.35g/cc.'
         Write (IOUT1,'(A23)') '  Execution terminated.'
         stop
      Endif

      If (PTerr.eq.'E'.or.PTerr.eq.'E') then
         Write (IOUT1,'(A32,A52)') '  Specified P and T fall too far',
     +         ' outside the P-T data array (.01-600 C, 1-5000 bar).'
         Write (IOUT1,'(A23)') '  Execution terminated.'
         stop
      Endif

      If (isoerr.eq.'W'.or.isoerr.eq.'w') then
         Write (*,*)
         Write (IOUT1,'(A47)')
     +         ' Warning: Specified P-T conditions are near the'
         Write (IOUT1,'(A41)')
     +         '  liquid-vapor saturation boundary or the'
         Write (IOUT1,'(A41)')
     +         '  critical isochore, rho(H2O) = 0.35g/cc.'
      Endif

      If (PTerr.eq.'W'.or.PTerr.eq.'w') then
         Write (IOUT1,*)
         Write (IOUT1,'(A38)')
     +         ' Warning: pressure and/or temperature:'
         Write (IOUT1,'(4x,F8.2,A8,F8.2,A5)') P,' bar,   ',T,' degC'
         Write (IOUT1,'(A31,A25)') ' is outside range of P-T array, ',
     +         '(1-5000 bar, .01-600 C).'
         Write (IOUT1,'(A28)') ' Log K will be extrapolated.'
      Endif

      Write (IOUT1,'(//)')

      goto 100
      end

*=======================================================================
*
*     Subroutine FindPT - Finds indices Pidx and Tidx for Log K array
*                         using specified P, T
*
*=======================================================================

      SUBROUTINE FindPT (P, T, Pidx, Tidx, PTfind, PTerr, isoerr)

      LOGICAL    PTfind

      INTEGER    Pidx, Tidx

      CHARACTER  PTerr, isoerr

      REAL*8     P, T

**** Assign array index integer for pressure (Pidx).

      If (P.le.0.0D0) then
         Pidx = 1
         PTerr = 'E'
      Else If (P.gt.0.0D0.and.P.lt.1.0D0) then
         Pidx = 1
         PTerr = 'W'
      Else If (P.ge.1.0D0.and.P.lt.1.013D0) then
         Pidx = 1
      Else If (P.ge.1.013D0.and.P.lt.4.758D0) then
         Pidx = 2
      Else If (P.ge.4.758D0.and.P.lt.15.537D0) then
         Pidx = 3
      Else If (P.ge.15.537D0.and.P.lt.39.737D0) then
         Pidx = 4
      Else If (P.ge.39.737D0.and.P.lt.85.839D0) then
         Pidx = 5
      Else If (P.ge.85.839D0.and.P.lt.165.212D0) then
         Pidx = 6
      Else If (P.ge.165.212D0.and.P.lt.200.0D0) then
         Pidx = 7
      Else If (P.ge.200.0D0.and.P.lt.400.0D0) then
         Pidx = 8
      Else If (P.ge.400.0D0.and.P.lt.600.0D0) then
         Pidx = 9
      Else If (P.ge.600.0D0.and.P.lt.800.0D0) then
         Pidx = 10
      Else If (P.ge.800.0D0.and.P.lt.1000.0D0) then
         Pidx = 11
      Else If (P.ge.1000.0D0.and.P.lt.1500.0D0) then
         Pidx = 12
      Else If (P.ge.1500.0D0.and.P.lt.2000.0D0) then
         Pidx = 13
      Else If (P.ge.2000.0D0.and.P.lt.2500.0d0) then
         Pidx = 14
      Else If (P.ge.2500.0D0.and.P.lt.3000.0D0) then
         Pidx = 15
      Else If (P.ge.3000.0D0.and.P.lt.3500.0D0) then
         Pidx = 16
      Else If (P.ge.3500.0D0.and.P.lt.4000.0D0) then
         Pidx = 17
      Else If (P.ge.4000.0D0.and.P.lt.5000.0D0) then
         Pidx = 18
      Else If (P.ge.5000.0D0.and.P.lt.6000.0D0) then
         Pidx = 18
         PTerr = 'W'
      Else If (P.ge.6000.0D0) then
         Pidx = 18
         PTerr = 'E'
      Endif

**** Assign array index integer for temperature (Tidx).

      If (T.lt.0.01D0) then
         PTerr = 'E'
         Tidx = 3
         If (Pidx.gt.17) Pidx = 17
      Else If (T.ge.0.01D0.and.T.le.50.0D0)  then
         Tidx = 3
         If (Pidx.gt.17) Pidx = 17
      Else If (T.gt.50.0D0.and.T.le.100.0D0) then
         Tidx = 4
         If (P.lt.1.013D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
      Else If (T.gt.100.0D0.and.T.le.150.0D0) then
         Tidx = 5
         If (P.lt.4.758D0)   Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.1.013D0)   isoerr = 'E'
      Else If (T.gt.150.0D0.and.T.le.200.0D0) then
         Tidx = 6
         If (P.lt.15.537D0)  Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.4.758D0)   isoerr = 'E'
      Else If (T.gt.200.0D0.and.T.le.250.0D0) then
         Tidx = 7
         If (P.lt.39.737D0)  Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.15.537D0)   isoerr = 'E'
      Else If (T.gt.250.0D0.and.T.le.300.0D0) then
         Tidx = 8
         If (P.lt.85.839D0)  Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.39.737D0)  isoerr = 'E'
      Else If (T.gt.300.0D0.and.T.le.350.0D0) then
         Tidx = 9
         If (P.lt.165.212D0) Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.85.839D0)  isoerr = 'E'
      Else If (T.gt.350.0D0.and.T.le.400.0D0) then
         Tidx = 10
         If (P.ge.165.212D0.and.P.lt.200.0D0) then
            Pidx = Pidx + 2
            isoerr = 'W'
         Endif
         If (P.ge.200.0D0.and.P.lt.400.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.165.212D0) isoerr = 'E'
      Else If (T.gt.400.0D0.and.T.le.410.0D0) then
         Tidx = 11
         If (P.ge.200.0D0.and.P.lt.400.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.200.0D0) isoerr = 'E'
      Else If (T.gt.410.0D0.and.T.le.420.0D0) then
         Tidx = 12
         If (P.ge.200.0D0.and.P.lt.400.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.200.0D0) isoerr = 'E'
      Else If (T.gt.420.0D0.and.T.le.430.0D0) then
         Tidx = 13
         If (P.ge.200.0D0.and.P.lt.400.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.200.0D0) isoerr = 'E'
      Else If (T.gt.430.0D0.and.T.le.450.0D0) then
         Tidx = 14
         If (P.ge.400.0D0.and.P.lt.600.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.400.0D0) isoerr = 'E'
      Else If (T.gt.450.0D0.and.T.le.460.0D0) then
         Tidx = 15
         If (P.ge.400.0D0.and.P.lt.600.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.400.0D0) isoerr = 'E'
      Else If (T.gt.460.0D0.and.T.le.470.0D0) then
         Tidx = 16
         If (P.ge.400.0D0.and.P.lt.600.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.400.0D0) isoerr = 'E'
      Else If (T.gt.470.0D0.and.T.le.480.0D0) then
         Tidx = 17
         If (P.ge.400.0D0.and.P.lt.600.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.400.0D0) isoerr = 'E'
      Else If (T.gt.480.0D0.and.T.le.490.0D0) then
         Tidx = 18
         If (P.ge.400.0D0.and.P.lt.600.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.400.0D0) isoerr = 'E'
      Else If (T.gt.490.0D0.and.T.le.500.0D0) then
         Tidx = 19
         If (P.ge.600.0D0.and.P.lt.800.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.600.0D0) isoerr = 'E'
      Else If (T.gt.500.0D0.and.T.le.550.0D0) then
         Tidx = 20
         If (P.ge.600.0D0.and.P.lt.800.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.600.0D0) isoerr = 'E'
      Else If (T.gt.550.0D0.and.T.le.600.0D0) then
         Tidx = 21
         If (P.ge.800.0D0.and.P.lt.1000.0D0)  then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.800.0D0) isoerr = 'E'
      Else If (T.gt.600.0D0.and.T.le.700.0D0) then
         Tidx = 21
         If (PTerr.eq.'N') PTerr = 'W'
         If (P.ge.800.0D0.and.P.lt.1000.0D0)  then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.800.0D0) isoerr = 'E'
      Else If (T.gt.700.0D0) then
         Tidx = 21
         PTerr = 'E'
      Endif

      PTfind = .true.

      Return
      End

*=======================================================================
*
*     Subroutine FindLK - Interpolates LogK from Log K array using
*                         specified P, T, and indices found using
*                         subroutine findPT, Pidx and Tidx
*
*=======================================================================

***   Variables in terms of X, Y, f(X,Y):
***                     (X, Y, Xindex, Yindex, table of f(Xindex,Yindex), f(X,Y))
      Subroutine FindLK (P, T, Pidx, Tidx, LKarra, LogK)

      INTEGER    Pidx, Tidx, i, j, k, l

      REAL*8     P, T, PVal(20), TVal(21), LKarra (20,21)
      REAL*8     sum, lk(3), LogK

      DATA PVal  / 1.0d0, 1.013d0, 4.758, 15.537, 39.737, 85.839,
     +             165.212, 200, 400, 600, 800, 1000, 1500, 2000,
     +             2500, 3000, 3500, 4000, 4500, 5000 /
      DATA TVal  / 0.01, 25, 50, 100, 150, 200, 250, 300, 350, 400, 410,
     +             420, 430, 450, 460, 470, 480, 490, 500, 550, 600 /

      lk(1) = 0.0d0
      lk(2) = 0.0d0
      lk(3) = 0.0d0
      logK  = 0.0d0

      l = 1
      Do 500, i = Tidx - 2, Tidx
         Do 400, j = Pidx + 2, Pidx, -1
            sum = LKarra(j,i)
            Do 300, k = Pidx + 2, Pidx, -1
               If (j.eq.k) goto 300
               sum = sum*(P-PVal(k))/(PVal(j)-PVal(k))
300        Continue
           lK(l) = lK(l) + sum
400      Continue
         l = l + 1
500   Continue

      k = 1
      Do 700, i = Tidx - 2, Tidx
         sum = lk(k)
         Do 600, j = Tidx - 2, Tidx
            If (i.eq.j) goto 600
            sum = sum*(T-TVal(j))/(TVal(i)-TVal(j))
600      Continue
         logK = logK + sum
         k = k + 1
700   Continue

      write (*,'(F8.3)') logk

      Return
      END

*=======================================================================