SUBROUTINE ACT_RHM(TK, XIL, XGK, XPY, S, T, U, ILMENITE, 1 GEIKIELITE, HEMATITE, PYROPHANITE, LSTOP) C C Program to calculate the activities of ilmenite, geikielite, C Hematite and pyrophanite according to the solution model of C Ghiorso and Sack (1990, Appendix) - extention of Ghiorso (1990) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C DOUBLE PRECISION ILMENITE, ILMENITE_PURE LOGICAL LSTOP, CONVERGED, DEBUG EXTERNAL ORDER DIMENSION WORK(27), AUXVAR(29) C P = 1.0 HIL = 17000.0 SIL = 0.0 HGK = 17000.0 SGK = 0.0 HPY = 17000.0 SPY = 0.0 DWILHM = -19000.0 DWGKHM = -19000.0 DWPYHM = -19000.0 WILHM = 5000.0 WGKHM = 5000.0 WPYHM = 5000.0 DWILGK = 12000.0 DWGKIL = 8000.0 DWILPY = 17000.0 DWPYIL = 17000.0 DWGKPY = 8000.0 DWPYGK = 12000.0 WILGKDIS = 25000.0 WILGKORD = 5000.0 WILPYDIS = 2200.0 WILPYORD = 2200.0 WGKPYDIS = 25000.0 WGKPYORD = 5000.0 C AUXVAR( 1) = TK AUXVAR( 2) = P AUXVAR( 3) = XIL AUXVAR( 4) = XGK AUXVAR( 5) = XPY AUXVAR( 6) = HIL AUXVAR( 7) = SIL AUXVAR( 8) = HGK AUXVAR( 9) = SGK AUXVAR(10) = HPY AUXVAR(11) = SPY AUXVAR(12) = DWILHM AUXVAR(13) = DWGKHM AUXVAR(14) = DWPYHM AUXVAR(15) = WILHM AUXVAR(16) = WGKHM AUXVAR(17) = WPYHM AUXVAR(18) = DWILGK AUXVAR(19) = DWGKIL AUXVAR(20) = DWILPY AUXVAR(21) = DWPYIL AUXVAR(22) = DWGKPY AUXVAR(23) = DWPYGK AUXVAR(24) = WILGKDIS AUXVAR(25) = WILGKORD AUXVAR(26) = WILPYDIS AUXVAR(27) = WILPYORD AUXVAR(28) = WGKPYDIS AUXVAR(29) = WGKPYORD C DGIL = HIL - TK*SIL DGGK = HGK - TK*SGK DGPY = HPY - TK*SPY C WORK(1) = XIL*0.9 ! Initial guess for s WORK(2) = XGK*0.9 ! Initial guess for t WORK(3) = XPY*0.9 ! Initial guess for u ITMAX = 250 DEBUG = .FALSE. CALL VARMET(ORDER, 3, WORK, GSOLN, WORK(4), WORK(7), 1 WORK(10), WORK(13), WORK(16), 3, AUXVAR, ITMAX, DEBUG, 2 MODE) S = WORK(1) T = WORK(2) U = WORK(3) IF(MODE.NE.0) WRITE(*,*) 'Mode (VARMET) = ', MODE C RTLNHM = WPYHM*XPY**2 + WPYHM*XIL*XPY - WILPYDIS*XIL*XPY 1 + WILHM*XIL*XPY + WPYHM*XGK*XPY - WGKPYDIS*XGK*XPY 2 + WGKHM*XGK*XPY - 4.0*DWPYHM*U**2*XPY - 4.0*DGPY*U**2*XPY 3 - 4.0*DWGKPY*T**2*XPY - 4.0*DWGKHM*T**2*XPY 4 - 4.0*DWILPY*S**2*XPY - 4.0*DWILHM*S**2*XPY + WILHM*XIL**2 5 + WILHM*XGK*XIL - WILGKDIS*XGK*XIL + WGKHM*XGK*XIL 6 - 4.0*DWPYIL*U**2*XIL - 4.0*DWPYHM*U**2*XIL 7 - 4.0*DWGKIL*T**2*XIL - 4.0*DWGKHM*T**2*XIL 8 - 4.0*DWILHM*S**2*XIL - 4.0*DGIL*S**2*XIL + WGKHM*XGK**2 9 - 4.0*DWPYHM*U**2*XGK - 4.0*DWPYGK*U**2*XGK A - 4.0*DWGKHM*T**2*XGK - 4.0*DGGK*T**2*XGK B - 4.0*DWILHM*S**2*XGK - 4.0*DWILGK*S**2*XGK C - S*U*WILPYORD + S*U*WILPYDIS - S*T*WILGKORD + S*T*WILGKDIS D - T*U*WGKPYORD + T*U*WGKPYDIS + 2.0*DWPYHM*U**2 E + 3.0*DGPY*U**2 + DWPYGK*T*U + DWGKPY*T*U + DWPYIL*S*U F + DWILPY*S*U + 2.0*DWGKHM*T**2 + 3.0*DGGK*T**2 + DWILGK*S*T G + DWGKIL*S*T + 2.0*DWILHM*S**2 + 3.0*DGIL*S**2 HEMATITE = (1.0-XIL-XGK-XPY)**2*EXP(RTLNHM/(8.3143*TK)) RTLNIL = WPYHM*XPY**2 + WPYHM*XIL*XPY - WILPYDIS*XIL*XPY 1 + WILHM*XIL*XPY + WPYHM*XGK*XPY - WGKPYDIS*XGK*XPY 2 + WGKHM*XGK*XPY - WPYHM*XPY + WILPYDIS*XPY - WILHM*XPY 3 - 4.0*DWPYHM*U**2*XPY - 4.0*DGPY*U**2*XPY 4 - 4.0*DWGKPY*T**2*XPY - 4.0*DWGKHM*T**2*XPY 5 - 4*DWILPY*S**2*XPY - 4.0*DWILHM*S**2*XPY 6 + 4.0*DWILPY*S*XPY + 4.0*DWILHM*S*XPY + WILHM*XIL**2 7 + WILHM*XGK*XIL - WILGKDIS*XGK*XIL + WGKHM*XGK*XIL 8 - 2.0*WILHM*XIL - 4.0*DWPYIL*U**2*XIL - 4.0*DWPYHM*U**2*XIL 9 - 4.0*DWGKIL*T**2*XIL - 4.0*DWGKHM*T**2*XIL A - 4.0*DWILHM*S**2*XIL - 4.0*DGIL*S**2*XIL B + 4.0*DWILHM*S*XIL + 4.0*DGIL*S*XIL + WGKHM*XGK**2 C - WILHM*XGK + WILGKDIS*XGK - WGKHM*XGK D - 4.0*DWPYHM*U**2*XGK - 4.0*DWPYGK*U**2*XGK E - 4.0*DWGKHM*T**2*XGK - 4.0*DGGK*T**2*XGK F - 4.0*DWILHM*S**2*XGK - 4.0*DWILGK*S**2*XGK G + 4.0*DWILHM*S*XGK + 4.0*DWILGK*S*XGK - S*U*WILPYORD H + U*WILPYORD + S*U*WILPYDIS - U*WILPYDIS + WILHM I - S*T*WILGKORD + T*WILGKORD + S*T*WILGKDIS -T*WILGKDIS J - T*U*WGKPYORD + T*U*WGKPYDIS + 2.0*DWPYIL*U**2 K + 4.0*DWPYHM*U**2 + 3.0*DGPY*U**2 + DWPYGK*T*U + DWGKPY*T*U L + DWPYIL*S*U + DWILPY*S*U - DWPYIL*U - DWILPY*U M + 2.0*DWGKIL*T**2 + 4.0*DWGKHM*T**2 + 3.0*DGGK*T**2 N + DWILGK*S*T + DWGKIL*S*T - DWILGK*T - DWGKIL*T O + 4.0*DWILHM*S**2 + 5.0*DGIL*S**2 - 4.0*DWILHM*S P - 6.0*DGIL*S + DGIL ILMENITE = 0.25*(XIL+S)*(XIL+S+XGK+T+XPY+U) 1 * EXP(RTLNIL/(8.3143*TK)) C GEIKIELITE = 0.0 PYROPHANITE = 0.0 C C Calculations for pure ilmenite C AUXVAR( 3) = 1.0 AUXVAR( 4) = 0.0 AUXVAR( 5) = 0.0 WORK(1) = 0.9 ! Initial guess for s WORK(2) = 0.0 WORK(3) = 0.0 ITMAX = 250 DEBUG = .FALSE. CALL VARMET(ORDER, 3, WORK, GSOLN, WORK(4), WORK(7), 1 WORK(10), WORK(13), WORK(16), 3, AUXVAR, ITMAX, DEBUG, 2 MODE) S_PURE = WORK(1) IF(MODE.NE.0) WRITE(*,*) 'Mode (VARMET) = ', MODE RTLNIL_PURE = DGIL*(1.0-S_PURE)**2 ILMENITE_PURE = 0.25*(1.0+S_PURE)**2 1 * EXP(RTLNIL_PURE/(8.3143*TK)) ILMENITE = ILMENITE/ILMENITE_PURE C LSTOP = .FALSE. RETURN END SUBROUTINE ORDER(B, P0, G, AUXVAR, LSTOP) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION AUXVAR(29), B(3), G(3) LOGICAL LSTOP C TK = AUXVAR( 1) P = AUXVAR( 2) XIL = AUXVAR( 3) XGK = AUXVAR( 4) XPY = AUXVAR( 5) HIL = AUXVAR( 6) SIL = AUXVAR( 7) HGK = AUXVAR( 8) SGK = AUXVAR( 9) HPY = AUXVAR(10) SPY = AUXVAR(11) DWILHM = AUXVAR(12) DWGKHM = AUXVAR(13) DWPYHM = AUXVAR(14) WILHM = AUXVAR(15) WGKHM = AUXVAR(16) WPYHM = AUXVAR(17) DWILGK = AUXVAR(18) DWGKIL = AUXVAR(19) DWILPY = AUXVAR(20) DWPYIL = AUXVAR(21) DWGKPY = AUXVAR(22) DWPYGK = AUXVAR(23) WILGKDIS = AUXVAR(24) WILGKORD = AUXVAR(25) WILPYDIS = AUXVAR(26) WILPYORD = AUXVAR(27) WGKPYDIS = AUXVAR(28) WGKPYORD = AUXVAR(29) C S = B(1) T = B(2) U = B(3) XHM = 1.0 - XIL - XGK - XPY IF((S.LT.0.0).OR.(S.GT.XIL).OR.(T.LT.0.0).OR.(T.GT.XGK).OR. 1 (U.LT.0.0).OR.(U.GT.XPY)) THEN LSTOP = .TRUE. RETURN END IF C IF((XHM.GT.0.0).AND.(XIL.GT.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.GT.0.0)) THEN DGIL = HIL - TK*SIL DGGK = HGK - TK*SGK DGPY = HPY - TK*SPY P0 = - WPYHM*XPY**2 + (-WPYHM+WILPYDIS-WILHM)*XIL*XPY 1 + (-WPYHM+WGKPYDIS-WGKHM)*XGK*XPY + (WPYHM+DGPY)*XPY 2 + (2.0*DWPYHM+2.0*DGPY)*U**2*XPY 3 + (2.0*DWGKPY+2.0*DWGKHM)*T**2*XPY 4 + (2.0*DWILPY+2.0*DWILHM)*S**2*XPY - WILHM*XIL**2 5 + (-WILHM+WILGKDIS-WGKHM)*XGK*XIL + (WILHM+DGIL)*XIL 6 + (2.0*DWPYIL+2.0*DWPYHM)*U**2*XIL 7 + (2.0*DWGKIL+2.0*DWGKHM)*T**2*XIL 8 + (2.0*DWILHM+2.0*DGIL)*S**2*XIL - WGKHM*XGK**2 9 + (WGKHM+DGGK)*XGK A + (2.0*DWPYHM+2.0*DWPYGK)*U**2*XGK B + (2.0*DWGKHM+2.0*DGGK)*T**2*XGK C + (2.0*DWILHM+2.0*DWILGK)*S**2*XGK D + S*U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) E + S*T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) F + T*U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) G + (-2.0*DWPYHM-3.0*DGPY)*U**2 H + (-2.0*DWGKHM-3.0*DGGK)*T**2 I + (-2.0*DWILHM-3.0*DGIL)*S**2 G(1) = 2.0*(2.0*DWILPY+2.0*DWILHM)*S*XPY 1 + 2.0*(2.0*DWILHM+2.0*DGIL)*S*XIL 2 + 2.0*(2.0*DWILHM+2.0*DWILGK)*S*XGK 3 + U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 4 + T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 5 + 2.0*(-2.0*DWILHM-3.0*DGIL)*S G(2) = 2.0*(2.0*DWGKPY+2.0*DWGKHM)*T*XPY 1 + 2.0*(2.0*DWGKIL+2.0*DWGKHM)*T*XIL 2 + 2.0*(2.0*DWGKHM+2.0*DGGK)*T*XGK 3 + S*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 4 + U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 5 + 2.0*(-2.0*DWGKHM-3.0*DGGK)*T G(3) = 2.0*(2.0*DWPYHM+2.0*DGPY)*U*XPY 1 + 2.0*(2.0*DWPYIL+2.0*DWPYHM)*U*XIL 2 + 2.0*(2.0*DWPYHM+2.0*DWPYGK)*U*XGK 3 + S*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 4 + T*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 5 + 2.0*(-2.0*DWPYHM-3.0*DGPY)*U C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XGK-T+XPY-U) 1 *LOG((XIL-S+XGK-T+XPY-U)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XGK+T+XPY+U) 1 *LOG((XIL+S+XGK+T+XPY+U)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.GT.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.GT.0.0)) THEN DGIL = HIL - TK*SIL DGGK = HGK - TK*SGK DGPY = HPY - TK*SPY P0 = WILPYDIS*XIL*XPY + WGKPYDIS*XGK*XPY + DGPY*XPY 1 + 2.0*DGPY*U**2*XPY + 2.0*DWGKPY*T**2*XPY 2 + 2.0*DWILPY*S**2*XPY + WILGKDIS*XGK*XIL + DGIL*XIL 3 + 2.0*DWPYIL*U**2*XIL + 2.0*DWGKIL*T**2*XIL 4 + 2.0*DGIL*S**2*XIL + DGGK*XGK + 2.0*DWPYGK*U**2*XGK 5 + 2.0*DGGK*T**2*XGK + 2.0*DWILGK*S**2*XGK 6 + S*U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 7 + S*T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 8 + T*U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 9 - 3.0*DGPY*U**2 - 3.0*DGGK*T**2 - 3.0*DGIL*S**2 G(1) = 4.0*DWILPY*S*XPY + 4.0*DGIL*S*XIL + 4.0*DWILGK*S*XGK 1 + U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 2 + T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 3 - 6.0*DGIL*S G(2) = 4.0*DWGKPY*T*XPY + 4.0*DWGKIL*T*XIL + 4.0*DGGK*T*XGK 1 + S*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 2 + U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 3 - 6.0*DGGK*T G(3) = 4.0*DGPY*U*XPY + 4.0*DWPYIL*U*XIL + 4.0*DWPYGK*U*XGK 1 + S*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 2 + T*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 4 - 6.0*DGPY*U C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XGK-T+XPY-U) 1 *LOG((XIL-S+XGK-T+XPY-U)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XGK+T+XPY+U) 1 *LOG((XIL+S+XGK+T+XPY+U)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T+XPY+U)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.LE.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.GT.0.0)) THEN DGGK = HGK - TK*SGK DGPY = HPY - TK*SPY P0 = - WPYHM*XPY**2 + (-WPYHM+WGKPYDIS-WGKHM)*XGK*XPY 1 + (WPYHM+DGPY)*XPY + (2.0*DWPYHM+2.0*DGPY)*U**2*XPY 2 + (2.0*DWGKPY+2.0*DWGKHM)*T**2*XPY - WGKHM*XGK**2 3 + (WGKHM+DGGK)*XGK + (2.0*DWPYHM+2.0*DWPYGK)*U**2*XGK 4 + (2.0*DWGKHM+2.0*DGGK)*T**2*XGK 5 + T*U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 6 + (-2.0*DWPYHM-3.0*DGPY)*U**2 7 + (-2.0*DWGKHM-3.0*DGGK)*T**2 G(1) = 0.0 G(2) = 2.0*(2.0*DWGKPY+2.0*DWGKHM)*T*XPY 1 + 2.0*(2.0*DWGKHM+2.0*DGGK)*T*XGK 2 + U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 3 + 2.0*(-2.0*DWGKHM-3.0*DGGK)*T G(3) = 2.0*(2.0*DWPYHM+2.0*DGPY)*U*XPY 1 + 2.0*(2.0*DWPYHM+2.0*DWPYGK)*U*XGK 2 + T*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 3 + 2.0*(-2.0*DWPYHM-3.0*DGPY)*U C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK-T+XPY-U)*LOG((XGK-T+XPY-U)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XGK-T+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK+T+XPY+U)*LOG((XGK+T+XPY+U)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XGK+T+XPY+U)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.GT.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.GT.0.0)) THEN DGIL = HIL - TK*SIL DGPY = HPY - TK*SPY P0 = - WPYHM*XPY**2 + (-WPYHM+WILPYDIS-WILHM)*XIL*XPY 1 + (WPYHM+DGPY)*XPY + (2.0*DWPYHM+2.0*DGPY)*U**2*XPY 2 + (2.0*DWILPY+2.0*DWILHM)*S**2*XPY - WILHM*XIL**2 3 + (WILHM+DGIL)*XIL + (2.0*DWPYIL+2.0*DWPYHM)*U**2*XIL 4 + (2.0*DWILHM+2.0*DGIL)*S**2*XIL 6 + S*U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 7 + (-2.0*DWPYHM-3.0*DGPY)*U**2 8 + (-2.0*DWILHM-3.0*DGIL)*S**2 G(1) = 2.0*(2.0*DWILPY+2.0*DWILHM)*S*XPY 1 + 2.0*(2.0*DWILHM+2.0*DGIL)*S*XIL 2 + U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 3 + 2.0*(-2.0*DWILHM-3.0*DGIL)*S G(2) = 0.0 G(3) = 2.0*(2.0*DWPYHM+2.0*DGPY)*U*XPY 1 + 2.0*(2.0*DWPYIL+2.0*DWPYHM)*U*XIL 2 + S*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 3 + 2.0*(-2.0*DWPYHM-3.0*DGPY)*U C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XPY-U)*LOG((XIL-S+XPY-U)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XIL-S+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XPY+U)*LOG((XIL+S+XPY+U)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XIL+S+XPY+U)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.GT.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.LE.0.0)) THEN DGIL = HIL - TK*SIL DGGK = HGK - TK*SGK P0 = - WILHM*XIL**2 + (-WILHM+WILGKDIS-WGKHM)*XGK*XIL 1 + (WILHM+DGIL)*XIL + (2.0*DWGKIL+2.0*DWGKHM)*T**2*XIL 2 + (2.0*DWILHM+2.0*DGIL)*S**2*XIL - WGKHM*XGK**2 3 + (WGKHM+DGGK)*XGK + (2.0*DWGKHM+2.0*DGGK)*T**2*XGK 4 + (2.0*DWILHM+2.0*DWILGK)*S**2*XGK 5 + S*T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 6 + (-2.0*DWGKHM-3.0*DGGK)*T**2 7 + (-2.0*DWILHM-3.0*DGIL)*S**2 G(1) = 2.0*(2.0*DWILHM+2.0*DGIL)*S*XIL 1 + 2.0*(2.0*DWILHM+2.0*DWILGK)*S*XGK 2 + T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 3 + 2.0*(-2.0*DWILHM-3.0*DGIL)*S G(2) = 2.0*(2.0*DWGKIL+2.0*DWGKHM)*T*XIL 1 + 2.0*(2.0*DWGKHM+2.0*DGGK)*T*XGK 2 + S*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) 3 + 2.0*(-2.0*DWGKHM-3.0*DGGK)*T G(3) = 0.0 C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XGK-T)*LOG((XIL-S+XGK-T)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XGK+T)*LOG((XIL+S+XGK+T)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.LE.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.GT.0.0)) THEN DGGK = HGK - TK*SGK DGPY = HPY - TK*SPY P0 = WGKPYDIS*XGK*XPY + DGPY*XPY + 2.0*DGPY*U**2*XPY 1 + 2.0*DWGKPY*T**2*XPY + DGGK*XGK + 2.0*DWPYGK*U**2*XGK 2 + 2.0*DGGK*T**2*XGK + T*U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) 3 - 3.0*DGPY*U**2 - 3.0*DGGK*T**2 G(1) = 0.0 G(2) = 4.0*DWGKPY*T*XPY + 4.0*DGGK*T*XGK - 6.0*DGGK*T 1 + U*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) G(3) = 4.0*DGPY*U*XPY + 4.0*DWPYGK*U*XGK - 6.0*DGPY*U 1 + T*(WGKPYORD-WGKPYDIS-DWPYGK-DWGKPY) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK-T+XPY-U)*LOG((XGK-T+XPY-U)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XGK-T+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK+T+XPY+U)*LOG((XGK+T+XPY+U)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XGK+T+XPY+U)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.GT.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.GT.0.0)) THEN DGIL = HIL - TK*SIL DGPY = HPY - TK*SPY P0 = WILPYDIS*XIL*XPY + DGPY*XPY + 2.0*DGPY*U**2*XPY 1 + 2.0*DWILPY*S**2*XPY + DGIL*XIL + 2.0*DWPYIL*U**2*XIL 2 + 2.0*DGIL*S**2*XIL + S*U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) 3 - 3.0*DGPY*U**2 - 3.0*DGIL*S**2 G(1) = 4.0*DWILPY*S*XPY + 4.0*DGIL*S*XIL - 6.0*DGIL*S 1 + U*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) G(2) = 0.0 G(3) = 4.0*DGPY*U*XPY + 4.0*DWPYIL*U*XIL - 6.0*DGPY*U 1 + S*(WILPYORD-WILPYDIS-DWPYIL-DWILPY) C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XPY-U)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XPY-U)*LOG((XIL-S+XPY-U)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*0.5*LOG((XIL-S+XPY-U)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XPY+U)*LOG((XIL+S+XPY+U)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*0.5*LOG((XIL+S+XPY+U)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.GT.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.LE.0.0)) THEN DGIL = HIL - TK*SIL DGGK = HGK - TK*SGK P0 = WILGKDIS*XGK*XIL + DGIL*XIL + 2.0*DWGKIL*T**2*XIL 1 + 2.0*DGIL*S**2*XIL + DGGK*XGK + 2.0*DGGK*T**2*XGK 2 + 2.0*DWILGK*S**2*XGK - 3.0*DGGK*T**2 - 3.0*DGIL*S**2 3 + S*T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) G(1) = 4.0*DGIL*S*XIL + 4.0*DWILGK*S*XGK - 6.0*DGIL*S 1 + T*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) G(2) = 4.0*DWGKIL*T*XIL + 4.0*DGGK*T*XGK - 6.0*DGGK*T 1 + S*(WILGKORD-WILGKDIS-DWILGK-DWGKIL) G(3) = 0.0 C P0 = P0 + 8.3143*TK*0.5*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XGK-T)/2.0) C P0 = P0 + 8.3143*TK*0.5*(XIL-S+XGK-T)*LOG((XIL-S+XGK-T)/2.0) G(1) = G(1) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*0.5*LOG((XIL-S+XGK-T)/2.0) P0 = P0 + 8.3143*TK*0.5*(XIL+S+XGK+T)*LOG((XIL+S+XGK+T)/2.0) G(1) = G(1) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*0.5*LOG((XIL+S+XGK+T)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.GT.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.LE.0.0)) THEN DGIL = HIL - TK*SIL P0 = - WILHM*XIL**2 + (WILHM+DGIL)*XIL 1 + (2.0*DWILHM+2.0*DGIL)*S**2*XIL 2 + (-2.0*DWILHM-3.0*DGIL)*S**2 G(1) = 2.0*(2.0*DWILHM+2.0*DGIL)*S*XIL 1 + 2.0*(-2.0*DWILHM-3.0*DGIL)*S G(2) = 0.0 G(3) = 0.0 C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*(XIL+S)*LOG((XIL+S)/2.0) G(1) = G(1) + 8.3143*TK*LOG((XIL+S)/2.0) P0 = P0 + 8.3143*TK*(XIL-S)*LOG((XIL-S)/2.0) G(1) = G(1) - 8.3143*TK*LOG((XIL-S)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.LE.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.LE.0.0)) THEN DGGK = HGK - TK*SGK P0 = - WGKHM*XGK**2 + (WGKHM+DGGK)*XGK 1 + (2.0*DWGKHM+2.0*DGGK)*T**2*XGK 2 + (-2.0*DWGKHM-3.0*DGGK)*T**2 G(1) = 0.0 G(2) = 2.0*(2.0*DWGKHM+2.0*DGGK)*T*XGK 1 + 2.0*(-2.0*DWGKHM-3.0*DGGK)*T G(3) = 0.0 C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*(XGK+T)*LOG((XGK+T)/2.0) G(2) = G(2) + 8.3143*TK*LOG((XGK+T)/2.0) P0 = P0 + 8.3143*TK*(XGK-T)*LOG((XGK-T)/2.0) G(2) = G(2) - 8.3143*TK*LOG((XGK-T)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.LE.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.GT.0.0)) THEN DGPY = HPY - TK*SPY P0 = - WPYHM*XPY**2 + (WPYHM+DGPY)*XPY 1 + (2.0*DWPYHM+2.0*DGPY)*U**2*XPY 2 + (-2.0*DWPYHM-3.0*DGPY)*U**2 G(1) = 0.0 G(2) = 0.0 G(3) = 2.0*(2.0*DWPYHM+2.0*DGPY)*U*XPY 1 + 2.0*(-2.0*DWPYHM-3.0*DGPY)*U C P0 = P0 + 8.3143*TK*2.0*XHM*LOG(XHM) C P0 = P0 + 8.3143*TK*(XPY+U)*LOG((XPY+U)/2.0) G(3) = G(3) + 8.3143*TK*LOG((XPY+U)/2.0) P0 = P0 + 8.3143*TK*(XPY-U)*LOG((XPY-U)/2.0) G(3) = G(3) - 8.3143*TK*LOG((XPY-U)/2.0) C ELSE IF((XHM.GT.0.0).AND.(XIL.LE.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.LE.0.0)) THEN P0 = 0.0 G(1) = 0.0 G(2) = 0.0 G(3) = 0.0 C ELSE IF((XHM.LE.0.0).AND.(XIL.GT.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.LE.0.0)) THEN DGIL = HIL - TK*SIL P0 = DGIL*(1.0-S**2) G(1) = - 2.0*DGIL*S G(2) = 0.0 G(3) = 0.0 C P0 = P0 + 8.3143*TK*(1.0+S)*LOG((1.0+S)/2.0) G(1) = G(1) + 8.3143*TK*LOG((1.0+S)/2.0) P0 = P0 + 8.3143*TK*(1.0-S)*LOG((1.0-S)/2.0) G(1) = G(1) - 8.3143*TK*LOG((1.0-S)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.LE.0.0).AND.(XGK.GT.0.0).AND. 1 (XPY.LE.0.0)) THEN DGGK = HGK - TK*SGK P0 = DGGK*(1.0-T**2) G(1) = 0.0 G(2) = - 2.0*DGGK*T G(3) = 0.0 C P0 = P0 + 8.3143*TK*(1.0+T)*LOG((1.0+T)/2.0) G(2) = G(2) + 8.3143*TK*LOG((1.0+T)/2.0) P0 = P0 + 8.3143*TK*(1.0-T)*LOG((1.0-T)/2.0) G(2) = G(2) - 8.3143*TK*LOG((1.0-T)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.LE.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.GT.0.0)) THEN DGPY = HPY - TK*SPY P0 = DGPY*(1.0-U**2) G(1) = 0.0 G(2) = 0.0 G(3) = - 2.0*DGPY*U C P0 = P0 + 8.3143*TK*(1.0+U)*LOG((1.0+U)/2.0) G(3) = G(3) + 8.3143*TK*LOG((1.0+U)/2.0) P0 = P0 + 8.3143*TK*(1.0-U)*LOG((1.0-U)/2.0) G(3) = G(3) - 8.3143*TK*LOG((1.0-U)/2.0) C ELSE IF((XHM.LE.0.0).AND.(XIL.LE.0.0).AND.(XGK.LE.0.0).AND. 1 (XPY.LE.0.0)) THEN STOP 'Absurd logic error in call to min(s,t,u) function.' C END IF C P0 = P0/1000.0 G(1) = G(1)/1000.0 G(2) = G(2)/1000.0 G(3) = G(3)/1000.0 C LSTOP = .FALSE. RETURN END SUBROUTINE VARMET(SUB, N, B, P0, G, X, C, T, BMAT, NB, AUXVAR, 1 ITMAX, DEBUG, MODE) C C Algorithm 21. Variable metric method. (page 159-160) C Nash, J.C. Compact Numerical Methods for Computers: Linear C Algebra and Function Minimisation. John Wiley and Sons, C New York, 227 pages (1979) C C SUB(B, P0, G, AUXVAR, LSTOP) SUBROUTINE which returns the value C of the function (P0) and its gradient (G) at the "point" B C subject to values of auxillary "known" parameters contained C in the array AUXVAR. LSTOP is returned with a value of .TRUE. C if the function or its gradient cannot be evaluated at B. C N Number of parameters in the function C B Parameter values C ON INPUT: Initial guess C ON OUTPUT: Final values of the function parameters which C minimize the function defined in S C P0 Function value at the minimum C G Gradient of the function evaluated at the minimum C X Working array of length N C C Working array of length N C T Working array of length N C BMAT Woorking square matrix of size N*N C NB Row dimension of BMAT C AUXVAR User supplied auxillary variables passed to S() in order C to evaluate the function at its gradient C ITMAX Maximum number of allowed calls to S() C DEBUG Logical variable, set to .TRUE. in DEBUG mode C MODE Status variable, On exit Mode equals C 0 when computation is successful C 1 when initial point is infeasible C 2 when iteration limit is exceeded C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION AUXVAR(*), B(N), G(N), X(N), C(N), T(N), BMAT(NB, N) INTEGER COUNT DOUBLE PRECISION K LOGICAL LSTOP, DEBUG, CONVERGED, RESTART C C Step 0 Enter n, the number of parameters in the function C Enter b[i], i=1,2,...,n, the starting values for the C parameters C Let ig = 0 (graient evaluation counter) C Let ifn = 0 (function evaluation counter) C Let w = 0.2, let tolerance = 0.0001 to define acceptable C point strategy C IG = 0 IFN = 0 W = 0.2 TOLERANCE = 0.0001 MODE = 0 CONVERGED = .FALSE. RESTART = .TRUE. C C Step 1 Compute P0 = S(b). If this is not possible, stop, as the C initial point is infeasible. Let ifn = ifn + 1 C Step 2 Compute g(b), the gradient at the current point. C Let ig = ig + 1 C CALL SUB(B, P0, G, AUXVAR, LSTOP) IF(LSTOP) THEN MODE = 1 RETURN END IF IFN = IFN + 1 IG = IG + 1 DO WHILE(.NOT.CONVERGED) C C Step 3 Set BMAT to the unit matrix 1n. C Let ilast = ig to mark point at which B most recently C set to unity C IF(RESTART) THEN DO I=1,N DO J=1,N BMAT(I,J) = 0.0 END DO BMAT(I,I) = 1.0 END DO ILAST = IG RESTART = .FALSE. END IF C C The top of the main iteration loop. To follow the progress of the C minimisation ig, ifn, and P0 could be printed here. C IF(DEBUG) THEN WRITE(*,'(1X,''P0, IG, IFN:'', G12.6,2I10)') P0, IG, IFN END IF C DO I=1,N X(I) = B(I) C(I) = G(I) END DO C C Step 5 Form t = -Bg and t^T g C D1 = 0.0 DO I=1,N S = 0.0 DO J=1,N S = S - BMAT(I,J)*G(J) END DO T(I) = S D1 = D1 - S*G(I) END DO C C Step 6 If D1>0, goto step 7 to perform linear search. Otherwise, C search direction is uphill, so if ilast=ig, goto step 18, C since the algorithm is taken to have converged. Otherwise, C goto step 3 to restart iteration. C IF(D1.GT.0.0) THEN C C Step 7 Let k=1, to intialise step length C K = 1.0 C C Step 8 Let count = 0 C LSTOP = .TRUE. DO WHILE(LSTOP) COUNT = 0 DO I=1,N B(I) = X(I) + K*T(I) IF(B(I).EQ.X(I)) COUNT = COUNT + 1 END DO C C Step 9 Convergence test C If Count < n, goto step 10 to evaluate the function C Otherwise, if ilast=ig, goto step 18, since the C algorithm is taken to have converged. Otherwise goto C step 3 to restart iteration C IF(COUNT.LT.N) THEN C C Step 10 Compute P=S(b), let ifn=ifn+1. If the function cannot be C computed, let k = w*k, then goto step 8 to repeat the C adjustment of parameters C CALL SUB(B, P, G, AUXVAR, LSTOP) IFN = IFN + 1 IF(IFN.GT.ITMAX) THEN MODE = 2 RETURN END IF IF(LSTOP) THEN K = W*K ELSE C C Step 11 Acceptance test on new point C If P < P0 - D1*k*tolerance, goto 12 to continue iteration. C Otherwise, let k=w*k, then goto step 8 C IF(P.LT.(P0-D1*K*TOLERANCE)) THEN C C Step 12 Let p0=P to save the function value. C Compute g at b, the gradient at the new point. C Let ig = ig + 1 C P0 = P IG = IG + 1 C C Step 13 Prepare for matrix update. C D1 = 0.0 DO I=1,N T(I) = K*T(I) C(I) = G(I) - C(I) D1 = D1 + T(I)*C(I) END DO C C Step 14 If D1 <= 0, goto step 3 to restart iteration since positive C definiteness of matrix cannot be guaranteed and there is C a danger of a division by zero C IF(D1.LE.0.0) THEN RESTART = .TRUE. ELSE C C Step 15 Computation of BMATy, y^T BMAT y C D2 = 0.0 DO I=1,N S = 0.0 DO J=1,N S = S + BMAT(I,J)*C(J) END DO X(I) = S D2 = D2 + S*C(I) END DO C C Step 16 Let D2 = 1 + D2/D1. Note this is (15.36) times d1. C D2 = 1.0 + D2/D1 DO I=1,N DO J=1,N BMAT(I,J) = BMAT(I,J) - (T(I)*X(J) 1 + X(I)*T(J) 2 - D2*T(I)*T(J)) / D1 END DO END DO C C This completes the update. C C Note that there are other ways to organise the update. The C rounding errors involved in some cases may cause different C schemes to have very different convergence patterns. So far no C one scheme seems to be better than the others, so the most direct C implementation has been chosen here. C C Step 17 Goto step 4 for another iteration. C END IF ELSE K = W*K LSTOP = .TRUE. END IF END IF ELSE IF(IG.EQ.ILAST) THEN CONVERGED = .TRUE. LSTOP = .FALSE. ELSE RESTART = .TRUE. LSTOP = .FALSE. END IF END DO ELSE IF(IG.EQ.ILAST) THEN CONVERGED = .TRUE. ELSE RESTART = .TRUE. END IF END DO C C Step 18 PRINT IG, IFN, P0, b[i], i=1,2,...,n. C STOP. C IF(DEBUG) THEN WRITE(*,'(1X,''IG, IFN, P0:'',2I10,G12.6)') IG, IFN, P0 END IF C RETURN END SUBROUTINE ACT_SPN(T, X2, X3, X4, X5, S1, S2, S3, S4, 1 XMG2TET, XFE2TET, XAL3TET, XFE3TET, XCR3TET, XMG2OCT, 2 XFE2OCT, XAL3OCT, XFE3OCT, XTI4OCT, XCR3OCT, 3 HERCYNITE, SPINEL, ULVOSPINEL, MAGNETITE, MGTITANATE, 4 MGFERRITE, CHROMITE, MGCHROMITE, LSTOP) C C T (input) Temperature (K) C X2 (input) mole fraction MgAl2O4 C X3 (input) mole fraction FeCr2O3 C X4 (input) mole fraction Fe2TiO4 C X5 (input) mole fraction Fe3O4 C S1 (output) ordering variable C S2 (output) ordering variable C S3 (output) ordering variable C S4 (output) ordering variable C XMG2TET (output) X Mg2+ on tetrahedral site C XFE2TET (output) X Fe2+ on tetrahedral site C XAL3TET (output) X Al3+ on tetrahedral site C XFE3TET (output) X Fe3+ on tetrahedral site C XCR3TET (output) X Cr3+ on tetrahedral site C XMG2OCT (output) X Mg2+ on octahedral site C XFE2OCT (output) X Fe2+ on octahedral site C XAL3OCT (output) X Al3+ on octahedral site C XFE3OCT (output) X Fe3+ on octahedral site C XTI4OCT (output) X Ti4+ on octahedral site C XCR3OCT (output) X Cr3+ on octahedral site C HERCYNITE (output) a of hercynite [ ref to Fe2+(Al,Al)2O4 ] C SPINEL (output) a of spinel [ ref to Mg(Al,Al)2O4 ] C ULVOSPINEL (output) a of ulvospinel [ ref to Fe2+(Fe2+,Ti)2O4 ] C MAGNETITE (output) a of magnetite [ ref to Fe3+(Fe2+,Fe3+)2O4 ] C MGTITANATE (output) a of Mg-titanate [ ref to Mg(Mg,Ti)2O4 ] C MGFERRITE (output) a of Mg-ferrite [ ref to Fe3+(Mg,Fe3+)2O4 ] C CHROMITE (output) a of chromite [ ref to Fe2+(Cr3+,Cr3+)2O4 ] C MGCHROMITE (output) a of Mg-chromite [ ref to Mg2+(Cr3+,Cr3+)2O4 ] C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C DOUBLE PRECISION MAGNETITE, MGTITANATE, MGFERRITE, MGCHROMITE LOGICAL LSTOP 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 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.40 S21 = -0.0098 W22 = 3.6 H24 = 6.55 S24 = 0.0 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 = 12.5 W3P4 = 10.0 W3PU4U = 10.4 W35 = 0.0 W3P5 = 10.0 W35P = 8.0 W3P5P = 0.0 C LSTOP = .FALSE. 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(*,1000) 1000 FORMAT(1X,78('*'),/, 1 1X,'%SPN_F_NCONV, Failed to converge in ordering parameter ', 2 'subroutine.',/, 2 1X,78('*')) HERCYNITE = 0.0 SPINEL = 0.0 ULVOSPINEL = 0.0 MAGNETITE = 0.0 MGTITANATE = 0.0 MGFERRITE = 0.0 CHROMITE = 0.0 MGCHROMITE = 0.0 RETURN END IF IF( ((XMG2TET.LT.0.0).AND.(ABS(XMG2TET).GT.1.0D-12)).OR. 1 ((XFE2TET.LT.0.0).AND.(ABS(XFE2TET).GT.1.0D-12)).OR. 2 ((XAL3TET.LT.0.0).AND.(ABS(XAL3TET).GT.1.0D-12)).OR. 3 ((XFE3TET.LT.0.0).AND.(ABS(XFE3TET).GT.1.0D-12)).OR. 4 ((XCR3TET.LT.0.0).AND.(ABS(XCR3TET).GT.1.0D-12)).OR. 5 ((XMG2OCT.LT.0.0).AND.(ABS(XMG2OCT).GT.1.0D-12)).OR. 6 ((XFE2OCT.LT.0.0).AND.(ABS(XFE2OCT).GT.1.0D-12)).OR. 7 ((XAL3OCT.LT.0.0).AND.(ABS(XAL3OCT).GT.1.0D-12)).OR. 8 ((XFE3OCT.LT.0.0).AND.(ABS(XFE3OCT).GT.1.0D-12)).OR. 9 ((XTI4OCT.LT.0.0).AND.(ABS(XTI4OCT).GT.1.0D-12)).OR. A ((XCR3OCT.LT.0.0).AND.(ABS(XCR3OCT).GT.1.0D-12)) ) THEN LSTOP = .TRUE. WRITE(*,1001) 1001 FORMAT(1X,78('*'),/, 1 1X,'%SPN_F_DEPN, Failed to converge. Ordering parameter ', 2 'relationship detected.',/, 3 1X,78('*')) HERCYNITE = 0.0 SPINEL = 0.0 ULVOSPINEL = 0.0 MAGNETITE = 0.0 MGTITANATE = 0.0 MGFERRITE = 0.0 CHROMITE = 0.0 MGCHROMITE = 0.0 IF(XMG2TET.LT.0.0) WRITE(*,1002) 'XMg2+ tet', XMG2TET 1002 FORMAT(14X,A9,' is equal to ',G20.13) IF(XFE2TET.LT.0.0) WRITE(*,1002) 'XFe2+ tet', XFE2TET IF(XAL3TET.LT.0.0) WRITE(*,1002) 'XAl3+ tet', XAL3TET IF(XFE3TET.LT.0.0) WRITE(*,1002) 'XFe3+ tet', XFE3TET IF(XCR3TET.LT.0.0) WRITE(*,1002) 'XCr3+ tet', XCR3TET IF(XMG2OCT.LT.0.0) WRITE(*,1002) 'XMg2+ oct', XMG2OCT IF(XFE2OCT.LT.0.0) WRITE(*,1002) 'XFe2+ oct', XFE2OCT IF(XAL3OCT.LT.0.0) WRITE(*,1002) 'XAl3+ oct', XAL3OCT IF(XFE3OCT.LT.0.0) WRITE(*,1002) 'XFe3+ oct', XFE3OCT IF(XTI4OCT.LT.0.0) WRITE(*,1002) 'XTi4+ oct', XTI4OCT IF(XCR3OCT.LT.0.0) WRITE(*,1002) 'XCr3+ oct', XCR3OCT RETURN END IF C G1 = 0.0 ! Standard state parameters in MACSYMA expressions G2 = 0.0 ! are zeroed out (this leaves them intact G3 = 0.0 G4 = 0.0 G5 = 0.0 R = 0.00198726 C G11 = H11 - T*S11 G33 = H33 - T*S33 G55 = H55 - T*S55 ! Not used in expressions below GX = HX - T*SX GEX = HEX - T*SEX G24 = H24 - T*S24 G23 = H23 - T*S23 G25 = H25 - T*S25 C C Activity expression generated by MACSYMA C HERCYNITE = 0.0 SPINEL = 0.0 C ULVOSPINEL = W15*X5**2-W45*X4*X5+W15*X4*X5+W14*X4*X5-W35*X3*X5+W15 1 *X3*X5+W13*X3*X5-W2P5U*X2*X5/2.0+W22*X2*X5/2.0+W1P5*X2*X5/2.0-W 2 11*X2*X5/2.0-G25*X2*X5-S4*W55*X5+W45*X5-S3*W3P5*X5+S3*W35*X5-S1 3 *W2P5U*X5/2.0+S1*W22*X5/2.0-S2*W1P5*X5+S1*W1P5*X5/2.0+S4*W15P*X 4 5-S4*W15*X5+S2*W15*X5-W15*X5-W14*X5+S3*W13P*X5-S3*W13*X5+S2*W11 5 *X5-S1*W11*X5/2.0+W14*X4**2-W34*X3*X4+W14*X3*X4+W13*X3*X4-WTET* 6 X2*X4/2.0+W24U*X2*X4/2.0-W14*X2*X4/2.0-GX*X2*X4/4.0+GEX*X2*X4/4 7 .0-G24*X2*X4/2.0-S1*WTET*X4/2.0-S4*W45P*X4+S4*W45*X4-S3*W3P4*X4 8 +S3*W34*X4+S1*W24U*X4/2.0-S2*W1P4*X4+S4*W15P*X4-S4*W15*X4+S2*W1 9 4*X4-S1*W14*X4/2.0-2*W14*X4+S3*W13P*X4 ULVOSPINEL = ULVOSPINEL -S3*W13*X4+S2*W11*X4+GX*S1 : *X4/4.0-GEX*S1*X4/4.0+G24*S1*X4/2.0+W13*X3**2-W2P3U*X2*X3/2.0+ ; W22*X2*X3/2.0+W1P3*X2*X3/2.0-W11*X2*X3/2.0-G23*X2*X3-S4*W35P*X3 < +S4*W35*X3+W34*X3-S3*W33*X3-S1*W2P3U*X3/2.0+S1*W22*X3/2.0-S2*W1 = P3*X3+S1*W1P3*X3/2.0+S4*W15P*X3-S4*W15*X3-W14*X3+S3*W13P*X3-S3* > W13*X3+S2*W13*X3-W13*X3+S2*W11*X3-S1*W11*X3/2.0+WTET*X2**2/4.0+ ? WOCT*X2**2/4.0+GX*X2**2/4.0-S4*WTET*X2/2.0-S3*WTET*X2/2.0-S2*WT @ ET*X2/2.0+S1*WTET*X2/2.0+WTET*X2/2.0+S4*WOCT*X2/2.0+S3*WOCT*X2/ 1 2.0+S2*WOCT*X2/2.0-S1*WOCT*X2/2.0-S4*W4U5PU*X2+S4*W45P*X2-S3*W3 2 PU4U*X2+S3*W3P4*X2+S4*W2P5U*X2/2.0-S2*W2P4U*X2+S3*W2P3U*X2/2.0+ 3 S4*W25PU*X2/2.0-W24U*X2/2.0 ULVOSPINEL = ULVOSPINEL +S3*W23PU*X2/2.0-S4*W22*X2/2.0-S3*W22 4 *X2/2.0+S2*W22*X2/2.0-S4*W1P5*X2/2.0+S2*W1P4*X2-S3*W1P3*X2/2.0 5 -S4*W15P*X2/2.0+W14*X2/2.0-S3*W13P*X2/2.0+S4*W11*X2/2.0+S3*W11* 6 X2/2.0-S2*W11*X2/2.0+GEX*S4*X2/2.0+G25*S4*X2-G24*S4*X2+GEX*S3*X 7 2/2.0-G24*S3*X2+G23*S3*X2+GEX*S2*X2/2.0-G24*S2*X2+GX*X2/4.0-GEX 8 *X2/4.0+G24*X2/2.0-S1*S4*WTET/2.0-S1*S3*WTET/2.0-S1*S2*WTET/2.0 9 +S1**2*WTET/4.0+S1*WTET/2.0-S1*S4*WOCT/2.0-S1*S3*WOCT/2.0-S1*S2 : *WOCT/2.0+S1**2*WOCT/4.0+S4**2*W55 ULVOSPINEL = ULVOSPINEL +S4*W45P-S4*W45-S3*S4*W3P5P+S3 ; *S4*W3P5+S3*W3P4+S3*S4*W35P-S3*S4*W35-S3*W34+S3**2*W33+S1*S4*W < 2P5U/2.0+S1*S3*W2P3U/2.0+S1*S4*W25PU/2.0-S1*W24U/2.0+S1*S3*W23P = U/2.0-S1*S4*W22/2.0-S1*S3*W22/2.0+S1*S2*W22/2.0-S2*S4*W1P5P+S2* > S4*W1P5-S1*S4*W1P5/2.0+S2*W1P4-S2*S3*W1P3P+S2*S3*W1P3-S1*S3*W1P ? 3/2.0+S2*S4*W15P-S1*S4*W15P/2.0-S4*W15P-S2*S4*W15+S4*W15-S2*W14 @ +S1*W14/2.0+W14+S2*S3*W13P-S1*S3*W13P/2.0-S3*W13P-S2*S3*W13+S3* 1 W13+S1*S4*W11/2.0+S1*S3*W11/2.0+S2**2*W11-S1*S2*W11/2.0-S2*W11+ 2 GX*S1*S4/2.0+GX*S1*S3/2.0+GX*S1*S2/2.0-GX*S1**2/4.0-GX*S1/4.0+G 3 EX*S1/4.0-G24*S1/2.0+G4 C IF((XFE2TET.GT.0.0).AND. 1 (XFE2OCT.GT.0.0).AND.(XTI4OCT.GT.0.0)) THEN ULVOSPINEL = ULVOSPINEL + R*T*(LOG(XFE2TET) + LOG(XFE2OCT) 1 + LOG(XTI4OCT)) C Pure Ulvospinel: [x2=0, x3=0, x4=1, x5=0, s1=0, s2=0, s3=0, s4=0] PURE = - 2.0*R*T*LOG(2.0) + G4 ULVOSPINEL = EXP((ULVOSPINEL-PURE)/(R*T)) ELSE ULVOSPINEL = 0.0 END IF C MAGNETITE = W15*X5**2-W45*X4*X5+W15*X4*X5+W14*X4*X5-W35*X3*X5+W15* 1 X3*X5+W13*X3*X5-W2P5U*X2*X5/2.0+W22*X2*X5/2.0+W1P5*X2*X5/2.0-W1 2 1*X2*X5/2.0-G25*X2*X5-S4*W55*X5-S3*W3P5*X5+S3*W35*X5-S1*W2P5U*X 3 5/2.0+S1*W22*X5/2.0-S2*W1P5*X5+S1*W1P5*X5/2.0+S4*W15P*X5-S4*W15 4 *X5+S2*W15*X5-2*W15*X5+S3*W13P*X5-S3*W13*X5+S2*W11*X5-S1*W11*X5 5 /2.0+W14*X4**2-W34*X3*X4+W14*X3*X4+W13*X3*X4-WTET*X2*X4/2.0+W24 6 U*X2*X4/2.0-W14*X2*X4/2.0-GX*X2*X4/4.0+GEX*X2*X4/4.0-G24*X2*X4/ 7 2.0-S1*WTET*X4/2.0-S4*W45P*X4+S4*W45*X4+W45*X4-S3*W3P4*X4+S3*W3 8 4*X4+S1*W24U*X4/2.0-S2*W1P4*X4+S4*W15P*X4-S4*W15*X4-W15*X4+S2*W 9 14*X4-S1*W14*X4/2.0-W14*X4+S3*W13P*X4 MAGNETITE = MAGNETITE -S3*W13*X4+S2*W11*X4+GX*S1 : *X4/4.0-GEX*S1*X4/4.0+G24*S1*X4/2.0+W13*X3**2-W2P3U*X2*X3/2.0+W ; 22*X2*X3/2.0+W1P3*X2*X3/2.0-W11*X2*X3/2.0-G23*X2*X3-S4*W35P*X3+ < S4*W35*X3+W35*X3-S3*W33*X3-S1*W2P3U*X3/2.0+S1*W22*X3/2.0-S2*W1P = 3*X3+S1*W1P3*X3/2.0+S4*W15P*X3-S4*W15*X3-W15*X3+S3*W13P*X3-S3*W > 13*X3+S2*W13*X3-W13*X3+S2*W11*X3-S1*W11*X3/2.0+WTET*X2**2/4.0+W ? OCT*X2**2/4.0+GX*X2**2/4.0-S4*WTET*X2/2.0-S3*WTET*X2/2.0-S2*WTE @ T*X2/2.0+S1*WTET*X2/2.0+S4*WOCT*X2/2.0+S3*WOCT*X2/2.0+S2*WOCT*X 1 2/2.0-S1*WOCT*X2/2.0-S4*W4U5PU*X2+S4*W45P*X2-S3*W3PU4U*X2+S3*W3 2 P4*X2+S4*W2P5U*X2/2.0+W2P5U*X2/2.0-S2*W2P4U*X2+S3*W2P3U*X2/2.0+ 3 S4*W25PU*X2/2.0+S3*W23PU*X2/2.0 MAGNETITE = MAGNETITE -S4*W22*X2/2.0-S3*W22*X2/2.0+S2* 4 W22*X2/2.0-W22*X2/2.0-S4*W1P5*X2/2.0-W1P5*X2/2.0+S2*W1P4*X2-S3* 5 W1P3*X2/2.0-S4*W15P*X2/2.0-S3*W13P*X2/2.0+S4*W11*X2/2.0+S3*W11* 6 X2/2.0-S2*W11*X2/2.0+W11*X2/2.0+GEX*S4*X2/2.0+G25*S4*X2-G24*S4* 7 X2+GEX*S3*X2/2.0-G24*S3*X2+G23*S3*X2+GEX*S2*X2/2.0-G24*S2*X2+G2 8 5*X2-S1*S4*WTET/2.0-S1*S3*WTET/2.0-S1*S2*WTET/2.0+S1**2*WTET/4. 9 0-S1*S4*WOCT/2.0-S1*S3*WOCT/2.0-S1*S2*WOCT/2.0+S1**2*WOCT/4.0+S : 4**2*W55+S4*W55-S3*S4*W3P5P+S3*S4*W3P5+S3*W3P5 MAGNETITE = MAGNETITE +S3*S4*W35P-S3*S4 ; *W35-S3*W35+S3**2*W33+S1*S4*W2P5U/2.0+S1*W2P5U/2.0+S1*S3*W2P3U/ < 2.0+S1*S4*W25PU/2.0+S1*S3*W23PU/2.0-S1*S4*W22/2.0-S1*S3*W22/2.0 = +S1*S2*W22/2.0-S1*W22/2.0-S2*S4*W1P5P+S2*S4*W1P5-S1*S4*W1P5/2.0 > +S2*W1P5-S1*W1P5/2.0-S2*S3*W1P3P+S2*S3*W1P3-S1*S3*W1P3/2.0+S2*S ? 4*W15P-S1*S4*W15P/2.0-S4*W15P-S2*S4*W15+S4*W15-S2*W15+W15+S2*S3 @ *W13P-S1*S3*W13P/2.0-S3*W13P-S2*S3*W13+S3*W13+S1*S4*W11/2.0+S1* 1 S3*W11/2.0+S2**2*W11-S1*S2*W11/2.0-S2*W11+S1*W11/2.0+GX*S1*S4/2 2 .0+GX*S1*S3/2.0+GX*S1*S2/2.0-GX*S1**2/4.0+G5 C IF((XFE3TET.GT.0.0).AND. 1 (XFE2OCT.GT.0.0).AND.(XFE3OCT.GT.0.0)) THEN MAGNETITE = MAGNETITE + R*T*(LOG(XFE3TET) + LOG(XFE3OCT) + 1 LOG(XFE2OCT)) C Pure Magnetite: [x2=0, x3=0, x4=0, x5=1, s1=0, s2=0, s3=0] DX2 = 0.0D0 DX3 = 0.0D0 DX4 = 0.0D0 DX5 = 1.0D0 CALL ORD_SPN(T, DX2, DX3, DX4, DX5, DS1, DS2, DS3, DS4, 1 DUMMY1, DUMMY2, DUMMY3, DUMMY4, DUMMY5, DUMMY6, 2 DUMMY7, DUMMY8, DUMMY9, DUMMY10, DUMMY11, LSTOP) IF(LSTOP) WRITE(*,'(1X,''Failed in Mt call to ORD_SPN.'')') PURE = R*T*(LOG(DS4+1.0)+2.0*LOG(1.0-DS4)-2.0*LOG(2.0)) PURE = PURE + DS4**2*W55 + G5 MAGNETITE = EXP((MAGNETITE-PURE)/(R*T)) ELSE MAGNETITE = 0.0 END IF C MGTITANATE = 0.0 MGFERRITE = 0.0 CHROMITE = 0.0 MGCHROMITE = 0.0 C RETURN END FUNCTION DIFF(X,Y) C C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUNE 7 C TO APPEAR IN ^SOLVING LEAST SQUARES PROBLEMS^, PRENTICE-HALL, 1974 C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIFF = X - Y RETURN END PROGRAM GEOTHERMOMETER C C Fe-Ti oxide geothermometer C C Implements the solution models of: C C Ghiorso and Sack (1990) for rhombohedral oxides and Sack C and Ghiorso (1990) for spinels along with the standard state C database of Berman (1988). C IMPLICIT NONE C C Data in these common blocks is set in GET_DATA.FOR C DOUBLE PRECISION PRES, XIL, XGK, XPY, XSP, XCR, XUV, XMG COMMON /VALUES/ PRES, XIL, XGK, XPY, XSP, XCR, XUV, XMG C C Variables needed by GIBBS.FOR C CHARACTER*20 PHASE C C Variables used by ACT_SPN.FOR C DOUBLE PRECISION S1, S2, S3, S4, A XMG2TET, XFE2TET, XAL3TET, XFE3TET, XCR3TET, B XMG2OCT, XFE2OCT, XAL3OCT, XFE3OCT, XTI4OCT, C XCR3OCT, D HERCYNITE, SPINEL, ULVOSPINEL, MAGNETITE, E MGTITANATE, MGFERRITE, CHROMITE, MGCHROMITE C C Variables used by ACT_RHM.FOR C DOUBLE PRECISION ILMENITE, GEIKIELITE, HEMATITE, PYROPHANITE C C Local variables C LOGICAL LOOP PARAMETER(LOOP = .TRUE.) DOUBLE PRECISION BIG ! Square root of the largest D-floating PARAMETER(BIG = 1.0D19) DOUBLE PRECISION R ! R in kJ PARAMETER(R = 8.3143D0/1000.0D0) DOUBLE PRECISION DG0, FO2, ST, T INTEGER I, ITMAX, MODE LOGICAL LSTOP C C Declared functions C DOUBLE PRECISION GIBBS, ZERO EXTERNAL ZERO C C Loop endlessly through calculation loop until end-of-data is C detected and program is halted in GET_DATA C DO WHILE(LOOP) C C Acquire new data point (passes information in common blocks VALUES) C CALL GET_DATA C C Calculate T by minimizing the Gibbs energy of the reaction: C C Fe2TiO4 + Fe2O3 = Fe3O4 + FeTiO3 C WRITE(*,'(/,1X,''*****CALCULATING THE TEMPERATURE*****'')') PRES = 1.0 ! Pressure in bars T = 1200.0 ! Initial guess at equilibration temperature (K) ST = 100.0 ! Initial step increment in T ITMAX = 250 ! Maximum number of iterations C CALL MINIMA(ZERO, T, ST, BIG, S1, ITMAX, MODE) IF(S1.GT.1.0D-12) MODE = 6 ! The true minimum must be zero C IF(MODE.EQ.0) THEN C C Calculate an log f O2 from the reaction: C C 6 Fe2O3 = 4 Fe3O4 + O2 C WRITE(*,'(/,1X,''*****CALCULATING THE FO2*****'')') PHASE = 'O2 ' DG0 = GIBBS(PHASE, T, PRES) PHASE = 'MAGNETITE ' DG0 = DG0 + 4.0*GIBBS(PHASE, T, PRES) PHASE = 'HEMATITE ' DG0 = DG0 - 6.0*GIBBS(PHASE, T, PRES) CALL ACT_SPN(T, XSP, XCR, XUV, XMG, S1, S2, S3, S4, 1 XMG2TET, XFE2TET, XAL3TET, XFE3TET, XCR3TET, XMG2OCT, 2 XFE2OCT, XAL3OCT, XFE3OCT, XTI4OCT, XCR3OCT, 3 HERCYNITE, SPINEL, ULVOSPINEL, MAGNETITE, MGTITANATE, 4 MGFERRITE, CHROMITE, MGCHROMITE, LSTOP) IF(.NOT.LSTOP) THEN CALL ACT_RHM(T, XIL, XGK, XPY, S1, S2, S3, ILMENITE, 1 GEIKIELITE, HEMATITE, PYROPHANITE, LSTOP) IF(.NOT.LSTOP) THEN FO2 = -DG0/(R*T)+6.0*LOG(HEMATITE)-4.0*LOG(MAGNETITE) FO2 = FO2/LOG(10.0D0) ELSE MODE = 5 FO2 = 0.0 END IF C ELSE MODE = 4 FO2 = 0.0 END IF ELSE T = 0.0 FO2 = 0.0 END IF C C Output results C IF(MODE.EQ.0) THEN WRITE(*,'(/,1X,''*******************************'', 1 /,1X,''T (C): '', F7.2, '' log f O2: '', F6.2, 2 /,1X,''*******************************'')') 3 T - 273.15, FO2 ELSE WRITE(*,'(1X,''*****CONVERGENCE ERROR*****'')') END IF C C Return for more data C END DO C STOP END SUBROUTINE GET_DATA C IMPLICIT NONE C C Subroutine to acquire data (either from a file or from C the terminal) for the Fe-Ti oxide geothermometer C program C C The following common block passes the necessary data back C to the main program and to the minimization routines C DOUBLE PRECISION PRES, ! Pressure estimate for assemblage A XIL, ! Mole fraction of ilmenite in Rhm oxide B XGK, ! Mole fraction of geikielite in Rhm oxide C XPY, ! Mole fraction of pyrophanite in Rhm oxide D XSP, ! Mole fraction of spinel in cubic oxide E XCR, ! Mole fraction of chromite in cubic oxide F XUV, ! Mole fraction of ulvospinel in cubic oxide G XMG ! Mole fraction of magnetite in cubic oxide COMMON /VALUES/ PRES, XIL, XGK, XPY, XSP, XCR, XUV, XMG C CHARACTER*5 NAMES(11) ! Oxide names DOUBLE PRECISION WEIGHTS(11) ! Oxide molecular weights DOUBLE PRECISION CHARGE(11) ! Cation charges DOUBLE PRECISION CATIONS(11) ! Cation numbers DOUBLE PRECISION WTCUBIC(11), WTRHOM(11) C C Local variables C DOUBLE PRECISION SUM_C_SP, SUM_C_RH, FE3_SP, FE2_SP, FE3_RH, A FE2_RH, SUM_CHARGE_SP, SUM_CHARGE_RH, RATIO, X_FE3_SP, B X_FE2_SP, SUM, TI DOUBLE PRECISION MCUBIC(11), MRHOM(11) INTEGER I LOGICAL FIRST CHARACTER*10 INPUT C DATA NAMES / 'SiO2 ', 'TiO2 ', 'Al2O3', 'V2O3 ', 'Cr2O3', A 'FeO ', 'MnO ', 'MgO ', 'CaO ', 'ZnO ', B 'NiO ' / DATA WEIGHTS / 60.0848D0, 79.8988D0, 101.9612D0, 149.8822D0, A 151.9902D0, 71.8464D0, 70.9374D0, 40.3114D0, B 56.0794D0, 81.3694D0, 74.7094D0 / DATA CHARGE / 4.0D0, 4.0D0, 3.0D0, 3.0D0, 3.0D0, 2.0D0, 2.0D0, A 2.0D0, 2.0D0, 2.0D0, 2.0D0 / DATA CATIONS / 1.0D0, 1.0D0, 2.0D0, 2.0D0, 2.0D0, 1.0D0, 1.0D0, A 1.0D0, 1.0D0, 1.0D0, 1.0D0 / DATA WTCUBIC / 11*0.0D0 /, WTRHOM / 11*0.0D0 / DATA FIRST / .TRUE. / C C Get input data - Spinel first, Rhombohedral phase second C IF(.NOT.FIRST) THEN WRITE(*,'(/,1X,''Do another calculation (Y or N)? '',$)') READ(*,'(A10)') INPUT IF((INPUT(1:1).NE.'Y').AND.(INPUT(1:1).NE.'y')) STOP END IF FIRST = .FALSE. C WRITE(*,'(1X,''***** CUBIC OXIDE PHASE*****'')') DO I=1,11 WRITE(*,'(1X,''Input Wt % '',A5,'' in spinel [default: '', 1 F6.2,'']: '',$)') NAMES(I), WTCUBIC(I) READ(*,'(A10)') INPUT IF(INPUT.NE.' ') READ(INPUT,'(F10.0)') WTCUBIC(I) END DO WRITE(*,'(/,1X,''***** RHOMBOHEDRAL OXIDE PHASE*****'')') DO I=1,11 WRITE(*,'(1X,''Input Wt % '',A5,'' in rhom ox [default: '', 1 F6.2,'']: '',$)') NAMES(I), WTRHOM(I) READ(*,'(A10)') INPUT IF(INPUT.NE.' ') READ(INPUT,'(F10.0)') WTRHOM(I) END DO C C Convert spinel and rhombohedral oxide analyses to 3 cations and C 4 oxygens or 2 cations and 3 oxygens, respectively C SUM_C_SP = 0.0 SUM_C_RH = 0.0 DO I=1,11 MCUBIC(I) = WTCUBIC(I)/WEIGHTS(I) SUM_C_SP = SUM_C_SP + MCUBIC(I)*CATIONS(I) MRHOM(I) = WTRHOM(I)/WEIGHTS(I) SUM_C_RH = SUM_C_RH + MRHOM(I)*CATIONS(I) END DO SUM_C_SP = 3.0/SUM_C_SP SUM_C_RH = 2.0/SUM_C_RH DO I=1,11 MCUBIC(I) = MCUBIC(I)*CATIONS(I)*SUM_C_SP MRHOM(I) = MRHOM(I)*CATIONS(I)*SUM_C_RH END DO SUM_CHARGE_SP = 0.0 SUM_CHARGE_RH = 0.0 DO I=1,11 SUM_CHARGE_SP = SUM_CHARGE_SP + CHARGE(I)*MCUBIC(I) SUM_CHARGE_RH = SUM_CHARGE_RH + CHARGE(I)*MRHOM(I) END DO SUM_CHARGE_SP = SUM_CHARGE_SP - 2.0*MCUBIC(6) SUM_CHARGE_RH = SUM_CHARGE_RH - 2.0*MRHOM(6) FE3_SP = 8.0 - SUM_CHARGE_SP - 2.0*MCUBIC(6) FE2_SP = MCUBIC(6) - FE3_SP FE3_RH = 6.0 - SUM_CHARGE_RH - 2.0*MRHOM(6) FE2_RH = MRHOM(6) - FE3_RH C C mole fractions for Ghiorso and Sack (1991) model. C RATIO = (FE2_SP+MCUBIC(7)+MCUBIC(8)+MCUBIC(9)+MCUBIC(10) 1 + MCUBIC(11))/(FE2_SP+MCUBIC(8)) XSP = MCUBIC(8)*RATIO XCR = MCUBIC(5)/2.0 XUV = MCUBIC(2) XMG = FE3_SP/2.0 C C mole fractions for Ghiorso and Sack (1991) model. C SUM = FE3_RH/2.0 + MRHOM(2) XIL = (1.0-FE3_RH/(2.0*SUM)) 1 * (FE2_RH)/(FE2_RH+MRHOM(7)+MRHOM(8)) XGK = (1.0-FE3_RH/(2.0*SUM)) 1 * MRHOM(8)/(FE2_RH+MRHOM(7)+MRHOM(8)) XPY = (1.0-FE3_RH/(2.0*SUM)) 1 * MRHOM(7)/(FE2_RH+MRHOM(7)+MRHOM(8)) C RETURN END DOUBLE PRECISION FUNCTION GIBBS(PHASE, T, P) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C PARAMETER(R=8.3143D0) PARAMETER(TR = 298.15D0) PARAMETER(PR = 1.0D0) C DOUBLE PRECISION K0, K1, K2, K3, CP_T, CP_H, L1, L2 CHARACTER*20 PHASE C IF(PHASE.EQ.'MAGNETITE ') THEN HS = -1117403.0 ! J Berman (1988) SS = 146.114 ! J Berman (1988) VS = 4.452 ! J Berman (1988) K0 = 207.93 ! J Berman and Brown (1985) K1 = 0.0 K2 = -72.433D5 ! J K3 = 66.436D7 ! J CP_T = 848.0 CP_H = 1565.0 ! J L1 = -19.502D-2 ! J L2 = 61.037D-5 ! J V1 = -0.582D-6 ! Berman (1988) V2 = 1.751D-12 ! Berman (1988) V3 = 30.291D-6 ! Berman (1988) V4 = 138.470D-10 ! Berman (1988) ELSE IF(PHASE.EQ.'ULVOSPINEL ') THEN HS = -1488500.0 ! -1496000.0 ! J SS = 185.447 ! 174.70 J/K Todd and King (1953) + 2 R ln 2 VS = 4.682 ! J Robie et al. (1978) K0 = 249.63 ! J Berman and Brown (1985) K1 = -18.174D2 ! J K2 = 0.0 ! J K3 = -5.453D7 ! J CP_T = 0.0 CP_H = 0.0 L1 = 0.0 L2 = 0.0 V1 = -0.582D-6 ! Berman (1988) - Assumed Magnetite Properties V2 = 1.751D-12 ! Berman (1988) - Assumed Magnetite Properties V3 = 30.291D-6 ! Berman (1988) - Assumed Magnetite Properties V4 = 138.470D-10 ! Berman (1988) - Assumed Magnetite Properties ELSE IF(PHASE.EQ.'HEMATITE ') THEN HS = -822000.0 ! -825627.0 ! (-824640 +/- 1.255) Berman (1988) SS = 87.40 ! 87.437 ! (87.40 +/- 0.08) Berman (1988) VS = 3.027 ! J Berman (1988) K0 = 146.86 ! J Berman and Brown (1985) K1 = 0.0 K2 = -55.768D5 ! J K3 = 52.563D7 ! J CP_T = 955.0 CP_H = 1287.0 ! J L1 = -7.403D-2 ! J L2 = 27.921D-5 ! J V1 = -0.479D-6 ! Berman (1988) V2 = 0.304D-12 ! Berman (1988) V3 = 38.310D-6 ! Berman (1988) V4 = 1.650D-10 ! Berman (1988) ELSE IF(PHASE.EQ.'ILMENITE ') THEN HS = -1231947.0 ! (-1236620.0+/-1.59 Robie) Berman(1988) SS = 108.628 ! (108.90 +/- 0.50 Grenvold) Berman(1988) VS = 3.170 ! J Berman (1988) K0 = 150.00 ! J Berman (1988) K1 = -4.416D2 ! J K2 = -33.237D5 ! J K3 = 34.815D7 ! J CP_T = 0.0 CP_H = 0.0 L1 = 0.0 L2 = 0.0 V1 = -0.584D-6 ! Berman (1988) V2 = 1.230D-12 ! Berman (1988) V3 = 27.248D-6 ! Berman (1988) V4 = 29.968D-10 ! Berman (1988) ELSE IF(PHASE.EQ.'O2 ') THEN HS = 0.0 SS = 0.0 VS = 0.0 K0 = 0.0 K1 = 0.0 K2 = 0.0 K3 = 0.0 CP_T = 0.0 CP_H = 0.0 L1 = 0.0 L2 = 0.0 V1 = 0.0 V2 = 0.0 V3 = 0.0 V4 = 0.0 ELSE WRITE(*,*) 'Illegal phase in call to GIBBS!' GIBBS = 0.0 RETURN END IF C IF(PHASE.EQ.'O2 ') THEN HS = 23.10248*(T-TR) + 2.0*804.8876*(SQRT(T)-SQRT(TR)) 1 - 1762835.0*(1.0/T-1.0/TR) 2 - 18172.91960*LOG(T/TR) + 0.5*0.002676*(T**2-TR**2) SS = 205.15 + 23.10248*LOG(T/TR) 1 - 2.0*804.8876*(1.0/SQRT(T)-1.0/SQRT(TR)) 2 - 0.5*1762835.0*(1.0/T**2-1.0/TR**2) 3 + 18172.91960*(1.0/T-1.0/TR) + 0.002676*(T-TR) C VS = R*T*LOG(P/PR) VS = 0.0 ! Standard state is any T and 1 bar for gases ELSE 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) IF(CP_T.NE.0.0) THEN IF(T.GT.CP_T) THEN HS = HS + CP_H 1 + 0.5*L1**2*(CP_T**2-TR**2) 2 + (2.0/3.0)*L1*L2*(CP_T**3-TR**3) 3 + 0.25*L2**2*(CP_T**4-TR**4) SS = SS + CP_H/CP_T 1 + L1**2*(CP_T-TR) 2 + L1*L2*(CP_T**2-TR**2) 3 + (1.0/3.0)*L2**2*(CP_T**3-TR**3) ELSE HS = HS + 0.5*L1**2*(T**2-TR**2) 1 + (2.0/3.0)*L1*L2*(T**3-TR**3) 2 + 0.25*L2**2*(T**4-TR**4) SS = SS + L1**2*(T-TR) 1 + L1*L2*(T**2-TR**2) 2 + (1.0/3.0)*L2**2*(T**3-TR**3) END IF END IF 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)) END IF C GS = HS - T*SS + VS GIBBS = GS/1000.0 ! Convert to kj C RETURN END SUBROUTINE MINIMA(S,B,ST,BIG,S1,ITMAX,MODE) C C IMPLEMENTATION OF ALGORITHM 17. C C SUCCESS-FAILURE LINEAR SEARCH WITH PARABOLIC INVERSE C INTERPOLATION C C FROM J.C. NASH (1979) C COMPACT NUMERICAL METHODS FOR COMPUTERS:LINEAR ALGEBRA C AND FUNCTION MINIMISATION. JOHN WILEY AND SONS, NEW C YORK, 227 PAGES. C C THIS PROCEDURE IS DESIGNED TO FIND A LOCAL MINIMUM OF THE C FUNCTION S(B) WITH RESPECT TO B. C C MARK S. GHIORSO C DEPARTMENT OF GEOLOGICAL SCIENCES, AJ-20 C UNIVERSITY OF WASHINGTON C SEATTLE, WA 98195 C C VERSION 1.0 JANUARY 1983 C IMPLEMENTED WITH A LOGICAL VARIABLE STEP8 WHICH C IS SET TO INCLUDE (.TRUE.) OR EXCLUDE (.FALSE.) C ALGORITHMIC STEP NUMBER EIGHT. C C ------------------------------------------------------------------ C C VARIABLES: C C S (B,LSTOP) FUNCTION TO BE MINIMIZED. MUST BE DECLARED EXTERNAL C IN THE CALLING PROGRAM. B IS THE FUNCTION PARAMETER. LSTOP C IS A LOGICAL VARIABLE WHICH ASSUMES THE VALUE .FALSE. IF C B IS A VALID PARAMETER VALUE AND .TRUE. IF B IS INVALID. S C MUST BE WRITTEN BY THE USER. C C B INPUT: INITIAL GUESS TO THE MINIMUM OF S(B) C OUTPUT: VALUE OF THE FUNCTION PARAMETER WHICH MINIMIZES S C IF MODE EQUALS ZERO. IF MODE EQUALS ONE B RETURNS C THE INITIAL GUESS. C C ST INPUT: INITIAL STEP SIZE. MODIFIED ON OUTPUT. C C BIG A NUMBER SUCH THAT -BIG IS CLEARLY LESS THAN ANY RESULT OF C THE COMPUTATION OF THE FUNCTION S(B). UNCHANGED ON OUTPUT. C C S1 OUTPUT: VALUE OF S(B) IF MODE EQUALS ZERO (I.E.B IS A C MINIMIZER). C C ITMAX INPUT: MAXIMUM NUMBER OF FUNCTION EVALUATIONS IN THE C MINIMIZATION ATTEMPT. C OUTPUT: NUMBER OF CALLS TO S USED IN THE MINIMIZATION C ATTEMPT IF MODE EQUALS ZERO, OTHERWISE UNCHANGED. C C MODE OUTPUT: 0 IF MINIMIZATION IS SUCCESSFUL. C 1 IF INITIAL POINT IS INFEASIBLE. C 2 IF ITERATION LIMIT IS EXCEEDED. C C ----------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) C PARAMETER(TAU = 1.0D-14) C LOGICAL LSTOP,STEP8,DEBUG C COMMON /STATUS/ DEBUG C C CHOSE TO AVOID ALGORITHM STEP NUMBER 8. C THIS ALWAYS ASSURES CONVERGENCE, BUT IT MAY SLOW DOWN THE RATE C STEP8 = .FALSE. C C SET DEFAULT MODE PARAMETER RETURN (ERROR CONDITION OF ITERATION C LIMIT EXCEEDED. C MODE = 2 C C STEP 0 ENTER B, THE INITIAL VALUE OF THE PARAMETER WITH RESPECT C TO WHICH THE FUNCTION IS TO BE MINIMIZED. C C B = C C ENTER ST, THE INITIAL STEP LENGTH (CAN BE NEGATIVE) C C ST = C C LET IFN = 0, TO INITIALIZE THE FUNCTION EVALUATION COUNTER C IFN = 0 C C LET A1 = 1.5, LET A2 = -0.25, THE SUCCESS AND FAILURE MULTIPLIERS C OF THE STEP LENGTH C A1 = 1.5 A2 = -0.25 C C ENTER BIG, A NUMBER SUCH THAT -BIG IS CLEARLY LESS THAN ANY RESULT C OF THE COMPUTATION OF THE FUNCTION S(B). C C BIG = C C STEP 1 LET IFN = IFN + 1, TO COUNT FUNCTION EVALUATIONS C IFN = IFN + 1 C C COMPUTE P=S(B). IF FUNCTION CANNOT BE COMPUTED, STOP. C P = S(B,LSTOP) IF(.NOT.LSTOP) GO TO 10 MODE = 1 RETURN C C STEP 2 LET S1 = P, THE MINIMUM FUNCTION VALUE SO FAR C 10 S1 = P C C LET S0 = -BIG, TO INDICATE THAT A VALUE FOR S0 HAS YET TO BE C CALCULATED C S0 = -BIG C C LET X1 = 0, LET BMIN = B, TO SAVE POSITION OF MINIMUM POINT SO FAR C X1 = 0.0 BMIN = B C C STEP 3 LET X2 = X1 + ST, TO ADJUST THE STEP C 15 X2 = X1 + ST C C LET B = BMIN + X2, TO CHANGE THE PARAMETER C B = BMIN + X2 C C STEP 4 IF B = BMIN + X1, STOP, AS THE VALUE OF THE PARAMETER IS C NOT ALTERED FROM THE MINIMUM POSITION. THE FUNCTION S(B) C HAS A (LOCAL) MINIMUM VALUE OF S1 AT B. THIS TEST MAY FAIL C IF BMIN IS VERY SMALL, IN WHICH CASE THE TEST IF ABS(B) + C EPS = ABS(BMIN) + ABS(X1) + EPS, MAY BE SUBSTITUTED WHERE C EPS IS THE MACHINE PRECISION. C C THE LATTER SUGGESTION IS IMPLEMENTED WITH LAWSON AND HANSON"S C SUBROUTINE DIFF. C C IF(DIFF(BMIN+X1,B).NE.0.0) GO TO 20 IF(ABS(DIFF(BMIN+X1,B)).GE.SQRT(TAU)) GO TO 20 P = S(B,LSTOP) IF(LSTOP) GO TO 30 S1 = P ! Modification 10/6/1989 for unusual minima case MODE = 0 ITMAX = IFN RETURN C C STEP 5 LET IFN = IFN + 1, TO COUNT FUNCTION EVALUATIONS. CHECK TO C MAKE SURE ITERATION LIMIT HAS NOT BEEN EXCEEDED. C 20 IFN = IFN + 1 IF(IFN.GT.ITMAX) RETURN C C COMPUTE P = S(B). IF THE FUNCTION CANNOT BE COMPUTED GO TO STEP 9 C (FAILURE). C P = S(B,LSTOP) IF(LSTOP) GO TO 30 C C STEP 6 IF P < S1, GO TO STEP 10 (SUCCESS). EQUALITY IS COUNTED AS C A FAILURE. C IF(P.LT.S1) GO TO 40 C C STEP 7 IF S0 >= S1, GO TO STEP 11 TO PERFORM THE INVERSE INTER- C POLATION SINCE THERE ARE NOW THREE POINTS AVAILABLE. THIS C IS THE REASON FOR SETTING S0 IN STEP 2, THOUGH OTHER C TECHNIQUES FOR INDICATING THAT THREE POINTS HAVE BEEN C FOUND IN A V-SHAPE MAY BE USED. C IF(S0.GE.S1) GO TO 50 C C STEP 8 SAVE POINT, RESETTING S0, THIS STEP MAY BE DELETED IF THE C (INITIAL POINT, FAILURE, FAILURE) CASE IS TO BE EXCLUDED C FROM THE SET OF CONDITIONS WHICH PRECEDE INVERSE INTER- C POLATION. NOTE THAT THIS STEP IS SKIPPED IF THE FAILURE C IS DUE TO NON-COMPUTABILITY OF THE FUNCTION, IN WHICH CASE C THE STEP SIZE IS SIMPLY REDUCED. C C LET S0 = P, LET X0 = X2 C IF(.NOT.STEP8) GO TO 30 S0 = P X0 = X2 C C STEP 9 ACTION IN CASE OF "FAILURE" C LET ST = A2*ST TO REDUCE SIZE AND CHANGE SIGN OF STEP LENGTH C 30 ST = A2*ST C C GO TO STEP 3 TO TRY ANOTHER POINT C GO TO 15 C C STEP 10 ACTION IN CASE OF "SUCCESS". C LET X0 = X1, LET S0 = S1, TO SAVE LAST POINT C 40 X0 = X1 S0 = S1 C C LET X1 = X2, LET S1 = P, TO SAVE NEW MINIMUM POINT C X1 = X2 S1 = P C C LET ST = A1*ST TO INCREASE THE STEP LENGTH C ST = A1*ST C C GO TO STEP 3 TO TRY ANOTHER POINT C GO TO 15 C C STEPP 11 PREPARE FOR PARABOLIC INVERSE INTERPOLATION C LET X0 = X0 - X1, TO CALCULATE X0 WHICH HAS BEEN C STORING B SUB 0. C 50 X0 = X0 - X1 C C LET S0 = (S0-S1)*ST, LET P = (P-S1)*X0/ C NOTE: P CONTAINS S SUB 2. C S0 = (S0-S1)*ST P = (P-S1)*X0 C C STEP 12 IF P = S0, GO TO STEP 18. SAFETY CHECK TO AVOID THE INVERSE C INTERPOLATION IF IT WOULD RESULT IN A DIVISION BY ZERO. C THIS CAN ONLY ARISE IN THE (INITIAL POINT,FAILURE,FAILURE) C CASE WHEN THE FUNCTION IS FLAT, THAT IS, IN THE EXCEPTION C MENTIONED AFTER EQUATION (13.9). C IF(P.EQ.S0) GO TO 60 C C STEP 13 LET ST = 0.5*(P*X0-S0*ST)/(P-S0) TO EVALUATE EQUATION (13.14) C ST = 0.5*(P*X0 - S0*ST) / (P-S0) C C STEP 14 LET X2 = X1 + ST. LET B = BMIN + X2, TO ADJUST STEP AND C PARAMETER. C X2 = X1 + ST B = BMIN + X2 C C STEP 15 IF B = BMIN + X1, GO TO STEP 20, SINCE INVERSE INTERPOLATION C DOES NOT ALTER PARAMETER FROM THE MINIMUM POSITION. C IF(DIFF(B,BMIN+X1).EQ.0.0) GOTO 80 C C STEP 16 LET IFN = IFN + 1, TO COUNT FUNCTION EVALUATIONS C IFN = IFN + 1 C C COMPUTE P = S(B). IF THE FUNCTION CANNOT BE COMPUTED, GO TO STEP 18 C P = S(B,LSTOP) IF(LSTOP) GO TO 60 C C STEP 17 IF P < S1, GO TO STEP 19, AS THE INVERSE INTERPOLATION C HAS FOUND A NEW MINIMUM POINT C IF(P.LT.S1) GO TO 70 C C STEP 18 LET B = BMIN + X1, LET P = S1, TO RESET THE MINIMUM VALUE C AND POSITION. C 60 B = BMIN + X1 P = S1 C C GO TO STEP 20. NOTE X1 IS ZERO IF ALL STEPS ARE FAILURES TO THIS C POINT. C GO TO 80 C C STEP 19 LET X1 = X2, SO THAT X1 GIVES THE STEP FROM BMIN TO THE C MINIMUM POSITION. C 70 X1 = X2 C C STEP 20 LET ST = A2*ST TO REDUCE THE STEP LENGTH (THE SIGN SHOULD C BE IRRELEVANT!). NOTE THAT THE NEW ST IS ZERO IF ALL STEPS C ARE FAILURES TO THIS POINT. THIS WILL LEAD TO CONVERGENCE C ON THE NEXT CYCLE. THUS THE SEQUENCE (INITIAL POINT, C FAILURE, FAILURE) FOLLOWED BY AN UNSUCCESSFUL INVERSE C INTERPOLATION WILL POSSIBLY CAUSE THE ALGORITHM TO STOP C BEFORE AN ACCURATE APPROXIMATION TO THE MINIMUM HAS C BEEN FOUND. DELETING STEP 8 WILL AVOID THIS EARLY CON- C VERGENCE BUT MAY INCREASE THE WORK TO LOCATE A MINIMUM C WHEN THE INITIAL POINT IS A GOOD APPROXIMATION TO A C MINIMUM. C 80 ST = A2*ST C C GO TO STEP 2 TO RESTART THE SUCCESS-FAILURE CYCLE C GO TO 10 C C NOTE THAT THE ALGORITHM TERMINATES VIA A SEQUENCE OF OPERATIONS C WHICH REDUCE THE STEP SIZE UNTIL THE TEST AT STEP 4 IS SATIS- C FIED. C 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) THEN WRITE(*,1000) 1000 FORMAT(1X,'Iteration limit (1000) exceeded in ORD_SPN!') LSTOP = .TRUE. END IF 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 C WRITE(*,'(1X,''s1, s1old: '', 3G20.13)') S1, S1OLD, C 1 ABS(S1-S1OLD) C WRITE(*,'(1X,''s2, s2old: '', 3G20.13)') S2, S2OLD, C 1 ABS(S2-S2OLD) C WRITE(*,'(1X,''s4, s4old: '', 3G20.13)') S4, S4OLD, C 1 ABS(S4-S4OLD) C END DO C C WRITE(*,*) 'Iterations = ', J 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 ZERO(T, LSTOP) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C COMMON /VALUES/ P, XIL, XGK, XPY, XSP, XCR, XUV, XMG CHARACTER*20 PHASE DOUBLE PRECISION MAGNETITE, MGTITANATE, MGFERRITE, ILMENITE, 1 MGCHROMITE LOGICAL LSTOP PARAMETER(R = 8.3143D0/1000.0D0) ! R in kJ C LSTOP = .TRUE. IF(P.EQ.0.0) P = 1.0 C C Catch run-away values of temperature C IF((T.LT.273.15D0).OR.(T.GT.1773.15D0)) RETURN C C Evaluate standard state Gibbs energy of the reaction: C C Fe2TiO4 + Fe2O3 = Fe3O4 + FeTiO3 C PHASE = 'MAGNETITE ' DG0 = GIBBS(PHASE, T, P) ! Gibbs energy is in kJ PHASE = 'ILMENITE ' DG0 = DG0 + GIBBS(PHASE, T, P) PHASE = 'ULVOSPINEL ' DG0 = DG0 - GIBBS(PHASE, T, P) PHASE = 'HEMATITE ' DG0 = DG0 - GIBBS(PHASE, T, P) C C Evaluate activities of components in phases C CALL ACT_SPN(T, XSP, XCR, XUV, XMG, S1, S2, S3, S4, 1 XMG2TET, XFE2TET, XAL3TET, XFE3TET, XCR3TET, XMG2OCT, 2 XFE2OCT, XAL3OCT, XFE3OCT, XTI4OCT, XCR3OCT, 3 HERCYNITE, SPINEL, ULVOSPINEL, MAGNETITE, MGTITANATE, 4 MGFERRITE, CHROMITE, MGCHROMITE, LSTOP) IF((MAGNETITE.EQ.0.0).OR.(ULVOSPINEL.EQ.0.0).OR.LSTOP) RETURN C CALL ACT_RHM(T, XIL, XGK, XPY, S1, S2, S3, ILMENITE, 1 GEIKIELITE, HEMATITE, PYROPHANITE, LSTOP) IF((ILMENITE.EQ.0.0).OR.(HEMATITE.EQ.0.0).OR.LSTOP) RETURN C C Evaluate total Gibbs energy change for the reaction C DG = R*T*LOG(MAGNETITE) DG = DG + R*T*LOG(ILMENITE) DG = DG - R*T*LOG(ULVOSPINEL) DG = DG - R*T*LOG(HEMATITE) C C Correction for excess volumes of Uv-Mg and Il-Hm C WV1 = -1.250/10000.0 ! kJ XMg^2*XUv term WV2 = 1.018/10000.0 ! kJ XMg*XUv^2 term WV = -0.414/10000.0 ! kJ Xil*Xhm term DGV = ( 3.0*(WV2-WV1)*XUV**2 - 2.0*(WV2-WV1)*XUV 1 + WV1*(2.0*XUV-1.0) )*(P-1.0) 2 + WV*(1.0-2.0*XIL)*(P-1.0) DG = DG + DGV C DG = DG0 + DG C C Square the result so that deviations from equilibrium are C indicated by a positive quadratic function centered on C T-equilibrium C ZERO = DG**2 LSTOP = .FALSE. C WRITE(*,1000) T-273.15, ZERO C1000 FORMAT(1X,'ZERO: T (C) = ',F8.2,5X,'ZERO = ',G12.6) C RETURN END