PROGRAM GEOTHERMOMETER C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C CHARACTER*5 LABELS(11) DATA LABELS 1 / 'TiO2 ', 'Al2O3', 'Cr2O3', 'V2O3 ', 'Fe2O3', 'FeO ', 2 'MgO ', 'MnO ', 'ZnO ', 'CoO ', 'NiO ' / DOUBLE PRECISION WT(11) DATA WT 1 / 79.8988D0, 101.9612D0, 151.9902D0, 149.8822D0, 2 159.6922D0, 71.8464D0, 40.3114D0, 70.9374D0, 3 81.3694D0, 74.9324D0, 74.7094D0 / DOUBLE PRECISION CATIONS(11) DATA CATIONS 1 / 1.0D0, 2.0D0, 2.0D0, 2.0D0, 2.0D0, 1.0D0, 1.0D0, 1.0D0, 2 1.0D0, 1.0D0, 1.0D0 / DOUBLE PRECISION CHARGE(11) DATA CHARGE 1 / 4.0D0, 3.0D0, 3.0D0, 3.0D0, 3.0D0, 2.0D0, 2.0D0, 2.0D0, 2 2.0D0, 2.0D0, 2.0D0 / DOUBLE PRECISION OXYGENS(11) DATA OXYGENS 1 / 2.0D0, 3.0D0, 3.0D0, 3.0D0, 3.0D0, 1.0D0, 1.0D0, 1.0D0, 2 1.0D0, 1.0D0, 1.0D0 / DOUBLE PRECISION MOLES(11) LOGICAL BIG_LOOP, FIRST, CONV_NEWTON, LSTOP CHARACTER*10 PHASE CHARACTER*20 INPUT_DATA C COMMON /PARAM/ H11, S11, H55, S55, HX, SX, HEX, SEX, H24, S24, 1 H25, S25, W11, W22, W55, W15, W1P5, W15P, W1P5P, W14, W1P4, 2 W24U, W2P4U, W2P5U, W25PU, WOCT, WTET, W45, W45P, W4U5PU C COMMON /PARAMCR/ H33, S33, H23, S23, W33, W13, W1P3, W13P, 1 W1P3P, W2P3U, W23PU, W34, W3P4, W3PU4U, W35, W3P5, W35P, 2 W3P5P C C Initialize constants C BIG_LOOP = .TRUE. R1 = 0.001987 P = 1.0 C C Solution parameters from Sack and Ghiorso (1990b) (kcals) C H11 = -8.7 S11 = 0.0 W11 = 4.5 W14 = 20.8 W1P4 = 12.4 W15 = 10.0 W1P5 = 14.4 W15P = 11.7 W1P5P = 7.0 H21 = -88.4 S21 = -0.0098 W22 = 3.6 H24 = 6.55 S24 = 0.0000 W24U = 12.6 W2P4U = 10.9 H25 = 8.05 S25 = 0.0 W25PU = 15.3 W2P5U = 14.4 W45 = 6.0 W45P = 2.0 W4U5PU = 1.8 H55 = 6.25 S55 = 0.0 W55 = 0.0 HEX = -3.6 SEX = 0.0 HX = 2.4 SX = 0.0 WOCT = 2.0 WTET = 2.0 C H33 = -20.0 S33 = 0.0 H23 = 0.0 S23 = 0.0 W33 = 7.0 W13 = 10.0 W1P3 = 15.2 W13P = 11.3 W1P3P = 5.9 W2P3U = 10.0 W23PU = 9.7 W34 = 11.5 W3P4 = 10.0 W3PU4U = 10.4 W35 = 0.0 W3P5 = 10.0 W35P = 8.0 W3P5P = 0.0 C C Input the composition of a coexisting spinel and olivine C SUM_OXY = 0.0 WRITE(*,1000) 1000 FORMAT(/,1X,'Input the composition of the spinel phase: ') DO I=1,11 WRITE(*,1001) LABELS(I) 1001 FORMAT(6X,'Wt % of ',A5,' [default is 0.00]: ',$) READ(*,'(A20)') INPUT_DATA IF(INPUT_DATA.NE.' ') THEN READ(INPUT_DATA,*) MOLES(I) MOLES(I) = MOLES(I)/WT(I) SUM_OXY = SUM_OXY + MOLES(I)*OXYGENS(I) ELSE MOLES(I) = 0.0 END IF END DO C SUM_CAT = 0.0 IF(SUM_OXY.GT.0.0) THEN DO I=1,11 SUM_CAT = SUM_CAT + MOLES(I)*CATIONS(I) END DO ELSE STOP '%GTH_F_ZEROSP, The spinel must contain something.' END IF C SUM_CAT = 3.0/SUM_CAT DO I=1,11 MOLES(I) = MOLES(I)*CATIONS(I)*SUM_CAT END DO C SUM_CHARGE = 0.0 DO I=1,11 SUM_CHARGE = SUM_CHARGE + CHARGE(I)*MOLES(I) END DO C SUM_CHARGE = SUM_CHARGE - 3.0*MOLES(5) - 2.0*MOLES(6) FE3 = 8.0 - SUM_CHARGE - 2.0*(MOLES(5)+MOLES(6)) FE2 = MOLES(5) + MOLES(6) - FE3 C IF(FE3.LT.0.0) STOP '%GTH_F_NEGFE3, Negative Fe3+ in spinel' IF(FE2.LT.0.0) STOP '%GTH_F_NEGFE2, Negative Fe2+ in spinel' MOLES(5) = FE3 MOLES(6) = FE2 C SUM = MOLES(6) + MOLES(7) + MOLES(8) + MOLES(9) + MOLES(10) + 1 + MOLES(11) X2 = MOLES(7)*SUM/(SUM-MOLES(8)-MOLES(9)-MOLES(10)-MOLES(11)) X3 = (MOLES(3)+MOLES(4))/2.0 X4 = MOLES(1) X5 = MOLES(5)/2.0 WRITE(*,1003) (1.0-X2-X3-X4-X5), X2, X3, X4, X5 1003 FORMAT(6X,'Mole fractions of FeAl2O4 (X1):',F10.6,/, 1 6X,'Mole fractions of MgAl2O4 (X2):',F10.6,/, 2 6X,'Mole fractions of FeCr2O4 (X3):',F10.6,/, 3 6X,'Mole fractions of Fe2TiO4 (X4):',F10.6,/, 4 6X,'Mole fractions of Fe3O4 (X5):',F10.6) C WRITE(*,1004) 1004 FORMAT(/,1X,'Input the composition of the olivine phase: ') SUM_CAT = 0.0 DO I=6,7 WRITE(*,1001) LABELS(I) READ(*,'(A20)') INPUT_DATA IF(INPUT_DATA.NE.' ') THEN READ(INPUT_DATA,*) MOLES(I) MOLES(I) = MOLES(I)/WT(I) SUM_CAT = SUM_CAT + MOLES(I) ELSE MOLES(I) = 0.0 END IF END DO IF(SUM_CAT.GT.0.0) THEN XFO = MOLES(7)/SUM_CAT ELSE STOP '%GTH_F_ZEROOL, The olivine must contain something.' END IF WRITE(*,1005) XFO, (1.0-XFO) 1005 FORMAT(6X,'Mole fractions of Mg2SiO4 (XFO):',F10.6,/, 1 6X,'Mole fractions of Fe2SiO4 (XFO):',F10.6) C C Calculate a temperature of coexistence C CONV_NEWTON = .FALSE. FIRST = .TRUE. T = 500.0 + 273.15 C DO WHILE(.NOT.CONV_NEWTON) C C Get ordering parameters at this T (K) C T, X2, X3, X4 and X5 are input, other variables are returned C CALL ORD_SPN(T, X2, X3, X4, X5, S1, S2, S3, S4, 1 XMG2TET, XFE2TET, XAL3TET, XFE3TET, XCR3TET, XMG2OCT, 2 XFE2OCT, XAL3OCT, XFE3OCT, XTI4OCT, XCR3OCT, LSTOP) IF(LSTOP) THEN WRITE(*,1006) 1006 FORMAT(1X,78('*'),/, 1 1X,'%SPN_F_NCONV, Failed to converge in ordering parameter ', 2 'subroutine.',/, 2 1X,78('*')) STOP END IF C G11 = H11 - T*S11 G33 = H33 - T*S33 G55 = H55 - T*S55 GX = HX - T*SX GEX = HEX - T*SEX G23 = H23 - T*S23 G24 = H24 - T*S24 G25 = H25 - T*S25 G0FO = GIBBS('FORSTERITE', T, P) G0FA = GIBBS('FAYALITE ', T, P) G21 = H21 - T*S21 C DG = 0.5*(G0FO - G0FA) + R1*T*LOG(XFO/(1.0-XFO)) 1 + 2.43*(1.0-2.0*XFO) - ( 2 + R1*T*LOG(XMG2TET/XFE2TET) + W2P5U*X5 - W22*X5 - W1P5*X5 3 + W11*X5 + G25*X5 + WTET*X4 - W24U*X4 + W14*X4 + W2P3U*X3 4 - W22*X3 - W1P3*X3 + W11*X3 + G23*X3 - WTET*X2 - GX*X2/2.0 5 + S4*WTET + S3*WTET + S2*WTET - S1*WTET + S4*W4U5PU 6 - S4*W45P + S3*W3PU4U - S3*W3P4 - S4*W2P5U + S2*W2P4U 7 - S3*W2P3U - S4*W25PU + W24U - S3*W23PU + S4*W22 + S3*W22 8 - S2*W22 + S4*W1P5 - S2*W1P4 + S3*W1P3 + S4*W15P - W14 9 + S3*W13P - S4*W11 - S3*W11 + S2*W11 - GX*S4/2.0 A - GEX*S4/2.0 - G25*S4 + G24*S4 - GX*S3/2.0 - GEX*S3/2.0 B + G24*S3 - G23*S3 - GX*S2/2.0 - GEX*S2/2.0 + G24*S2 C + GX*S1/2.0 + GX/2.0 + GEX/2.0 + G24 + G21 ) IF(FIRST) THEN FOLD = DG TOLD = T T = T + 50.0 FIRST = .FALSE. ELSE F = DG DF = (F-FOLD)/(T-TOLD) FOLD = F TOLD = T T = T - F/DF CONV_NEWTON = ABS(T-TOLD).LT.1.0D-4 WRITE(*,1007) T-273.15, F 1007 FORMAT(1X,'%GTH_I_NEWTON, Newton iteration: T , delta G: ', 1 F12.6, F14.6) END IF END DO WRITE(*,1008) T-273.15 1008 FORMAT(1X,/,78('*'),/, 1 1X,'Final equilibration temperature: ',F8.2,/, 2 1X,78('*'),/) WRITE(*,1009) S1, S2, S3, S4, XMG2TET, XFE2TET, XAL3TET, 1 XFE3TET, XCR3TET, XMG2OCT, XFE2OCT, XAL3OCT, 2 XFE3OCT, XCR3OCT, XTI4OCT 1009 FORMAT(1X,'Spinel ordering parameters and site mole fractions:' 1 ,/,1X,'s1,s2,s3,s4: ',4F10.6 2 ,/,1X,'X Mg2+ tet: ', F10.6 3 ,/,1X,'X Fe2+ tet: ', F10.6 4 ,/,1X,'X Al3+ tet: ', F10.6 5 ,/,1X,'X Fe3+ tet: ', F10.6 6 ,/,1X,'X Cr3+ tet: ', F10.6 7 ,/,1X,'X Mg2+ oct: ', F10.6 8 ,/,1X,'X Fe2+ oct: ', F10.6 9 ,/,1X,'X Al3+ oct: ', F10.6 A ,/,1X,'X Fe3+ oct: ', F10.6 B ,/,1X,'X Cr3+ oct: ', F10.6 C ,/,1X,'X Ti4+ oct: ', F10.6) C STOP END SUBROUTINE ORD_SPN(T, X2, X3, X4, X5, S1, S2, S3, S4, 1 XMG2TET, XFE2TET, XAL3TET, XFE3TET, XCR3TET, XMG2OCT, 2 XFE2OCT, XAL3OCT, XFE3OCT, XTI4OCT, XCR3OCT, LSTOP) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C COMMON /PARAM/ H11, S11, H55, S55, HX, SX, HEX, SEX, H24, S24, 1 H25, S25, W11, W22, W55, W15, W1P5, W15P, W1P5P, W14, W1P4, 2 W24U, W2P4U, W2P5U, W25PU, WOCT, WTET, W45, W45P, W4U5PU C COMMON /PARAMCR/ H33, S33, H23, S23, W33, W13, W1P3, W13P, 1 W1P3P, W2P3U, W23PU, W34, W3P4, W3PU4U, W35, W3P5, W35P, 2 W3P5P C LOGICAL CONVERGED, LSTOP C LSTOP = .FALSE. C C T is input in Kelvins, X2, X4 and X5 are input compositions C Everything else is output C G11 = H11 - T*S11 G33 = H33 - T*S33 G55 = H55 - T*S55 GX = HX - T*SX GEX = HEX - T*SEX G23 = H23 - T*S23 G24 = H24 - T*S24 G25 = H25 - T*S25 C H5PU5U = H55 + 2.0*H24 - H25 + (W24U-W14) + (W4U5PU-W45P) 1 - (W25PU-W15P) S5PU5U = S55 + 2.0*S24 - S25 W5U5PU = W55 + (W2P5U+W25PU+W11) - (W1P5+W15P+W22) W4U5U = W45 + H25 + (W24U+W2P5U+W11) - (W14+W1P5+W22) G5PU5U = H5PU5U - T*S5PU5U W4U5U = W4U5U - T*S25 H22 = 2.0*H24+H11+(W24U+W2P4U+W11)-(W14+W1P4+W22) S22 = 2.0*S24+S11 G22 = H22 - T*S22 C GS1 = 0.5*(-WOCT + W24U - W14 + 0.5*GX + 0.5*GEX + G24) GS2 = W11 + G11 GS3 = W13P - W13 + G33 GS4 = W15P - W15 + G55 GX2X2 = -0.25*(WTET + WOCT + GX) GX3X3 = - W13 GX4X4 = -W14 GX5X5 = -W15 GS1S1 = 0.25*(-WTET - WOCT + GX) GS2S2 = -W11 GS3S3 = -W33 GS4S4 = -W55 GS1S2= 0.5*(WTET - W22 + W11 + WOCT -GX) GS1S3 = 0.5*(WTET + WOCT - W2P3U - W23PU + W22 + W1P3 + W13P 1 - W11 - GX) GS1S4 = 0.5*(WTET + WOCT + W22 - W11 - W2P5U - W25PU + W1P5 + W15P 1 - GX) GS2S3 = W1P3P - W1P3 - W13P + W13 GS2S4 = W1P5P - W1P5 - W15P + W15 GS3S4 = W3P5P - W3P5 - W35P + W35 GX2X3 = 0.5*(W2P3U - W22 -W1P3 + W11 + 2.0*G23) GX2X4 = 0.5*(WTET - W24U + W14 + 0.5*GX - 0.5*GEX + G24) GX2X5 = 0.5*(-W22 + W11 + W2P5U - W1P5 + 2.0*G25) GX3X4 = W34 - W14 - W13 GX3X5 = W35 - W15 -W13 GX4X5 = W45 - W15 - W14 GX2S1 = 0.5*(WOCT - WTET) GX2S2 = 0.5*(WTET - W22 + W11 - WOCT + 2.0*W2P4U - 2.0*W1P4 - GEX 1 + 2.0*G24) GX2S3 = 0.5*(WTET - WOCT + 2.0*W3PU4U - 2.0*W3P4 - W2P3U - W23PU 1 + W22 + W1P3 + W13P - W11 - GEX + 2.0*G24 - 2.0*G23) GX2S4 = 0.5*(WTET - WOCT - W11 + W22 + 2.0*W4U5PU - 2.0*W45P 1 - W2P5U - W25PU + W1P5 + W15P - GEX - 2.0*G25 + 2.0*G24) GX3S1 = 0.5*(W2P3U - W22 - W1P3 + W11) GX3S2 = W1P3 - W13 - W11 GX3S3 = W33 - W13P + W13 GX3S4 = W35P - W35 - W15P + W15 GX4S1 = 0.5*(WTET - W24U + W14 - 0.5*GX + 0.5*GEX - G24) GX4S2 = -W11 + W1P4 - W14 GX4S3 = W3P4 - W34 - W13P + W13 GX4S4 = W45P - W45 - W15P + W15 GX5S1 = 0.5*(-W22 + W11 + W2P5U - W1P5) GX5S2 = -W11 + W1P5 - W15 GX5S3 = W3P5 - W35 - W13P + W13 GX5S4 = W55 - W15P + W15 IF((X2.EQ.1.0D0).AND.(X3.EQ.0.0D0).AND.(X4.EQ.0.0D0).AND. 1 (X5.EQ.0.0D0)) THEN S2 = 1.0 S1 = 2.0*S2 - 1.0 S3 = 0.0 S4 = 0.0 ELSE IF((X2.EQ.1.0D0).AND.(X3.EQ.1.0D0).AND.(X4.EQ.0.0D0).AND. 1 (X5.EQ.0.0D0)) THEN S2 = 0.0 S3 = 1.0 S4 = 0.0 S1 = 2.0*S3 - 1.0 ELSE IF((X2.EQ.1.0D0).AND.(X3.EQ.0.0D0).AND.(X4.EQ.0.0D0).AND. 1 (X5.EQ.1.0D0)) THEN S2 = 0.0 S3 = 0.0 S4 = 0.0 S1 = 2.0*S4 - 1.0 ELSE S1 = 0.5*X2 S2 = X2 S3 = X3 S4 = 0.0 END IF R1 = .001987 J = 0 CONVERGED = .FALSE. RELERR = 1.0D-6 SMALL = 1.0D-12 DO WHILE(.NOT.CONVERGED.AND.(J.LE.1000)) S1OLD = S1 S2OLD = S2 S4OLD = S4 C IF((X2.EQ.1.0D0).AND.(X3.EQ.0.0D0).AND.(X4.EQ.0.0D0).AND. 1 (X5.EQ.0.0D0)) THEN S1 = 2.0*S2 - 1.0 ELSE IF((X2.EQ.1.0D0).AND.(X3.EQ.1.0D0).AND.(X4.EQ.0.0D0).AND. 1 (X5.EQ.0.0D0)) THEN S1 = 2.0*S3 - 1.0 ELSE IF((X2.EQ.1.0D0).AND.(X3.EQ.0.0D0).AND.(X4.EQ.0.0D0).AND. 1 (X5.EQ.1.0D0)) THEN S1 = 2.0*S4 - 1.0 ELSE A = 2.0*GS1 + 4.0*GS1S1*S1 + 2.0*GX2S1*X2 + 2.0*GX4S1*X4 1 + 2.0*GX5S1*(X5) + 2.0*GS1S2*S2 + 2.0*GS1S4*S4 2 + 2.0*GX3S1*X3 + 2.0*GS1S3*S3 A = EXP(A/(R1*T)) C1 = (2.0*X2*X4 + 2.0*X2*S2 + 2.0*X2*S3 + 2.0*X2*S4 - X2**2) 1 - A*(2.0*X2 - X2**2 - 2.0*X2*S2 - 2.0*X2*S3 - 2.0*X2*S4) B1 = (-2.0*X4 - 2.0*S2 - 2.0*S3 - 2.0*S4) 1 - A*(2.0 - 2.0*S2 - 2.0*S3 - 2.0*S4) A1 = 1.0 - A Q = SQRT(B1**2 -4.0*A1*C1) IF(A1.NE.0.0) THEN S1 = (-B1 -Q) / (2.0*A1) ELSE S1 = -C1/B1 END IF END IF C IF((X2.EQ.1.0D0).AND.(X3.EQ.0.0D0).AND.(X4.EQ.0.0D0).AND. 1 (X5.EQ.0.0D0)) THEN A = G22 + W22*(1.0-2.0*S2) A = EXP(A/(R1*T)) A1 = (1.0 - A) B1 = -2.0 - A C1 = 1.0 Q = SQRT(B1**2-4.0*A1*C1) IF(A1.NE.0.0) THEN S2 = (-B1-Q)/(2.0*A1) ELSE S2 = -C1/B1 END IF S1 = 2.0*S2 - 1.0 ELSE A = GS2 + 2.0*GS2S2*S2 + GX2S2*X2 + GX3S2*X3 + GX4S2*X4 1 + GX5S2*(X5) + GS1S2*S1 +GS2S3*S3 + GS2S4*S4 A = EXP(A/(R1*T)) C1 = (1.0 - 0.5*X2 + 0.5*S1 - S3 - S4 -X3 + 0.5*X2*X3 1 - 0.5*S1*X3 + S3*X3 ) C1 = C1 + (S4*X3 - X4 + 0.5*X2*X4 - 0.5*S1*X4 + S3*X4 1 + S4*X4 - X5 + 0.5*X2*X5) C1 = C1 + (-0.5*S1*X5 + S3*X5 + S4*X5) C1 = C1 - A*(X4 - 0.5*X2 - 0.5*S1 + S3 + S4) C1 = C1 - A*(-X4*X3 + 0.5*X2*X3 + 0.5*S1*X3 - S3*X3 - S4*X3) C1 = C1 - A*(-X4*X4 + 0.5*X2*X4 + 0.5*S1*X4 - S3*X4 - S4*X4) C1 = C1 - A*(-X4*X5 + 0.5*X2*X5 + 0.5*S1*X5 - S3*X5 - S4*X5) B1 = (-1.0 + X3 + X4 + X5 -1.0 + 0.5*X2 - 0.5*S1 + S3 + S4) B1 = B1 - A*(1.0 - X3 - X4 - X5 + X4 - 0.5*X2 - 0.5*S1 1 + S3 + S4) A1 = 1.0 - A Q = SQRT(B1**2-4.0*A1*C1) IF(A1.NE.0.0) THEN S2 = (-B1-Q) / (2.0*A1) ELSE S2 = -C1/B1 END IF END IF C IF((X2.EQ.1.0D0).AND.(X3.EQ.0.0D0).AND.(X4.EQ.0.0D0).AND. 1 (X5.EQ.1.0D0)) THEN A = G5PU5U + W5U5PU*(1.0-2.0*S4) A = EXP(A/(R1*T)) A1 = (1.0 - A) B1 = -2.0 - A C1 = 1.0 Q = SQRT(B1**2-4.0*A1*C1) IF(A1.NE.0.0) THEN S4 = (-B1-Q)/(2.0*A1) ELSE S4 = -C1/B1 END IF S1 = 2.0*S4 - 1.0 ELSE A = GS4 + 2.0*GS4S4*S4 + GX2S4*X2 + GX3S4*X3 + GX4S4*X4 1 + GX5S4*(X5) + GS1S4*S1 + GS2S4*S2 + GS3S4*S3 A = EXP(A/(R1*T)) C1 = (X5 - 0.5*X2*X5 + 0.5*S1*X5 - S2*X5 -S3*X5) -A*(X4*X5 1 - 0.5*X2*X5 - 0.5*S1*X5 + S2*X5 + S3*X5) B1 = (-X5 -1.0 + 0.5*X2 - 0.5*S1 + S2 + S3) 1 - A*(X5 + X4 - 0.5*X2 - 0.5*S1 + S2 + S3) A1 = 1.0- A Q = SQRT(B1**2 -4.0*A1*C1) IF(A1.NE.0.0) THEN S4 = (-B1 -Q) / (2.0*A1) ELSE S4 = -C1/B1 END IF END IF C J = J + 1 IF(J.GT.1000) LSTOP = .TRUE. CONVERGED = (ABS(S1-S1OLD) .LE. RELERR*ABS(S1)) .OR. 1 (ABS(S1-S1OLD) .LE. SMALL) CONVERGED = CONVERGED .AND. ( 1 (ABS(S2-S2OLD) .LE. RELERR*ABS(S2)) .OR. 2 (ABS(S2-S2OLD) .LE. SMALL) ) CONVERGED = CONVERGED .AND. ( 1 (ABS(S4-S4OLD) .LE. RELERR*ABS(S4)) .OR. 2 (ABS(S4-S4OLD) .LE. SMALL) ) C END DO C XMG2TET = 0.5*(X2 + S1) XFE2TET = (X4 - 0.5*X2 - 0.5*S1 + S2 + S3 + S4) XAL3TET = (1.0 - X3 - X4 - X5 - S2) XFE3TET= X5 - S4 XCR3TET = X3 - S3 XMG2OCT = 0.25*(X2 - S1) XFE2OCT = 0.25*(2.0 - X2 + S1 - 2.0*S2 - 2.0*S3 - 2.0*S4) XAL3OCT = 0.5*(1.0 - X3 - X4 - X5 + S2) XFE3OCT = 0.5*(X5 + S4) XCR3OCT = 0.5*(X3 + S3) XTI4OCT = 0.5*X4 C RETURN END DOUBLE PRECISION FUNCTION GIBBS(PHASE, T, P) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C DOUBLE PRECISION K0,K1,K2,K3 CHARACTER*10 PHASE C IF(PHASE.EQ.'FORSTERITE') THEN HS = -2174420.0 ! J Berman (1988) SS = 94.010 ! J Berman (1988) VS = 4.366 ! J Berman (1988) K0 = 238.64 ! J Berman and Brown (1985) K1 = -20.013D2 ! J K2 = 0.0 K3 = -11.624D7 ! J V1 = -0.791D-6 ! Berman (1988) To 99 Kbars at V2 = 1.351D-12 ! Berman (1988) about 1000 C V3 = 29.464D-6 ! Berman (1988) V4 = 88.633D-10 ! Berman (1988) ELSE IF(PHASE.EQ.'FAYALITE ') THEN HS = -1479360. ! J Berman (1988) SS = 150.930 ! J Berman (1988) VS = 4.630 ! J Berman (1988) K0 = 248.93 ! J Berman and Brown (1985) K1 = -19.239D2 ! J K2 = 0.0 K3 = -13.910D7 ! J V1 = -0.730D-6 ! Berman (1988) To 73 Kbars at V2 = 0.0 ! Berman (1988) about 1225 C V3 = 26.546D-6 ! Berman (1988) V4 = 79.482D-10 ! Berman (1988) ELSE STOP 'Incorrect PHASE designator in call to GIBBS.' END IF C TR = 298.15D0 PR = 1.0D0 C HS = HS + K0*(T-TR) + 2.0*K1*(SQRT(T)-SQRT(TR)) 1 - K2*(1.0/T-1.0/TR) - 0.5*K3*(1.0/T**2-1.0/TR**2) SS = SS + K0*LOG(T/TR) 1 - 2.0*K1*(1.0/SQRT(T)-1.0/SQRT(TR)) 2 - 0.5*K2*(1.0/T**2-1.0/TR**2) 3 - (1.0/3.0)*K3*(1.0/T**3-1.0/TR**3) VS = VS*((V1/2.0-V2)*(P**2-PR**2)+V2*(P**3-PR**3)/3.0 1 +(1.0-V1+V2+V3*(T-TR)+V4*(T-TR)**2)*(P-PR)) C GS = HS - T*SS + VS GIBBS = GS/4184.0 ! Convert to kcals C RETURN END