Program Fugco
* ---------------------------------------------------------------------|
* Calculation of fugacity coefficients for mixtures of H2O, CO2,
* and NaCl. Warning葉his program should not be used at temperatures
* below 350 deg C, pressures below 500 bars, or pseudobinaries of
* NaCl above 35 weight percent. The upper temperature limit is
* roughly 600 deg C.
* Bowers, T.S., Helgeson, H.C., 1985. Fortran programs for generating
* fluid inclusion isochores and fugacity coefficients for the system
* H2O-CO2-NaCl at high pressures and temperatures. Comput. Geosci. 11,
* 203-213.
* RELATED: Bowers, T.S., Helgeson, H.C., 1983. Calculation of the
* thermodynamic and geochemical consequences of nonideal mixing in the
* system H2O-CO2-NaCl on phase relations in geologic systems; equation
* of state for H2O-CO2-NaCl fluids at high pressures. Geochim.
* Cosmochim. Acta 47, 1247-1275.
* Bowers, T.S., Helgeson, H.C., 1983. Calculation of the thermodynamic
* and geochemical consequences of nonideal mixing in the system
* H2O-CO2-NaCl on phase relations in geologic systems; metamorphic
* equilibria at high pressures and temperatures. Am. Min. 68, 1059-1075.
* X = Mole fraction, GAM = Fugacity coefficient,
* R = Fugacity = X*GAM*PBAR. Species order is 1=CO2, 2=H2O, 3=NaCl
* ---------------------------------------------------------------------|
real*8 x(5), gam(5), r(3), tk, pat
character*12 blank
data blank/' '/
in1 = 10
iout1 = 11
open ( in1, file = 'furun.dat', status = 'old')
open (iout1, file = 'fuout.dat', status = 'unknown')
* Input minimum, maximum, and incremental values of temperature and
* pressure.
write (*,*)
write (*,*)
write (*,'(a15)') " Program Fugco"
write (*,'(a54)')
+ " Computes fugacity coefficients at P > 500 bar and"
write (*,'(a47)')
+ " T from 350-600 C in the sytem H2O-CO2-NaCl"
write (*,*)
write (*,'(a16)') " I/O file names"
write (*,'(a31)') " -Input file = ""furun.dat"""
write (*,'(a31)') " -Output file = ""fuout.dat"""
write (*,*)
write (*,'(a19)') " Input file format"
write (*,'(a41)') " -Line 1: T(C), P(bars) [free format]"
write (*,'(a39,a35)')
+ " -Additional lines: X(CO2), X(H2O),",
+ " and X(NaCl) triplets [free format]"
write (*,*)
write (*,*)
read(in1,*) tc, pbar
write(iout1,1000) tc, pbar
write(*,1000) tc, pbar
1000 format (1X,F5.0,' DegC',2X,F6.0,' Bars',' ',/)
write (iout1,1100)
write (*,1100)
1100 format (4x, 'XCO2',12x, 'phiCO2',10x, 'fCO2',
+ 12x, 'XH2O',12x, 'phiH2O',10x, 'fH2O',
+ 12x,'XNaCl',11x,'phiNaCl', 9x,'fNaCl',/,
+ 4x,'------',10x,'------',10x,'------',
+ 10x,'------',10x,'------',10x,'------',
+ 10x,'------',10x,'------',10x,'------')
tk = tc + 273.15d0
pat = pbar/1.01325d0
* X(1) = X(CO2), X(2) = X(H2O), X(3) = X(NaCl)
60 read(in1,*,end=999) x(1), x(2), x(3)
call trkmix(tk,pat,x,gam)
write( *,1110) (x(m),gam(m),x(m)*gam(m)*pbar,
+ m = 1,3)
write(iout1,1110) (x(m),gam(m),x(m)*gam(m)*pbar,
+ m = 1,3)
1110 format(9(2x,e14.8))
goto 60
999 stop
end
SUBROUTINE TRKMIX(T,P,Y,GAM)
*======================================================================|
* Calculation of fluid mixture properties for the ternary system
* using De Santis et al. (1974)
* Modifications of the Redlich揖wong equation together with
* A and B parameters for H2O. CO2, and NaCl from Bowers and
* Helgeson (1983).
* P is in atmospheres, T is in deg kelvin, y(I) is the
* mole fraction. Gam(I) is the fucacity coefficient for
* the components in the mixture.
* Written by T. Bowers (1982)
*
* To use subroutine trkmix with program denfind (1) change GAM to
* VOL in subroutine card and (2) leave out lines 0740 through 1310
* ---------------------------------------------------------------------|
real * 8 Y(5), GAM(5)
real * 8 AI(3), BI(3)
real * 8 AZCO2, BCO2, T, P, R, TCEL, R2T, RT, W, BH, BN, PHI
real * 8 AH, AN1, AN2, AN3, BH2O, AZH2O, AH2O, ANACL, AH2OM
real * 8 ACO2M, XK, CO2H2O, BSUM, ASUM, A2BMIX, BMIX, Z
real * 8 CHI1, CHI2, CHI3, CHI5, CHI6, CHI7CW, DW
*======================================================================|
* CO2 parameters from De Santis et al. (1974) and Redlich and Kwong
* (1949)
AZCO2 = 46.d6
BCO2 = 29.7d0
* Define the gas constant in cm3*atm*deg-1*mole-1.
R = 82.05d0
TCEL = T - 273.15d0
R2T = R*R*T**2.5d0
RT = R * T**1.5d0
IF(Y(3) .LE. 1.d-8) GO TO 11
* Calculate weight percent NaCl relative to H2O + NaCl
W = 100.d0*Y(3)*58.448d0/(Y(3)*58.448d0 + Y(2)*18.016d0)
GO TO 12
11 W = 0.d0
* Begin calculating water a and b parameters
12 BH = 14.6d0
BN = - .04420283d0
AH = 4.881243d0 + .1823047d-2*TCEL - .1712269d-4*TCEL*TCEL
+ + .6479419d-8*TCEL*TCEL*TCEL
AN1 = .02636494d0 - .536994d-3*TCEL + .2687074d-5*TCEL*TCEL
+ - .4321741d-8*TCEL*TCEL*TCEL
AN2 = .6802827d-2 - .948023d-4*TCEL + .3770339d-6*TCEL*TCEL
+ - .5075318d-9*TCEL*TCEL*TCEL
AN3 = .5235827d-4 - .3505272d-7*TCEL
BH2O = BH + W*BN
AZH2O = EXP(AH + W*AN1 + W*W*AN2 + W*W*W*AN3)*10d5
AH2O = 111.3057d0 + 50.70033d0*EXP(-.982646d-2*TCEL)
ANACL = -8.05658d0*EXP(-.982646d-2*TCEL)
20 AH2O = AH2O * 10.d5
ANACL = ANACL * 10.d5
AH2OM = AH2O + W*ANACL
* CO2 parameter from Holloway (1977)
* Interaction coefficient from De Santis et al. (1974)
ACO2M = 73.03d0 - .0714*TCEL + 2.157d-5* TCEL*TCEL
ACO2M = ACO2M * 10.d5
XK=EXP(-11.071d0+(5953.d0/T)-(2.746d6/(T*T))+(4.646d8/(T*T*T)))
CO2H2O = XK * 0.5d0 * R2T
CO2H2O = CO2H2O + SQRT(AZCO2 * AZH2O)
Y(5) = Y(2) + Y(3)
BSUM = BCO2*Y(1) + BH2O*Y(5)
ASUM = Y(1)*Y(1)*ACO2M + Y(5)*Y(5)*AH2OM + 2.d0*Y(1)*Y(5)*CO2H2O
A2BMIX = ASUM / (BSUM * RT)
BMIX = BSUM / (R * T)
CALL REDKW(T,P,BMIX,A2BMIX,Z)
* Calculate V
VOL = R*T*Z/P
* For use of this subroutine with program denfind. Lines 0740
* through 1310 need not be included.
* Begin calculation of fugacity coefficients of the components
CHI1 = DLOG(VOL/(VOL-BSUM))
CHI2 = DLOG((VOL+BSUM)/VOL)
CHI3 = BSUM/ (VOL + BSUM)
DW = 58.448d0/18.016d0/100.d0
CW = 18.016d0/58.448d0/100.d0
IF(ABS(W-100.d0).LT.1.d-6) GO TO 16
GO TO 17
16 BI(1) = BCO2
BI(2) = BH2O - CW*BN*W*W*Y(5)/Y(3)
BI(3) = BH2O
AI(1) = 2.d0*Y(1)*ACO2M + 2.d0*Y(5)*CO2H2O
AI(2) = 2.d0*Y(5)*AH2OM + 2.d0*Y(1)*CO2H2O -
+ CW*W*W*Y(5)*(ANACL)
+ *y(5)/y(3) - Y(1)*Y(5)*SQRT(AZCO2*AZH2O)*(AN1 + 2.d0*W*AN2
+ + 3.d0*W*W*AN3)*CW*W*W/y(3)
AI(3) = 2.d0*Y(3)*AH2OM
GO TO 86
17 BI(1) = BCO2
IF(W.LE.1.d-8) GO TO 85
BI(2) = BH2O - CW*BN*W*W*Y(5)/Y(3)
BI(3) = BH2O + DW*BN*(100.d0 - W)*(100.d0 - W)*Y(5)/Y(2)
AI(1) = 2.d0*Y(1)*ACO2M + 2.d0*Y(5)*CO2H2O
AI(2) = 2.d0*Y(5)*AH2OM + 2.d0*Y(1)*CO2H2O -
+ CW*W*W*Y(5)*(ANACL)
+ *Y(5)/Y(3) - y(1)*Y(5)*SQRT(AZCO2*AZH2O)*(AN1 + 2.d0*W*AN2
+ + 3.d0*W*W*AN3)*CW*W*W/Y(3)
AI(3) = 2.d0*Y(5)*AH2OM + 2.d0*Y(1)*CO2H2O +
+ DW*(100.d0 - W)*(ANACL)
+ *(100.d0 - W)*Y(5)*Y(5)/Y(2) + Y(1)*Y(5)*SQRT(AZCO2*AZH2O)
+ *(AN1 + 2.d0*W*AN2 + 3.d0*W*W*AN3)*DW*(100.d0-W)
+ *(100.d0-W)/Y(2)
GO TO 86
85 AI(2) = 2.d0*Y(5)*AH2OM + 2.d0*Y(1)*CO2H2O
AI(1) = 2.d0*Y(1)*ACO2M + 2.d0*Y(5)*CO2H2O
IF(Y(2).LT.1.d-8) GO TO 91
AI(3) = 2.d0*Y(5)*AH2OM + 2.d0*Y(1)*CO2H2O +
+ DW*(100.d0 - W)*(ANACL)
+ *(100.d0 - w)*y(5)*y(5)/y(2) + y(1)*Y(5)*SQRT(AZCO2*AZH2O)
+ *(AN1 + 2.d0*W*AN2 + 3.d0*W*W*AN3)*DW*(100.d0 - W)
+ *(100.d0 - W)/Y(2)
BI(3) = BH2O + DW*BN*(100.d0 - W)*(100.d0 - W)*Y(5)/Y(2)
GO TO 93
91 AI(3) = 2.d0*y(1)*CO2H2O
BI(3) = BH2O
93 BI(2) = BH2O
86 DO 202 J = 1, 3
CHI5 = BI(J)/(VOL-BSUM)
CHI6 = (ASUM*BI(J))/(RT*BSUM**2)
CHI7 = AI(J)/(RT*BSUM)
PHI = CHI1+CHI5-CHI7*CHI2+CHI6*(CHI2-CHI3)-DLOG(Z)
PHI = EXP(PHI)
201 GAM(J) = PHI
202 CONTINUE
1000 CONTINUE
RETURN
END
SUBROUTINE REDKW(T,P,B,A2B,Z)
*======================================================================|
* A subroutine to calculate compressibility factor
* with the Redlich揖wong equation following the method
* of Edmister(1968). T is in kelvins, P is in atmospheres
* written by John R. Holloway(1973)
*======================================================================|
* Hmmmm. Not a single variable declaration. Let's add those...
Real * 8 T, B, P, BP, TH, A2B, RR, QQ, XN, XM, ARG, FP, Z
Real * 8 TANPHI, COSPHI, THETA
Real * 8 FAC, PHI, RH, R1, R2, R3, X, XN2, XMM, XNN
BP = B * P
TH = 1.d0 / 3.d0
IF(A2B.LE.0.0d0) A2B = 0.001d0
RR = -A2B * BP * BP
QQ = BP * (A2B - BP - 1.0d0)
XN = (QQ/3.d0) + RR - (2.d0/27.d0)
XM = QQ - 1.d0 /3.d0
ARG = ((XN*XN) / 4.d0) + (XM*XM*XM/27.d0)
IF(ARG.GT.0.0d0) GO TO 5
IF(ARG.LT.0.0d0) GO TO 15
FP = 1.0d0
Z = 1.0d0
RETURN
15 COSPHI = SQRT((XN*XN/4.d0) / (-XM*XM*XM/27.d0))
IF(XN.GT.0.0d0) COSPHI = -COSPHI
TANPHI = SQRT(1.d0 - COSpHI*COSPHI) / COSPHI
PHI =ATAN(TANPHI)/ 3.0d0
FAC = 2.d0 * SQRT(-XM/3.d0)
R1 = FAC * COS(PHI)
R2 = FAC * COS(PHI + 2.094395103d0)
R3 = FAC * COS(PHI + 4.188790206d0)
PHI = 3.d0 * PHI
* Sort for largest root
RH = R2
IF(R1.GT.R2) RH = R1
IF(R3.GT.RH) RH = R3
* Determine the value of theta
IF (RH.EQ.R1) THETA = 0.0d0
IF (RH.EQ.R2) THETA = 2.094395103d0
IF (RH.EQ.R3) THETA = 4.188790206d0
Z = RH + TH
GO TO 20
5 X = SQRT(ARG)
F = 1.0d0
XN2 = -XN / 2.d0
XMM = XN2 + X
IF(XMM.LT.0.0d0) F = -1.0
XMM = F * ((F * XMM) ** TH)
F = 1.0d0
XNN = XN2 - X
IF(XNN.LT.0.0d0) F = -1.0d0
XNN = F * ((F * XNN) ** TH)
Z = XMM + XNN + TH
20 RETURN
END