SUBROUTINE CA_OLIVINE_CALC(TK,PP,P,R,MU_FO,MU_FA,MU_MGFE,MU_MO, 1 ACT_FO,ACT_FA,ACT_MO) C *************************************************************** C SUBROUTINE TO CALCULATE CHEMICAL POTENTIALS FOR CA-MG-FE OLIVINES C M. HIRSCHMANN AMERICAN MINERALOGIST 1991 76:1232-1248 C *************************************************************** C INPUTS: C TK (REAL*8) T IN KELVINS C PP (REAL*8) P IN BARS C P (REAL*8) 2* X CA2SIO4 C R (REAL*8) 2* X FE2SIO4 - 1 C OUTPUTS: (ALL IN JOULES/GFW) C MU_FO (REAL*8) CHEMICAL POTENTIAL OF Mg2SiO4 (2 ATOM BASIS) C MU_FA (REAL*8) CHEMICAL POTENTIAL OF Fe2SiO4 (2 ATOM BASIS) C MU_MO (REAL*8) CHEMICAL POTENTIAL OF CaMgSiO4 (2 ATOM BASIS) C MU_MGFE (REAL*8) MgFe EXCHANGE POTENTIAL (1 ATOM BASIS) C ACT_FO (REAL*8) ACTIVITY OF MgSi(0.5)O2 (1 ATOM BASIS) C ACT_FA (REAL*8) ACTIVITY OF FeSi(O.5)O2 (1 ATOM BASIS) C ACT_FO (REAL*8) ACTIVITY OF CaMgSiO4 (2 ATOM BASIS) C *************************************************************** IMPLICIT DOUBLE PRECISION(A-H,O-Z) IMPLICIT INTEGER(I-N) C REAL*8 MU_FO,MU_FA,MU_MGFE,MU_MO REAL*8 K0,K1,K2,K3,L1,L2 REAL*8 MD(7) DIMENSION XM(2,3),STOICH(20) CHARACTER*12 SILICA_PHASE RR=8.31437 C ***************************************************************** C SOLUTION PARAMETERS IN MODULAR FORM: C MD(1)=(GX+2WFEMG)/4 C MD(2)=(GX-2WFEMG)/4 C MD(3)=(GEX-W1FEMG+W2FEMG)/4 C MD(4)=(GEX+W1FEMG-W2FEMG)/4 C MD(5)=(F0+2.0*(WCAMG+WCAFE))/4 C MD(6)=(F0-2.0*(WCAMG+WCAFE))/4 C MD(7)=WCAMG MD(1)=0.050836E+05 MD(2)=0.000000E+05 MD(3)=0.000100E+05 MD(4)=0.000100E+05 F0= 0.095000E+05 WCAMG=0.345000E+05 WCAFE=0.219000E+05 DCAMGDP=0.350000E+00 DCMD1DP =0.007500E+00 WCAMG=WCAMG+DCAMG_DP*PP MD(1)=MD(1)+DMD1_DP*PP MD(5)=0.25*(F0+2.0*(WCAMG+WCAFE)) MD(6)=0.25*(F0-2.0*(WCAMG+WCAFE)) MD(7)=WCAMG C ****************************************************************** C CALCULATE END MEMBER STANDARD STATE FREE ENERGIES C (ALL FROM BERMAN, 1988) C ** FORSTERITE (BERMAN, 1988) H= -2174.420*1000.0 S= 94.010 V= 4.366 CP_T= 0.0 CP_H= 0.0 K0= 238.64 K1= -20.013 D+02 K2= 0.0 D+05 K3= -11.624 D+07 L1= 0.0 D-02 L2= 0.0 D-05 V1= -.791 D-06 V2= 1.351 D-12 V3= 29.464 D-06 V4= 88.633 D-10 CALL G_INTEGRAL(G_FO,PP,TK,H,S,V,K0,K1 1 ,K2,K3,CP_H,CP_T,L1,L2 2 ,V1,V2,V3,V4) C FAYALITE (FROM BERMAN, 1988) H= -1479360. S= 150.930 V= 4.630 K0= 248.93 K1= -19.239 D+02 K2= 0.0 K3= -13.910 D+07 CP_T= 0.0 CP_H= 0.0 L1 = 0.0 L2 = 0.0 V1 = -0.730 D-06 V2 = 0.0 V3 = 26.546 D-06 V4 = 79.482 D-10 CALL G_INTEGRAL(G_FA,PP,TK,H,S,V,K0,K1 1 ,K2,K3,CP_H,CP_T,L1,L2 2 ,V1,V2,V3,V4) C ** MONTICELLITE (BERMAN, 1988) H= -2250.027*1000.0 S= 108.300 V= 5.148 CP_T= 0.0 CP_H= 0.0 K0= 226.34 K1= -15.427 D+02 K2= -11.797 D+05 K3= -2.329 D+07 L1= 0.0 D-02 L2= 0.0 D-05 V1= -.904 D-06 V2= 2.000 D-12 V3= 27.863 D-06 V4= 76.339 D-10 CALL G_INTEGRAL(G_MO,PP,TK,H,S,V,K0,K1 1 ,K2,K3,CP_H,CP_T,L1,L2 2 ,V1,V2,V3,V4) C ****************************************************************** C CALL SUBROUTINE TO CALCULATE ORDERING STATE AT GIVEN T CALL CA_ORDER(TK,P,R,MD,S,XM) C XM(1,1)=M1 FE C XM(2,1)=M2 FE C XM(1,2)=M1 MG C XM(2,2)=M2 MG C XM(1,3)=M1 CA C XM(2,3)=M2 CA DO I=1,2 DO J=1,3 IF (XM(I,J).EQ.0.0) XM(I,J) = 1.0E-12 END DO END DO C CALCULATE CHEMICAL POTENTIALS IF(P.NE.0) THEN MU_MO=G_MO+RR*TK*DLOG(XM(1,2)*XM(2,3))+MD(3)-MD(6) 1 +P*(-MD(7)+MD(6)-MD(3)+MD(1))+MD(7)*P**2 2 +R*(MD(1)+MD(3)-MD(6)-MD(7)+P*(MD(7)+MD(6)-MD(3)+MD(1))) 3 +MD(1)*R**2+S*(MD(2)-MD(3)+MD(5)-MD(7) 4 +P*(MD(7)-MD(5)+MD(4)-MD(2))+R*(MD(4)-MD(3)))-MD(2)*S**2 ACT_MO=DEXP((MU_MO-G_MO)/(RR*TK)) ELSE MU_MO=0.0 ACT_MO=0.0 END IF XFO=1.0-0.5*((R+1.0)+P) IF (DABS(XFO).LT.1.0D-10) XFO=0.0 IF(XFO.GT.0.0) THEN MU_FO=G_FO+RR*TK*DLOG(XM(1,2)*XM(2,2))+MD(1) 1 +P*(MD(1)-MD(3)+MD(6)+MD(7))+P**2*MD(7) 2 +R*(2.0*MD(1)+P*(MD(1)-MD(3)+MD(6)+MD(7)))+MD(1)*R**2 3 +S*(-MD(3)+MD(4)+P*(-MD(2)+MD(4)-MD(5)+MD(7)) 4 +R*(-MD(3)+MD(4)))-MD(2)*S**2 C TWO-ATOM BASIS ACTIVITY: ACT_FO=DEXP((MU_FO-G_FO)/(RR*TK)) C ONE ATOM BASIS ACTIVITY: ACT_FO=DSQRT(ACT_FO) ELSE MU_FO=0.0 ACT_FO=0.0 END IF IF (R.NE.-1.0) THEN MU_FA=G_FA+RR*TK*DLOG(XM(1,1)*XM(2,1))+MD(1) 1 +P*(-MD(1)+MD(3)-MD(6)-MD(7))+MD(7)*P**2 2 +R*(-2.0*MD(1)+P*(MD(1)-MD(3)+MD(6)+MD(7)))+MD(1)*R**2 3 +S*(MD(3)-MD(4)+P*(-MD(2)+MD(4)-MD(5)+MD(7))) 4 +R*(-MD(3)+MD(4))-MD(2)*S**2 C TWO-ATOM BASIS ACTIVITY: ACT_FA=DEXP((MU_FA-G_FA)/(RR*TK)) C ONE ATOM BASIS ACTIVITY: ACT_FA=DSQRT(ACT_FA) ELSE MU_FA=0.0 ACT_FA=0.0 END IF IF(R.NE.-1.0.AND.XFO.NE.0.0) THEN MU_MGFE=0.5*(G_FO-G_FA) 1 +RR*TK*DLOG(XM(1,2)*XM(2,2)/XM(1,1)/XM(2,1)) 2 +P*(MD(1)-MD(3)+MD(6)+MD(7))+2.0*R*MD(1)+S*(MD(4)-MD(3)) ELSE MU_MGFE=0.0 END IF RETURN END SUBROUTINE CA_ORDER(TK,P,R,MD,S,XM) C SUBROUTINE TO CALCULATE HOMOGENEOUS EQUILIBRIUM IN C CA-MG-FE OLIVINES C C HIRSCHMANN, M. AMERICAN MINERALOGIST 1991 C C WRITTEN AND MODIFIED 1990, 1991 IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) DIMENSION XM(2,3) REAL*8 MD(7),K1,K2 LOGICAL CONVERGE,PASS,SA,SB RR=8.31437 TOLERANCE=1.0E-10 NMAX=200 C ** CALCULATE EQUILIBRIUM ORDERING STATE: C ******************************************************************** C ERROR TRAP FOR COMPOSITIONS OUT OF RANGE: IF (P.GT.1.0.OR.P.LT.0.0) THEN WRITE(*,50) P 50 FORMAT(1X,'P OUT OF RANGE, = ',E16.10) STOP END IF IF (R.GT.1.0.OR.R.LT.-1.0) THEN WRITE(*,50) R 51 FORMAT(1X,'R OUT OF RANGE, = ',E16.10) STOP END IF C ******************************************************************** B1=2.0 B2=2.0*P-2.0 C1=1.0-R**2 C2=1.0-R**2-2.0*P*R-2.0*P C ******************************************************************** C ERROR TRAP FOR TRIVIAL CASES XFA=(R+1.0)/2.0 XCA=P/2.0 XFO=1.0-XFA-XCA IF (DABS(XFO).LT.1.0D-12) XFO=0.0 IF(XFA.EQ.0.0) THEN S=0.0 XM(1,1)=1.0E-12 !M1 FE XM(2,1)=1.0E-12 !M2 FE XM(1,2)=1.0 !M1 MG XM(2,2)=1.0-P !M2 MG XM(1,3)=1.0E-12 !M1 CA XM(2,3)=P !M2 CA RETURN END IF IF(XFO.EQ.0.0) THEN S=-P XM(1,1)=1.0 !M1 FE XM(2,1)=1.0-P !M2 FE XM(1,2)=1.0E-12 !M1 MG XM(2,2)=1.0E-12 !M2 MG XM(1,3)=1.0E-12 !M1 CA XM(2,3)=P !M2 CA RETURN END IF IF (P.EQ.1.0) THEN S=-1.0-R RETURN END IF C ******************************************************************** C ** FIRST GUESS FOR S: SMAX=1.0 IF (SMAX.GT.(1.0-R-2.0*P)) SMAX = (1.0-R-2.0*P) IF (SMAX.GT.R+1.0) SMAX=R+1.0 SMIN=-1.0 IF (SMIN.LT.(-R-1.0)) SMIN=-R-1.0 IF (SMIN.LT.(R-1.0)) SMIN=R-1.0 IF (SMIN.GE.SMAX) THEN WRITE(*,20)P,R,SMIN,SMAX 20 FORMAT(1X,'ERROR IN CA-OLIVINE ORDERING ROUTINE',/, 1 1X,'SMIN IS NOT LESS THAN SMAX! ',4(E12.6,1X)) STOP END IF !SMIN.GE.SMAX S=(SMIN+SMAX)/2.0 C ******************************************************************** C ** ITERATIVE LOOP: CONVERGE=.FALSE. NCOUNT=0 DO WHILE(CONVERGE.EQ.FALSE) K1=2.0*MD(2)*S+(MD(3)-MD(4))*R+(-MD(7)+MD(5)-MD(4)+MD(2))*P 1 +MD(4)+MD(3) K2=DEXP(-2.0*K1/(RR*TK)) A=1.0-K2 B=B1-K2*B2 C=C1-K2*C2 ROOT=B**2-4.0*A*C IF (R00T.LT.0.0) THEN WRITE(*,1001) 1001 FORMAT(1X,'IMAGINARY ROOT IN CA-OLIVINE ORDERING ROUTINE. STOP.') STOP END IF S1A=(-B+DSQRT(ROOT))/(2.0*A) S1B=(-B-DSQRT(ROOT))/(2.0*A) C ******************************************************************** C ALGORITHM TO CHOOSE APPROPRIATE ROOT: PASS=.FALSE. DO WHILE (.NOT.PASS) SA=.TRUE. SB=.TRUE. IF (S1A.GT.SMAX) SA=.FALSE. IF (S1B.GT.SMAX) SB=.FALSE. IF (S1A.LT.SMIN) SA=.FALSE. IF (S1B.LT.SMIN) SB=.FALSE. IF (SA.OR.SB) THEN PASS=.TRUE. ELSE IF (.NOT.SA) S1A=(S+S1A)/2.0 IF (.NOT.SB) S1B=(S+S1B)/2.0 END IF END DO IF(SA.AND.SB) THEN TEST1=DABS(S1A-S) TEST2=DABS(S1B-S) IF (TEST1.LT.TEST2) THEN S1=S1A ELSE S1=S1B ENDIF ELSE IF(SA) S1=S1A IF(SB) S1=S1B END IF C ******************************************************************** C ** TEST FOR CONVERGENCE: DIFF=DABS(S1-S) IF (DIFF.LT.TOLERANCE) THEN CONVERGE=.TRUE. ELSE S=S1 NCOUNT=NCOUNT+1 IF (NCOUNT.GT.NMAX) THEN IF (NMAX.LT.450) THEN TOLERANCE=TOLERANCE*10.0 NMAX=NMAX+50 ELSE WRITE(*,1002)NCOUNT,DIFF 1002 FORMAT(1X,'CA_ORDERING ROUTINE WLL NOT CONVERGE. '/, 1 1X,I3,' LOOPS. DIFF = ',E12.6) STOP END IF !NMAX IF-THEN END IF !NCOUNT IF-THEN END IF !DIFF IF-THEN END DO ! END OF CONVERGE DO-WHILE C ******************************************************************** XM(1,1)=(R-S+1.0)/2.0 !M1 FE XM(2,1)=(R+S+1.0)/2.0 !M2 FE XM(1,2)=(S-R+1.0)/2.0 !M1 MG XM(2,2)=(1.0-R-S)/2.0-P !M2 MG XM(1,3)=0.0 !M1 CA XM(2,3)=P !M2 CA DO J=1,2 DO I=1,2 IF (DABS(XM(I,J)).LT.1.0E-12) XM(I,J)=1.0E-12 IF (XM(I,J).LE.0.0) THEN WRITE(*,400)I,J,XM(I,J) 400 FORMAT(1X,'ERROR IN CA-OLIVINE ORDERING ROUTINE',/, 1 1X,'XM',I1,',',I1,' = ',E12.6) STOP END IF END DO END DO RETURN END SUBROUTINE G_INTEGRAL(G_0,P,T,HS,SS,VS,K0,K1,K2,K3, 1 CP_H,CP_T,L1,L2,V1,V2,V3,V4) C CALCULATES STD_STATE G IMPLICIT REAL*8(A-H,O-Z) REAL*8 K0,K1,K2,K3,L1,L2 TR=298.15 PR=1.0 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)) G_0 = HS - T*SS + VS RETURN END