*=======================================================================
*
*     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

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

      PARAMETER (Pdim = 20)
      PARAMETER (Tdim = 21)

      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,20
         Read (inp1,*) (LKarra (i,j), j = 1,21)
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,'(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.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.le.5000.0D0) then
         Pidx = 18
      Else If (P.gt.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)

****     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.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    Pdim, Tdim, Pidx, Tidx, i, j, k, l

      PARAMETER (Pdim = 20)
      PARAMETER (Tdim = 21)

      REAL*8     P, T, PVal(Pdim), TVal(Tdim), LKarra (Pdim,Tdim)
      REAL*8     sum, lk(3), LogK

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

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