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

      LOGICAL  TWARN

      INTEGER  i, iout1

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

**** All temperatures considered (not passed to subroutine, for ouput here only)

      DATA TARRAY / 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/

**** Quartz data to 350°C:
      DATA LVARRA /-4.204D0, -3.790D0, -3.493D0, -3.252D0, -3.045D0,
     +             -2.861D0, -2.694D0, -2.542D0, -2.403D0, -2.277D0,
     +             -2.165D0, -2.068D0, -1.992D0, -1.951D0, -1.984D0/

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

**** Test at every 5°C from 10 to 405°C
      DO 100, I = 5,405,5
         TESTT (I/5) = I
         CALL LKLVEX(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,15)
      WRITE (IOUT1,113) 'LOGK,EXPT.  ', (LVARRA(IWRITE), IWRITE = 1,15)
113   FORMAT (A12,15(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 LKLVEX. Interpolates log K(s) at LV saturation
*
*=======================================================================

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

      LOGICAL    TWARN

      INTEGER    I, IT, J, MAXTIX, iout1

      REAL*8     TEMPC, LGKDAT(15), TDATA(15), LOGK, MAXT, SUM

      DATA TDATA  / 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/

**Max TempC value and array index for that value (Tdata(15) = 350 currently)
      iout1 = 11
      MAXTIX = 15
      MAXT = 350.D0
      LOGK = 0.D0
      IT   = 0

**Compute desired logK value (log K) for the given TempC-value (T).
** Find where the desired TempC value lies in the TempC data array.
      DO 150, I = 1,14
         IF (TEMPC.GT.TDATA(I).AND.TEMPC.LE.TDATA(I+1)) IT = I
150   CONTINUE
      IF (TEMPC.GT.TDATA(14)) IT = 14

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

      RETURN
      END