*======================================================================= * * Program to test subroutine INTERP * *======================================================================= LOGICAL TWARN INTEGER i 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 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) 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