*=======================================================================
*
*                Program to test subroutine INTERP
*
*=======================================================================

      LOGICAL  TWARN

      INTEGER  i, iout1

      REAL*8 TEMPC, TARRAY(8), LVARRA(8), LOGKAR(100), TESTT(100), MAXT

**** All temperatures considered (not passed to subroutine, for ouput here only)
**** DATA TARRAY  /  25.D0,   50.D0,  100.D0,  150.D0,  200.D0,      250.D0,     300.D0,     350.D0/
**** DATA LVARRA /-4.047D0,-3.644D0,-3.117D0,-2.752D0,-2.465D0, 99999.999D0,99999.999D0,99999.999D0/

      DATA TARRAY  /25.D0, 50.D0, 100.D0, 150.D0, 200.D0, 250.D0,
     +             300.D0, 350.D0/

**** Quartz data to 350°C:
*      DATA LVARRA /-4.047D0,-3.644D0,-3.117D0,-2.752D0,-2.465D0,
*     1           -2.228D0,-2.034D0,-1.913D0/

**** Quartz data up to 300°C.
      DATA LVARRA /-4.047D0,-3.644D0,-3.117D0,-2.752D0,-2.465D0,
     1           -2.228D0,-2.034D0,99999.999D0/

**** Quartz data up to 250°C.
*      DATA LVARRA /-4.047D0,-3.644D0,-3.117D0,-2.752D0,-2.465D0,
*     1           -2.228D0,99999.999D0,99999.999D0/

**** Quartz data up to 200°C.
*      DATA LVARRA /-4.047D0,-3.644D0,-3.117D0,-2.752D0,-2.465D0,
*     1           99999.999D0,99999.999D0,99999.999D0/

**** Quartz data up to 150°C.
*      DATA LVARRA /-4.047D0,-3.644D0,-3.117D0,-2.752D0,99999.999D0,
*     1           99999.999D0,99999.999D0,99999.999D0/

**** Quartz data up to 100°C.
*      DATA LVARRA /-4.047D0,-3.644D0,-3.117D0,99999.999D0,99999.999D0,
*     1           99999.999D0,99999.999D0,99999.999D0/

**** Quartz data up to 50°C.
*      DATA LVARRA /-4.047D0,-3.644D0,99999.999D0,99999.999D0,
*     1           99999.999D0,
*     1           99999.999D0,99999.999D0,99999.999D0/

**** Quartz data up to 25°C.
*      DATA LVARRA /-4.047D0,99999.999D0,99999.999D0,99999.999D0,
*     1           99999.999D0,
*     1           99999.999D0,99999.999D0,99999.999D0/

      IOUT1 = 11
      OPEN (IOUT1,  FILE = 'LOGK.DAT',  STATUS = 'UNKNOWN')
      TWARN = .FALSE.

**** Test at every 10°C from 10 to 350°C
      DO 100, I = 5,460,5
         TESTT (I/5) = I
         CALL LOGKLV(TESTT(I/5), TWARN, MAXT, LVARRA, LOGKAR(I/5))
100   CONTINUE

**** Echo T, log K data
      WRITE (IOUT1,113) 'TEMP(C)     ', (TARRAY(IWRITE), IWRITE = 1,8)
      WRITE (IOUT1,113) 'LOGK,EXPT.  ', (LVARRA(IWRITE), IWRITE = 1,8)
113   FORMAT (A12,9(2X,F9.3))
      WRITE (IOUT1,*)

**** Warn if data not extend to 350°C

      WRITE (IOUT1,133) 'SOLTHERM LOG K DATA PRESENT UP TO',MAXT, '°C.'
133   FORMAT (A33, F6.1, A3)
      WRITE (IOUT1,*)

**** Output results.
      WRITE (IOUT1,'(A8,2X,A8)') 'TEMP(C) ','LOGK    '

      DO 200, I = 5,460,5
         WRITE (IOUT1,153) I, LOGKAR(I/5)
153      FORMAT (I4,2X,F9.3)
200   CONTINUE

      CLOSE (IOUT1)

      END

*=======================================================================
*
*     Subroutine logKlv. Interpolates log K(s) at LV saturation
*
*=======================================================================

      SUBROUTINE LOGKLV(TEMPC, TWARN, MAXT, LGKDAT, LOGK)

      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.LE.75.D0.AND.TEMPC.GE.0.1D0) THEN
            LOGK = LGKDAT(1) + (TEMPC - 25.D0)*
     +             (LGKDAT(2)-LGKDAT(1))/25.D0
         ELSE
            LOGK = 99999.999D0
         ENDIF
      ELSE

**Compute desired logK value (log K) 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