*=======================================================================
*
*     Program TestLK - tests subroutines FindPT and FindLK
*
*=======================================================================

**** PTfind = .false. if Pidx, Tidx not yet found using FindPT (Pidx = Tidx = -100)
**** PTfind = .true.  if Pidx, Tidx already found using FindPT

      LOGICAL    PTfind

      CHARACTER  PTerr, isoerr
      character*20 char

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

      PARAMETER (Pdim = 37)
      PARAMETER (Tdim = 36)

      REAL*8     P, T, LKarra(Pdim, Tdim), 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,Pdim
         Read (inp1,*) (LKarra (i,j), j = 1,Tdim)
10    Continue

      PTfind = .false.

      write (IOUT1,'(4(a6,4x),2x,A6,2(4x,a6))') '    P ', '    T ',
     +       '  Pidx', '  Tidx', 'LogK  ', 'isoerr  ', ' PTerr'

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

      read  (INP2,*,end = 200) P,T

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

      write (IOUT1,'(1x,2(F8.3,2x),2(4x,i2,4x),F12.3,2(7x,A1,2x),F8.3)')
     +      P, T, Pidx, Tidx, logk, isoerr, pterr

      goto 100
200   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).

****     PTerr = 'N' if no warning
****     PTerr = '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
****     PTerr = 'E' if P <= 0.0     or
****                 if P >= 6000.0  or
****                 if T > 700.0    or
****                 if T < 0.01

      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.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
         PTerr = 'W'
      Else If (P.ge.6000.0D0) then
         Pidx = 35
         PTerr = 'E'
      Endif

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

****     isoerr = 'N' if no warning
****     isoerr = 'W' if P/T near L-V saturation or
****                  if P/T near critical isochore, rho H2O < 0.35g/cc
****     isoerr = 'E' if P/T fully in the two-phase region or
****                  if P/T fully above critical isochore.

      If (T.lt.0.01D0) then
         PTerr = '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
            isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.1.013D0)   isoerr = 'E'
      Else If (T.gt.125.0D0.and.T.le.150.0D0) then
         Tidx = 7
         If (P.lt.4.758D0)   Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.2.321D0)   isoerr = 'E'
      Else If (T.gt.150.0D0.and.T.le.175.0D0) then
         Tidx = 8
         If (P.lt.8.919D0)  Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.4.758D0)   isoerr = 'E'
      Else If (T.gt.175.0D0.and.T.le.200.0D0) then
         Tidx = 9
         If (P.lt.15.537D0)  Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.8.919D0)   isoerr = 'E'
      Else If (T.gt.200.0D0.and.T.le.225.0D0) then
         Tidx = 10
         If (P.lt.25.479D0)  Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.15.537D0)   isoerr = 'E'
      Else If (T.gt.225.0D0.and.T.le.250.0D0) then
         Tidx = 11
         If (P.lt.39.737D0)  Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.25.479D0)   isoerr = 'E'
      Else If (T.gt.250.0D0.and.T.le.275.0D0) then
         Tidx = 12
         If (P.lt.59.432D0)  Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.39.737D0)  isoerr = 'E'
      Else If (T.gt.275.0D0.and.T.le.300.0D0) then
         Tidx = 13
         If (P.lt.85.839D0)  Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.59.432D0)  isoerr = 'E'
      Else If (T.gt.300.0D0.and.T.le.325.0D0) then
         Tidx = 14
         If (P.lt.120.458D0) Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.85.839D0)  isoerr = 'E'
      Else If (T.gt.325.0D0.and.T.le.350.0D0) then
         Tidx = 15
         If (P.lt.165.212D0) Then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.120.458D0)  isoerr = '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
            isoerr = 'W'
         Endif
         If (P.ge.200.0D0.and.P.lt.250.0D0) then
            Pidx = Pidx + 2
            isoerr = 'W'
         Endif
         If (P.ge.250.0D0.and.P.lt.300.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 = 17
         If (P.ge.200.0D0.and.P.lt.350.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 = 18
         If (P.ge.350.0D0.and.P.lt.400.0D0)   then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.350.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.350.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.400.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.450.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.450.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.500.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.550.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.550.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.600.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.600.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.650.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.700.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.700.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.700.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.800.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.800.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.850.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.900.0D0) isoerr = '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
            isoerr = 'W'
         Endif
         If (P.lt.900.0D0) isoerr = 'E'
      Else If (T.gt.600.0D0.and.T.le.700.0D0) then
         Tidx = 36
         If (PTerr.eq.'N') PTerr = 'W'
         If (P.ge.900.0D0.and.P.lt.950.0D0)  then
            Pidx = Pidx + 1
            isoerr = 'W'
         Endif
         If (P.lt.800.0D0) isoerr = 'E'
      Else If (T.gt.700.0D0) then
         Tidx = 36
         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    Pdim, Tdim, Pidx, Tidx, i, j, k, l

      PARAMETER (Pdim = 37)
      PARAMETER (Tdim = 36)

      REAL*8     P, T, PVal(Pdim), TVal(Tdim), LKarra (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

      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

      Return
      END

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