c----------------------------------------------------------------------- c FLIB - fluid phase subroutines common to FRENDLY, VERTEX, COHSRK, c RK, and BUILD. c Unless otherwise noted, the subroutines herein were written by c J. A. D. Connolly. c Please do not distribute this source. c----------------------------------------------------------------------- subroutine cfluid (fo2,fs2) c----------------------------------------------------------------------- c subroutine cfluid call fluid equations of state depending on c the value of ifug. The GCOH eos return ln(fo2) as fo2. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n) parameter (h5=7) double precision xf(2) common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst11 /f(2) * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn c----------------------------------------------------------------------- if (xco2.gt.1.d0) then xco2 = 1.d0 else if (xco2.lt.0.) then xco2 = 0. end if xf(2) = xco2 xf(1) = 1.d0 - xco2 goto (1,2,3,4,5,6,7,8,9,10,11,12,13, * 14,15,16,17,18,19,19,21,22),ifug call mrk goto 99 1 call hsmrk goto 99 2 call qrkmrk goto 99 3 call saxfei goto 99 4 call brmrk goto 999 5 call hprk goto 99 6 call trkmrk goto 99 7 call cohgra (fo2) goto 999 8 call cohhyb (fo2) goto 999 9 call cohfit (fo2) goto 999 10 call hocgra (fo2) goto 999 11 call hocmrk (fo2) goto 999 12 call cohsgr (fo2,fs2) goto 999 13 call hh2ork (fo2) goto 999 14 call error (999,xco2,ibuf,'NOGGY') goto 999 15 call lohork (fo2) goto 999 16 call homrk (fo2) goto 999 17 call hosmrk (fo2,fs2) goto 999 18 call dhhsrk goto 999 19 call xoxsrk (fo2,fs2) goto 999 21 call hcrk goto 99 22 call dhcork goto 99 99 do 90 i = 1, 2 if (xf(i).gt.1.d-38) goto 90 f(i) = -9.9d9 90 continue 999 end subroutine rfluid (irk,ifug,lun) c--------------------------------------------------------------------- c irk = 1 - write/read prompt for fluid equations of state c irk = 2 - write fluid equation of state for outtit c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n) parameter (k6=7,l2=5,nrk=22) character rkname(0:nrk)*54, y*1,vname*8, gname*78, xname*18 common / cst41 /n1,n2,n3,n4,n5,n6,n7,n8,n9,io3,io4,io5,io9 * / csta2 /xname(k6),vname(l2),gname(l2) * / cst24 /ipot,jv(l2),iv(l2)/ cst112 /buf(5) * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx save rkname data (rkname(i), i = 0, 11)/ * 'X(CO2) Modified Redlich-Kwong (MRK/DeSantis/Holloway)', * 'X(CO2) Kerrick & Jacobs 1981 (HSMRK)', * 'X(CO2) Hybrid MRK/HSMRK', * 'X(CO2) Saxena & Fei 1987 pseudo-virial expansion', * 'Bottinga & Richet 1981 (CO2 RK)', * 'X(CO2) Holland & Powell 1991, 1995 (CORK)', * 'X(CO2) Hybrid Haar et al 1979/CORK (TRKMRK)', * 'f(O2/CO2)-f(S2) Graphite buffered COHS MRK fluid', * 'f(O2/CO2)-f(S2) Graphite buffered COHS hybrid-EoS fluid', * 'Max X(H2O) GCOH fluid Cesare & Connolly 1993', * 'X(O) GCOH-fluid hybrid-EoS Connolly & Cesare 1993', * 'X(O) GCOH-fluid MRK Connolly & Cesare 1993'/ data (rkname(i), i = 12, nrk)/ * 'X(O)-f(S2) GCOHS-fluid hybrid-EoS Connolly & Cesare 1993', * 'X(H2) H2-H2O hybrid-EoS', * 'EoS Birch & Feeblebop (1993)', * 'X(H2) low T H2-H2O hybrid-EoS', * 'X(O) H-O HSMRK/MRK hybrid-EoS', * 'X(O) H-O-S HSMRK/MRK hybrid-EoS', * 'X(CO2) Delany/HSMRK/MRK hybrid-EoS, for P > 10 kb', * 'X(O)-X(S) COHS hybrid-EoS Connolly & Cesare 1993', * 'X(O)-X(C) COHS hybrid-EoS Connolly & Cesare 1993', * 'X(CO2) Halbach & Chatterjee 1982, P > 10 kb, hybrid-Eos', * 'X(CO2) DHCORK, hybrid-Eos'/ if (irk.eq.2) then write (lun,1060) rkname(ifug) goto 99 end if elag = 0.d0 ibuf = 1 dlnfo2 = 0.d0 10 write (*,1000) do 11 i = 0, nrk 11 write (*,1070) i,rkname(i) read (*,*,iostat=ier) ifug if (ifug.gt.nrk) ier = 1 call rerror (ier,*10) if (ifug.eq.12.or.ifug.eq.17.or.ifug.eq.19.or.ifug.eq.20) then c COHS & HOS equations of state c get sulfur fugacity constraint: vname(3) = 'X(O)' 12 write (*,1090) read (*,*,iostat=ier) ibuf if (ibuf.gt.3.or.ibuf.lt.1) ier = 1 call rerror (ier,*12) if (ibuf.eq.2) then c if ibuf = 2 dlnfo2 is the fe/s c of pyrrhotite 13 write (*,1100) read (*,*,iostat=ier) dlnfo2 call rerror (ier,*13) else if (ibuf.eq.3) then c if ibuf = 2 dlnfo2 is the fs2 14 write (*,1110) read (*,*,iostat=ier) dlnfo2 call rerror (ier,*14) dlnfo2 = 2.302585093d0 * dlnfo2 end if else if (ifug.gt.9.and.ifug.lt.13.or.ifug.gt.13) then if (ifug.ne.18) vname(3) = 'X(O)' else if (ifug.eq.13) then vname(3) = 'X(H2)' else if (ifug.ge.7.and.ifug.le.9) then c chosen COH speciation option c check that XCO2 isn't a independent c variable: if (iv(1).eq.3.or.iv(2).eq.3) then call warn (172,dlnfo2,ifug,'RFLUID') goto 10 end if if (ifug.eq.9) goto 99 c change default buffer 20 write (*,1020) read (*,1030) y ibuf = 2 dlnfo2 = 0.0 if (y.eq.'y'.or.y.eq.'Y') then write (*,1010) read (*,*,iostat=ier) ibuf call rerror (ier,*20) if (ibuf.gt.5.or.ibuf.lt.1) then call warn (173,dlnfo2,ifug,'RFLUID') goto 20 end if c ibuf = 5, define own buffer if (ibuf.eq.5) then 45 write (*,1180) read (*,*,iostat=ier) buf call rerror (ier,*45) goto 30 c ibuf = 3, constant fo2. else if (ibuf.eq.3) then 40 write (*,1140) read (*,*,iostat=ier) dlnfo2 call rerror (ier,*40) dlnfo2 = 2.302585093 * dlnfo2 goto 50 end if c ibuf 2 or 1, permit del(fo2) 30 write (*,1040) read (*,1030) y if (y.eq.'y'.or.y.eq.'Y') then write (*,1050) read (*,*,iostat=ier) dlnfo2 call rerror (ier,*30) dlnfo2 = 2.302585093 * dlnfo2 end if end if end if c for graphite EoS's allow c reduced gph activity: 50 if (ifug.ge.7.and.ifug.le.12.and.ifug.ne.9) then write (*,1170) read (*,1030) y hu = 0 if (y.eq.'y'.or.y.eq.'Y') hu = 1 write (*,1120) read (*,1030) y if (y.eq.'y'.or.y.eq.'Y') then write (*,1130) read (*,*,iostat=ier) elag call rerror (ier,*50) elag = dlog (elag) end if else if (ifug.eq.19) then c for XO-XS EoS get S/C write (*,1150) read (*,*,iostat=ier) elag call rerror (ier,*50) else if (ifug.eq.20) then c for XO-XC EoS get X(C) write (*,1160) read (*,*,iostat=ier) elag call rerror (ier,*50) end if 1000 format (' Select fluid equation of state:') 1010 format (' Select buffer: ',/, * ' 1 - aQFM, Holland & Powell, 298-1200K',/, * ' 2 - Maximum H2O content, 523-1273K, .5-30kbar',/, * ' 3 - constant f(O2)',/, * ' 4 - aQ-Ru-Cc-Tn-Gph, Holland & Powell',/, * ' 5 - ln(f(O2))= a + (b + c*p)/t + d/t**2 + e/t**3 ') 1020 format (' Modify default buffer (max H2O) (Y/N)? ') 1030 format (a1) 1040 format (' Modify calculated fO2 by a constant (Y/N)? ') 1050 format (' Enter constant in units of log10(fO2): ') 1060 format (' Fluid equation of state: ',a53) 1070 format (2x,i2,' - ',a54) 1080 format (/) 1090 format (' Choose a buffer:',/ * ,' 1 - Pyrite + Pyrrhotite', * /,' 2 - Pyrrhotite',/,' 3 - f(S2) ') 1100 format (' Enter atomic Fe/S of pyrrhotite: ') 1110 format (' Enter log10[f(S2)]: ') 1120 format (' Reduce graphite activity (Y/N)? ') 1130 format (' Enter activity of graphite: ') 1140 format (' Enter log10[f(O2)]: ') 1150 format (' Enter X(S), i.e., {n(S)/[n(S) + n(C)]}: ') 1160 format (' Enter X(C), i.e., {n(C)/[n(S)+n(C)+n(O)+n(H)]}: ') 1170 format (' Compute f(H2) & f(O2) as the dependent fugacities', * /,' (do not unless you project through carbon) (Y/N)? ') 1180 format (' Enter a-e :') 99 end subroutine cohsgr (fo2,fs2) c---------------------------------------------------------------------- c program to calculate graphite saturated C-H-O-S speciation as c a function of XO using an MRK/HSMRK hybrid. c Species are CO2 CH4 CO H2 H2O H2S O2 SO2 COS. The latter 3 species c are only significant for reduced graphite activities. c---------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10) integer ins(nsp) common / cst11 /fh2o,fco2/ cst5 /p,t,xo,u1,u2,tr,pr,rj,ps * / cstcoh /xh2o,xco2,xco,xch4,xh2,xh2s,xo2,xso2,xcos,xet, * gh2o,gco2,gco,gch4,gh2,gh2s,go2,gso2,gcos,get,v(nsp) * / cstchx /gmh2o,gmco2,gmch4,vm(3)/ cst26 /vol * / cst100 /rat,elag,gz,gy,gx,ibuf,hu,hv,hw,hx save tol, ins data tol, ins/ .0001d0, 1,2,3,4,5,6,7,8,9,0/ nit = 0 oh2o = 2.d0 xh2 = 0.00001d0 c fs2 = 1/2 ln (fs2), c k's are ln(k) call setfs2 (fs2,kh2s,kso2,kcos) call setup (ghh2o,ghco2,ghch4,kh2o,kco2,kco,kch4) c3 = dexp (kch4) * p c1 = dexp (kco2 - 2.d0*kco) * p c2 = dexp (kh2o - kco) * p c4 = dexp (kh2s + fs2) c5 = dexp (kcos + fs2) c6 = dexp (kso2 - 2.d0*kco + fs2) * p c7 = dexp (-2.d0*kco) * p c outer iteration loop: 10 k1 = c1 * gco**2/gco2 k2 = c2 * gco * gh2/gh2o k3 = c3 * gh2**2/gch4 k4 = c4 * gh2/gh2s k5 = c5 * gco/gcos k6 = c6 * gco**2/gso2 k7 = c7 * gco**2/go2 c solve for xh2, xco 15 call evlxh1 (k1,k2,k3,k4,k5,k6,k7,xo,xh2,xco,ier) if (ier.ne.0) call warn (501,xo,ier,'COHSGR') xch4 = k3 * xh2**2 xh2o = k2 * xh2 * xco xco2 = k1 * xco**2 xh2s = k4 * xh2 xso2 = k6 * xco**2 xcos = k5 * xco xo2 = k7 * xco**2 nit = nit + 1 if (nit.gt.100) call warn (502,xo,ier,'COHSGR') if (dabs(xh2o-oh2o).lt.tol*xh2o) goto 90 oh2o = xh2o call mrkmix (ins, 9) gh2o = ghh2o * gh2o gco2 = ghco2 * gco2 gch4 = ghch4 * gch4 goto 10 90 vol = vol + xh2o * vm(1) * + xco2 * vm(2) * + xch4 * vm(3) goto (91), hu fh2o = dlog(gh2o*p*xh2o) fco2 = dlog(gco2*p*xco2) fo2 = 2.d0 * (dlog(gco*p*xco) - kco) goto 99 91 fh2o = dlog(gh2*p*xh2) fco2 = 2.d0 * (dlog(gco*p*xco) - kco) 99 end subroutine evlxh1 (k1,k2,k3,k4,k5,k6,k7,xo,xh2,xco,ier) implicit double precision (a-h,k-z) it = 0 ier = 0 f0 = 2.d0*(k7 + k6 + k1) e0 = 2.d0*xo e1 = 1/f0 e2 = 1.d0 + k5**2 + 2.d0*(k5 + f0) e3 = 2.d0*k2*(1.d0 + k5) - 2.d0*f0*(k4 + 1.d0) e4 = k2**2 - 2.d0*k3*f0 e5 = e0 + e0*k4 e6 = 4.d0*xo*k3 e7 = xo - k5 - 1.d0 + xo*k5 e8 = f0 * (xo - 1.d0) e9 = k2*(3.d0* xo - 1.d0) 10 r0 = xh2 t2 = xh2**2 t10 = e2 + e3*xh2 + e4*t2 if (t10.lt.0.d0) then c if t10 < 0, bad guess c for xh2, find roots: c1 = dsqrt (e3**2 - 4.d0*e4*e2) xh2 = 0.9*(-c1 - e3/2.d0/ e4) r0 = xh2 t2 = xh2**2 t10 = e2 + e3*xh2 + e4*t2 end if t10 = dsqrt (t10) t11 = t10 - 1.d0 - k2*xh2 - k5 xco = e1*t11 g = e5*xh2 + e6*t2 + (e7 + e8*xco + e9*xh2)*xco t15 = (e3+2.d0*e4*xh2)/2.d0/t10 - k2 dg = e5 + 2.d0*e6*xh2 + e1*t15*(e9*xh2 + e7) * + t11*(2.d0*e8*e1**2*t15 + e9*e1) xh2 = r0 - g/dg if (xh2.lt.0.d0) xh2 = r0/2.d0 c converged: if (dabs((xh2-r0)/xh2).lt.0.1d-6) goto 999 it = it + 1 if (it.gt.1000) then ier = 2 goto 99 end if goto 10 999 xco = e1*(dsqrt(e2 + (e3 + e4*xh2)*xh2) - 1.d0 - k2*xh2 - k5) 99 end subroutine setup (ghh2o,ghco2,ghch4,kh2o,kco2,kco,kch4) c---------------------------------------------------------------------- c program to setup C-O-H-S speciation calculations c---------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10,h5=7) integer ins(nsp), jns(3) common / cst5 /p,t,xo,u1,u2,tr,pr,rj,ps * / cstcoh /xh2o,xco2,xco,xch4,xh2,xh2s,xo2,xso2,xcos,xet, * gh2o,gco2,gco,gch4,gh2,gh2s,go2,gso2,gcos,get,v(nsp) * / cstchx /gmh2o,gmco2,gmch4,vm(3) * / cst100 /rat,xs,ag,gy,gx,ibuf,hu,hv,hw,hx * / cst10 /is(2),i3(h5),ifug,i1,i2 save jns, ins data ins, jns/ 1,2,3,4,5,6,7,8,9,0,1,2,3/ t2 = t * t t3 = t2 * t c check if xo is <1, >0, c reset if necessary if (xo.gt.0.9999999999d0) then xo = 0.9999999999d0 else if (xo.lt.0.0000000001d0) then xo = 0.0000000001d0 end if c get pure species fugacities call hsmrkp (ins, 9, jns, 3) ghh2o = gh2o/gmh2o ghco2 = gco2/gmco2 ghch4 = gch4/gmch4 c graphite pressure effect: dg = p*( 1.8042d-06 + (0.058345d0 - 8.42d-08*p)/t ) c graphite activity effect: if (ifug.ne.20) dg = dg + xs c ln k's fitted from robie: kh2o = -7.028214449d0 + 30607.34044d0/t - 475034.4632d0/t2 * + 50879842.55d0/t3 kco2 = .04078341613d0 + 47681.676177d0/t - 134662.1904d0/t2 * + 17015794.31d0/t3 + dg kco = 10.32730663d0 + 14062.7396777d0/t- 371237.1571d0/t2 * + 53515365.95d0/t3 + dg kch4 = -13.86241656d0 + 12309.03706d0/t - 879314.7005d0/t2 * + .7754138439d8/t3 + dg end subroutine xoxsrk (fo2,fs2) c---------------------------------------------------------------------- c program to calculate C-H-O-S speciation as a function of XO and XS c or XC using an MRK/HSMRK hybrid. c---------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10,h5=7) integer ins(nsp) common / cst11 /fh2o,fco2/ cst5 /p,t,xo,u1,u2,tr,pr,rj,ps * / cstcoh /xh2o,xco2,xco,xch4,xh2,xh2s,xo2,xso2,xcos,xet, * gh2o,gco2,gco,gch4,gh2,gh2s,go2,gso2,gcos,get,v(nsp) * / cstchx /gmh2o,gmco2,gmch4,vm(3)/ cst26 /vol * / cst100 /rat,xs,ag,gy,gx,ibuf,hu,hv,hw,hx * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn save tol, ins data tol, ins/ .0001, 1,2,3,4,5,6,8,9,2*0/ nit = 0 oh2o = 2.d0 xh2 = 0.00001d0 xh2o = 0.1d0 c this fs2 = 1/2 ln (fs2), c k's are ln(k) call setfs2 (fs2,kh2s,kso2,kcos) call setup (ghh2o,ghco2,ghch4,kh2o,kco2,kco,kch4) c c1 = dexp (kcos + fs2) c2 = dexp (kco2 - kco - kh2o) c3 = dexp (kh2s + fs2) c5 = dexp (kso2 + 4.d0*(kco-kco2) + 2.d0*kh2o + fs2)/p c6 = dexp (kch4 + kh2o - kco) * p * p c outer iteration loop: 10 k1 = c1 * gco/gcos k2 = c2 * gco * gh2o/gh2/gco2 k3 = c3 * gh2/gh2s k5 = c5 * gco2**4 * gh2**2/gco**4/gso2/gh2o/gh2o k6 = c6 * gh2**3 * gco/gh2o/gch4 c solve for xh2, xco 15 if (ifug.eq.19) then call evlxh2 (k1,k2,k3,k5,k6,xo,xs,xh2,xco,xh2o,ier) else call evlxh3 (k1,k2,k3,k5,k6,xo,xs,xh2,xco,xh2o,ier) end if if (ier.ne.0) call warn (501,xh2,ier,'XOXSRK') xch4 = k6 * xco * xh2**3/xh2o xco2 = k2 * xh2o * xco/xh2 xh2s = k3 * xh2 xso2 = k5 * xh2o**2/xh2**2 xcos = k1 * xco nit = nit + 1 if (nit.gt.100) then call warn (175,xh2,ier,'XOXSRK') goto 90 end if if (dabs(xh2o-oh2o).lt.tol*xh2o) goto 90 oh2o = xh2o call mrkmix (ins, 8) gh2o = ghh2o * gh2o gco2 = ghco2 * gco2 gch4 = ghch4 * gch4 goto 10 90 fh2o = dlog(gh2o*p*xh2o) fco2 = dlog(gco2*p*xco2) fo2 = 2.d0 * (fh2o - dlog(gh2*p*xh2) - kh2o) c compute graphite activity: ag = fco2 - fo2 - kco2 vol = vol + xh2o * vm(1) * + xco2 * vm(2) * + xch4 * vm(3) 99 end subroutine evlxh2 (c1,c2,c3,c5,c6,xo,xs,xh2,xco,xh2o,ier) implicit double precision (a-h,k-z) jt = 0 ier = 0 c7 = 2.d0*c5 d2 = c5 - xs * c5 d3 = c3 - xs * c3 d4 = c1 - 2.d0*xs*c1 - xs d6 = 2.d0 * c2 d7 = xs * c2 d8 = xs * c6 d9 = 4.d0 * c6 d10 = 2.d0 * c3 d11 = 12.d0 * c6 d12 = 5.d0 * d8 d25 = d10 + 2.d0 f10 = 3.d0*d2 100 r1 = xh2o it = 0 ier = 0 d1 = xh2o**2 d5 = c7 * d1 d13 = d7 * d1 d14 = 2.d0*d4*xh2o d15 = d2*d1 d16 = d3*xh2o d17 = d4*xh2o d18 = -3.d0* d16 d19 = -4.d0*c5*d1 d20 = 3.d0*xh2o d21 = d9/xh2o d22 = d11*d3 d23 = d11/xh2o d24 = d15*xh2o g1 = -6.d0*c2*d1*d3 g2 = d6*xh2o c solve g for xh2, guessed xh2o 10 r0 = xh2 c evaluate g: e1 = xh2**2 e2 = e1*xh2 e3 = e2*xh2 e4 = e3*xh2 t14 = d24 + d16*e2 t31 = d17*e1 - xh2*d13 - d8*e4 t39 = -t14/t31 t37 = g2*t39/xh2 t43 = d5/e1 t45 = c1*t39 g = (t37 + t39 + t43 + xh2o + t45)/ * (t37 + t39 + t43 + d20 + t45 + 2.d0*xh2 * - d21*t14/t31*e2 + d10*xh2) - xo c evaluate dg: t27 = g1*xh2/t31 t37 = t31**2 t49 = d14*xh2 - d13 - d12*e3 t55 = g2*t14/t37/xh2*t49 t61 = -g2*t39/e1 t65 = d18*e1/t31 t67 = t14/t37*t49 t71 = d19/e2 t73 = c1*t65 t74 = c1*t67 t80 = -g2*t14/t31/xh2 t88 = c1*t39 t98 = t80 + t39 + t43 + d20 + t88 + 2.d0*xh2 * + d21*t39*e2 + d10*xh2 dg = (t27 + t55 + t61 + t65 + t67 + t71 + t73 + t74)/t98 * - (t80 + t39 + t43 + xh2o + t88)/t98/t98 * *(t27 + t55 + t61 + t65 + t67 + t71 + t73 + t74 * - d22*e4/t31 + d21*e2*t67 + d23*t39*e1 + d25) xh2 = r0 - g/dg if (xh2.lt.0.) xh2 = r0/2.d0 c converged: if (dabs((xh2-r0)/xh2).lt.0.1d-6) goto 40 it = it + 1 if (it.gt.100) then ier = 2 goto 40 end if goto 10 40 it = 0 c use xh2 to refine guesssed xh2o 20 t10 = xh2**2 t11 = xh2**3 t25 = xh2**5 f2 = 2.d0*d7*xh2 f3 = c2/xh2 f4 = xh2 - 1.d0 + c3*xh2 f7 = c2*xh2 f8 = d3*t11 f9 = d7*xh2 f11 = c5/t10 f12 = d4*t10 30 r0 = xh2o c evaluate f: t4 = xh2o**2 f13 = d2*t4 f1 = f12*xh2o - t4*f9 - d8*t25 t14 = f13*xh2o + xh2o*f8 t32 = t14/f1 f5 = c6*t11 f6 = f5/xh2o f = -t32 - f3*xh2o*t32 - f5*t32/xh2o - c1*t32 * + t4*f11 + xh2o + f4 c evaluate df: t13 = f10*t4 + f8 t31 = t13/f1 t37 = d2*t4*xh2o + xh2o*f8 t38 = f1**2 t45 = d4*t10 - f2*xh2o t47 = t37/t38*t45 t49 = -f3*f1 df = t47 - t31 + t37*t49 + xh2o*t13*t49 * + f7*xh2o*t47 - f6*t31 * + f6*t47 + f5*t37/f1/t4 * - c1*t31 + c1*t47 + c7*xh2o/t10 + 1.d0 xh2o = r0 - f/df if (xh2o.lt.0.d0) then xh2o = r0/2.d0 else if (xh2o.ge.1.d0) then xh2o = r0 + (1.d0 - r0)/2.d0 end if c converged: if (dabs((xh2o-r0)/xh2o).lt.0.1d-6) goto 99 it = it + 1 if (it.gt.100) then ier = 2 goto 99 end if goto 30 99 t4 = xh2o**2 xco = -(d2*t4*xh2o + d3*xh2o*t11)/ * (d4*t10*xh2o - d7*xh2*t4 - d8*t25) c is new xh2o same as guess?: if (dabs((xh2o-r1)/xh2o).lt.0.1d-6) goto 9999 jt = jt + 1 if (jt.gt.100) then ier = 2 goto 9999 end if goto 100 9999 end subroutine evlxh3 (c1,c2,c3,c5,c6,xo,xc,xh2,xco,xh2o,ier) implicit double precision (a-h,k-z) jt = 0 ier = 0 g0 = 3.d0*c5 g1 = 2.d0*c5 g2 = xo*g1 g3 = 3.d0*xo g4 = 2.d0*xo g5 = g4*c3 g6 = 2.d0*c2 g7 = g4*c2 g8 = xo*c1 g9 = 4.d0*xo*c6 g10 = 6.d0*xo g11 = 3.d0*g4 g12 = 3.d0*g5 g13 = 6.d0*c5 c11 = 3.d0*c3 c15 = c11 + 2.d0 c16 = 2.d0*g3 c17 = 3.d0*g2 c18 = 3.d0*g1 c19 = 2.d0*g6 c20 = 2.d0*g7 c25 = 2.d0*c1 c26 = 2.d0*g8 c27 = 5.d0*g9 100 r1 = xh2o t1 = xh2o**2 t2 = t1*xh2o c7 = g1*t2 c8 = g2*t2 c9 = g0*t1 c10 = 3.d0*xh2o c12 = c2*xh2o c14 = g13*t1 c21 = c6/xh2o c22 = g7*t1 c23 = g6*t1 it = 0 ier = 0 c solve c for xh2, guessed xh2o 10 r0 = xh2 c evaluate c: t4 = xh2**2 t5 = t4*t1 t21 = t4**2 t10 = t4*xh2 c13 = t10/xh2o t11 = t10*xh2o t16 = c7 + t5 - c8 - g3*t5 - g4*t11 - g5*t11 t17 = t1*xh2 t19 = t4*xh2o t32 = (g6-g7)*t17+t19+(c1-xo-g8)*t19-g9*t21*xh2 t33 = t16/t32 t39 = c12*t33/xh2 t45 = c21*t33*t10 t47 = c1*t33 c = (-t33 - t39 - t45 - t47)/ * (-3.d0*t39 - 2.d0*t33 + c9/t4 + c10 - 3.d0*t47 + 2.d0*xh2 * -5.d0*t45 + c11*xh2) - xc c evaluate dc: t12 = 2.d0*t17 - g10*t17 - g11*t19 - g12*t19 t27 = t12/t32 t42 = c7 + t5 - c8 -g3*t5 - g4*t11 - g5*t11 t43 = t32**2 t101 = t42/t32 t46 = xh2*xh2o t58 = c23 + 2.d0*t46 + c25*t46 -c22 - g4*t46 - c26*t46 - c27*t21 t60 = t42/t43*t58 t62 = t32*xh2 t65 = c12*t12/t62 t71 = c12/xh2*t60 t76 = c12*t101/t4 t79 = c6*c13/t32 t81 = t12*t79 t87 = c6*c13*t60 t91 = c21*t101*t4 t93 = c1*t27 t95 = c1*t60 t99 = c12*t42/t62 t107 = c1*t101 t111 = t42*t79 t115 = -3.d0*(t99 + t107) + 2.d0*(xh2 - t101) + c9/t4 + c10 * -5.d0*t111 + c11*xh2 dc = (-t27+t60-t65+t71+t76-t81+t87-3*t91-t93+t95) * /t115 - (-t101 - t99 - t111 - t107)/t115**2 * *(+ 3.d0*(t71 - t65 + t76 - t93 + t95) + 2.d0*(t60 - t27) * -c14/t10 + 5.d0*(t87- t81) -15.d0*t91 + c15) xh2 = r0 - c/dc if (xh2.lt.0.) xh2 = r0/2.d0 c converged: if (dabs((xh2-r0)/xh2).lt.0.1d-6) goto 40 it = it + 1 if (it.gt.100) then ier = 2 goto 40 end if goto 10 40 it = 0 c use xh2 to refine guesssed xh2o 20 t4 = xh2**2 t10 = t4*xh2 t27 = t4**2 f1 = c5/t4 f2 = c3*xh2 + xh2 -1.d0 f3 = g9*xh2 f4 = xo*t4 f5 = c2/xh2 f6 = + t4 + c1*t4 - f4 - g8*t4 f7 = g1/t4 f8 = c6*t10 30 r0 = xh2o c evaluate f: t1 = xh2o**2 t2 = t1*xh2o t17 = t1*xh2 t5 = t4*t1 t11 = t10*xh2o t19 = t4*xh2o t7 = g1*t2 + t5 - g2*t2 - g3*t5 - g4*t11 - g5*t11 t3 = g6*t17 + t19 + c1*t19 - g7*t17 - xo*t19 - g8*t19 - t27*f3 t33 = t7/t3 f = -t33 - f5*xh2o*t33 - f8*t33/xh2o - c1*t33 * + f1*t1 + xh2o + f2 c evaluate df: t16 = c18*t1 + 2.d0*t19- c17*t1- c16*t19 - g4*t10 - g5*t10 t32 = t16/t3 t49 = xh2*xh2o t59 = c19*t49 - c20*t49 - f6 t61 = t7/t3**2*t59 t63 = c2*t3/xh2 df = -t32 + t61 - t7*t63 - xh2o*t16*t63 + f5*xh2o*t61 * - f8*t32/xh2o + f8/xh2o*t61 + f8*t7/t3/t1 - c1*t32 * + c1*t61 + f7*xh2o + 1.d0 xh2o = r0 - f/df if (xh2o.lt.0.d0) then xh2o = r0/2.d0 else if (xh2o.ge.1.d0) then xh2o = r0 + (1.d0 - r0)/2.d0 end if c converged: if (dabs((xh2o-r0)/xh2o).lt.0.1d-6) goto 99 it = it + 1 if (it.gt.100) then ier = 2 goto 99 end if goto 30 99 xco = -((g1-g2)*xh2o**2+((1.d0-g3)*xh2o-(g4+g5)*xh2)*xh2**2)*xh2o/ * ((((g6-g7)*xh2o+(1.d0+c1-xo-g8)*xh2)*xh2o-g9*xh2**4)*xh2) c is new xh2o same as guess?: if (dabs((xh2o-r1)/xh2o).lt.0.1d-6) goto 9999 jt = jt + 1 if (it.gt.100) then ier = 2 goto 9999 end if goto 100 9999 end subroutine hh2ork (fo2) c---------------------------------------------------------------------- c program to calculate fh2o, fo2 for H2-H2O mixtures using hsmrk/mrk c hybrid EoS c---------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10) integer ins(nsp), jns(3) common / cst11 /fh2o,fh2/ cst5 /p,t,xv,u1,u2,tr,pr,rj,ps * / cstcoh /xh2o,xco2,xco,xch4,xh2,x(5), * gh2o,gco2,gco,gch4,gh2,g(5),v(nsp) * / cstchx /gmh2o,gmco2,gmch4,vm(3)/ cst26 /vol save ins, jns data ins, jns/1, 5, 8*0, 1, 2*0/ xh2 = xv xh2o = 1.d0 - xh2 c check if xo2 is <1, >0, c reset if necessary. if (xh2.gt.0.9999999999d0) then xh2 = 0.9999999999d0 else if (xh2.lt.0.0000000001d0) then xh2 = 0.0000000001d0 end if c get pure species fugacities call hsmrkp (ins, 2, jns, 1) ghh2o = gh2o/gmh2o c get mrk fugacities: call mrkmix (ins, 2) c evaluate lnk's kh2o = -7.028214449d0 + 30607.34044d0/t - 475034.4632d0/t/t * + 50879842.55d0/t/t/t c gh2o = ghh2o * gh2o fh2o = dlog(gh2o*p*xh2o) fh2 = dlog(gh2*p*xh2) fo2 = 2.d0 * (fh2o - fh2 - kh2o) vol = vol + xh2o * vm(1) end subroutine cohfit (fo2) c---------------------------------------------------------------------- implicit double precision (a-z) c subroutine which returns ln(f2o,fco2) from functions fit to the c values obtained from cohhyp for fo2 => xh2o max. c the functions are of the form: c ln(f) = a*t + b*p + c/t/t + d*p*t + e*p*p + f*t*t + g*sqrt(p*t) c + h*t**3 + i*p**3 + j*p/t/t + k ln t + l ln p c + m/p**2 + n*p*p/t + o*p/t + pp*t/p + q*t*t/p c + rp*t*ln(t) + s*p*ln(t) + tp common / cst11 /fh2o,fco2/ cst5 /p,t,xc,u1,u2,tr,pr,rj,ps save ah,bh,ch,dh,eh,fh,gh,hh,ih,jh,kh,lh,mh,nh,oh,ph,qh,rh,sh,th save ac,bc,cc,dc,ec,fc,gc,hc,ic,jc,kc,lc,mc,nc,oc,pc,qc,rc,sc,tc save ao,bo,co,do,eo,fo,go,ho,io,jo,ko,lo,mo,no,oo,po,qo,ro,so,to data th,ah,bh,ch,dh,eh,fh,gh,hh,ih,jh,kh,lh,mh,nh,oh,ph,qh,rh,sh/ *-121.5649 , -.5498738d-01, -.8052569d-02, -169883.7 , *-.5150895d-06, 0.2909217d-08, 0.1484028d-04, 0.1595388d-02, *-.2286740d-08, -.4611758d-14, -118.8630 , 24.90956 , *-1.283717 , 126339.8 , -.1875289d-05, 1.141895 , *-1.661033 , 0.7550874d-03, 0.9051986d-03, 0.1093571d-02/ data tc,ac,bc,cc,dc,ec,fc,gc,hc,ic,jc,kc,lc,mc,nc,oc,pc,qc,rc,sc/ *-68.24622 , -.2797826d-01, -.5658539d-02, -221752.4 , *-.2227993d-06, -.4785067d-08, -.2949820d-05, -.3942711d-02, *0.8136084d-09, 0.6607593d-13, -126.2944 , 12.42835 , *-.1328584 , -168530.0 , -.1849930d-05, 1.058393 , * 2.131351 , -.9849674d-03, 0.3118428d-02, 0.8233771d-03/ data to,ao,bo,co,do,eo,fo,go,ho,io,jo,ko,lo,mo,no,oo,po,qo,ro,so/ *-804.2316 , -.1652445 , -.5376252d-02, -4037433.d0 , *-.2091203d-06, -.4638105d-08, 0.3753368d-04, -.3853404d-02, *-.5442896d-08, 0.6484263d-13, -121.6754 , 127.5998 , *-.1486220 , -164866.6 , -.1863209d-05, 0.9622612 , * 2.097447 , -.9838123d-03, 0.3077560d-02, 0.7829503d-03/ lp = dlog(p) lt = dlog(t) spt = dsqrt(p*t) p2 = p*p t2 = t*t fh2o = th + t*(ah + dh*p + t*(fh + hh*t) + (ph+qh*t)/p + * rh*lp) * + p*(bh + p*(eh + ih*p) + sh*lt) * + p/t*(jh/t + nh*p + oh) * + kh*lt + lh*lp + ch/t2 + gh*spt + mh/p2 fco2 = tc + t*(ac + dc*p + t*(fc + hc*t) + (pc+qc*t)/p + * rc*lp) * + p*(bc + p*(ec + ic*p) + sc*lt) * + p/t*(jc/t + nc*p + oc) * + kc*lt + lc*lp + cc/t2 + gc*spt + mc/p2 fo2 = to + t*(ao + do*p + t*(fo + ho*t) + (po+qo*t)/p + * ro*lp) * + p*(bo + p*(eo + io*p) + so*lt) * + p/t*(jo/t + no*p + oo) * + ko*lt + lo*lp + co/t2 + go*spt + mo/p2 end subroutine hocgra (fo2) c---------------------------------------------------------------------- c program to calculate GCOH fluid properties as a function of XO c using HSMRK/MRK hybrid see Connolly and Cesare (1992) for details. c---------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10) integer ins(nsp) common / cst11 /fh2o,fco2/ cst5 /p,t,xo2,u1,u2,tr,pr,rj,ps * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx * / cstcoh /xh2o,xco2,xco,xch4,xh2,x(5), * gh2o,gco2,gco,gch4,gh2,g(5),v(nsp) * / cstchx /gmh2o,gmco2,gmch4,vm(3)/ cst26 /vol save tol, ins data tol, ins/ .0001, 1,2,3,4,5,5*0/ nit = 0 oh2o = 2. call setup (ghh2o,ghco2,ghch4,kh2o,kco2,kco,kch4) xt = 1.d0 - xo2 k30 = dexp (kch4) * p k10 = dexp (kco2 - 2.d0*kco) * p k20 = dexp (kh2o - kco) * p c solve 10 k1 = k10 * gco**2/gco2 k2 = k20 * gco * gh2/gh2o k3 = k30 * gh2**2/gch4 c solve for xh2 call evalxh (k1,k2,k3,xt,xh2,ier) c solve for xco call evalxc (k1,k2,k3,xt,xh2,xco) xch4 = k3 * xh2**2 xh2o = k2 * xh2 * xco xco2 = k1 * xco**2 nit = nit + 1 if (nit.gt.100) then write (*,*) ' not converging in hocgra ',xh2o,oh2o goto 99 end if if (dabs(xh2o-oh2o).lt.tol*xh2o) goto 90 oh2o = xh2o call mrkmix (ins, 5) gh2o = ghh2o * gh2o gco2 = ghco2 * gco2 gch4 = ghch4 * gch4 goto 10 90 goto (91), hu fh2o = dlog(gh2o*p*xh2o) fco2 = dlog(gco2*p*xco2) fo2 = 2.d0*(dlog(gco*p*xco) - kco) goto 99 91 fh2o = dlog(gh2*p*xh2) fco2 = 2.d0*(dlog(gco*p*xco) - kco) 99 vol = vol + xh2o * vm(1) * + xco2 * vm(2) * + xch4 * vm(3) end subroutine nogcoh (fo2) c---------------------------------------------------------------------- c program to calculate COH fluid properties for a specified speciation c using HSMRK/MRK hybrid see Connolly and Cesare (1992) for details. c---------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10) integer ins(nsp) common / cst11 /fh2o,fco2/ cst5 /p,t,xo2,u1,u2,tr,pr,rj,ps * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx * / cstcoh /xh2o,xco2,xco,xch4,xh2,x(5), * gh2o,gco2,gco,gch4,gh2,g(5),v(nsp) * / cstchx /gmh2o,gmco2,gmch4,vm(3)/ cst26 /vol save ins data ins/ 1,2,3,4,5,5*0/ call setup (ghh2o,ghco2,ghch4,kh2o,kco2,kco,kch4) call mrkmix (ins, 5) gh2o = ghh2o * gh2o gco2 = ghco2 * gco2 gch4 = ghch4 * gch4 90 goto (91), hu fh2o = dlog(gh2o*p*xh2o) fco2 = dlog(gco2*p*xco2) fo2 = 2.d0*(dlog(gco*p*xco) - kco) goto 99 91 fh2o = dlog(gh2*p*xh2) fco2 = 2.d0*(dlog(gco*p*xco) - kco) 99 vol = vol + xh2o * vm(1) * + xco2 * vm(2) * + xch4 * vm(3) end subroutine hocmrk (fo2) c---------------------------------------------------------------------- c program to calculate GCOH fluid properties as a function of XO using c MRK see Connolly and Cesare (1991) for details. c---------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10) integer ins(nsp) common / cst11 /fh2o,fco2/ cst5 /p,t,xo2,u1,u2,tr,pr,rj,ps * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx * / cstcoh /xh2o,xco2,xco,xch4,xh2,x(5), * gh2o,gco2,gco,gch4,gh2,g(5),v(nsp) * / cst26 /vol save tol, ins data tol, ins/0.0001, 1, 2, 3, 4, 5, 5*0/ t2 = t*t t3 = t2 * t nit = 0 oh2o = 2.d0 c check if xo2 is <1, >0, c reset if necessary. xt = 1.d0 - xo2 if (xt.gt.0.9999999999d0) then xt = 0.9999999999d0 else if (xt.lt.0.0000000001d0) then xt = 0.0000000001d0 end if c graphite pressure effect: dg = p*( 1.8042d-06 + (0.058345 - 8.42d-08*p)/t ) + elag c get pure species fugacities call mrkpur (ins, 5) c ln k's fitted from robie: kh2o = -7.028214449 + 30607.34044/t - 475034.4632/t2 * + 50879842.55/t3 kco2 = .04078341613 + 47681.676177/t - 134662.1904/t2 * + 17015794.31/t3 + dg kco = 10.32730663 + 14062.7396777/t - 371237.1571/t2 * + 53515365.95/t3 + dg c kch4 * p k30 = dexp (-13.86241656d0 + 12309.03706d0/t - 879314.7005d0/t2 * + .7754138439d8/t3 + dg) * p k10 = dexp (kco2 - 2.d0*kco) * p k20 = dexp (kh2o - kco) * p c solve 10 k1 = k10 * gco**2/gco2 k2 = k20 * gco * gh2/gh2o k3 = k30 * gh2**2/gch4 c solve for xh2 call evalxh (k1,k2,k3,xt,xh2,ier) c solve for xco call evalxc (k1,k2,k3,xt,xh2,xco) xch4 = k3 * xh2**2 xh2o = k2 * xh2 * xco xco2 = k1 * xco**2 nit = nit + 1 if (nit.gt.100) then write (*,*) ' not converging in hocmrk ',xh2o,oh2o stop end if if (dabs(xh2o-oh2o).lt.tol*xh2o) goto 90 oh2o = xh2o call mrkmix (ins, 5) goto 10 90 goto (91), hu fh2o = dlog(gh2o*p*xh2o) fco2 = dlog(gco2*p*xco2) fo2 = 2.d0*(dlog(gco*p*xco) - kco) goto 99 91 fh2o = dlog(gh2*p*xh2) fco2 = 2.d0*(dlog(gco*p*xco) - kco) 99 end subroutine evalxh (k1,k2,k3,xt,xh,ier) implicit double precision (a-h,k-z) sign = 1.d0 20 xh = 0.01d0 it = 0 ier = 0 10 r0 = xh call evalg (k1,k2,k3,xt,xh,g,dg,sign) xh = r0 - g/dg if (dabs((xh-r0)/xh).lt.0.1d-6) then c converged: if (xh.gt.0.d0) goto 999 c xh < 0. if (sign.lt.0.) then write (*,*) ' evalxh failed' stop else write (*,*) ' trying 2nd root in evalxh' sign = -1.d0 goto 20 end if end if it = it + 1 if (it.gt.100) then if (sign.gt.0.) then write (*,*) ' evalxh, did not converge on root 1' sign = -1.d0 goto 20 else write (*,*) ' evalxh, neither root converged' stop end if end if goto 10 999 end subroutine evalg (k1,k2,k3,xt,xh,g,dg,sign) implicit double precision (a-z) u1 = xt*k1 k2x = k2*xh xk2x = xt*k2x k22 = k2**2 t8 = xh**2 t9 = k22*t8 t15 = xt**2 t20 = k3*t8 t21 = k1*t20 t24 = k1*xh t32 = dsqrt (4.d0*(t9 - xk2x) * + t15*(6.d0*k2x + 9.d0*t9 + 1.d0 - 32.d0*t21 - 16.d0*t24) * + xt*(32.d0*t21 + 16.d0*t24 - 12.d0*t9)) t34 = - 2.d0*k2x + 3.d0*xk2x + xt - sign * t32 t36 = t34/u1 t40 = xt*k2 t41 = t34**2 t43 = k22*xh t53 = k3*xh t54 = k1*t53 t66 = - 2.d0*k2 + 3.d0*t40 - (4.d0*t43 * + xt*(32.d0*t54-12.d0*t43) - 2.d0*t40 * + t15*(9.d0*t43 + 3.d0*k2 - 32.d0*t54 - 8.d0*k1) * + 8.d0*u1 )/t32 t68 = t66/u1 g = -k2x*t36/4.d0 + t41/k1/t15/16.d0 + t20 - t36/4.d0 + xh - 1.d0 dg = (-k2/u1*t34 - k2x*t68 + t34/k1/t15*t66/2.d0)/4.d0 * + 2.d0*t53 - t68/4.d0 + 1.d0 end subroutine evalxc (k1,k2,k3,xt,xh,xco) implicit double precision (a-z) k2x = k2*xh xk2x = xt*k2x k22 = k2**2 t8 = xh**2 t9 = k22*t8 t15 = xt**2 t21 = k1*k3*t8 t24 = k1*xh t32 = dsqrt (4.d0*(t9 - xk2x) * + t15*(9.d0*t9 + 6.d0*k2x + 1.d0 - 32.d0*t21 - 16.d0*t24) * + xt*(32.d0*t21 - 12.d0*t9 + 16.d0*t24)) xco = (2.d0*k2x - 3.d0*xk2x - xt + t32)/xt/k1/4.d0 if (xco.le.0.0) then write (*,*) ' changed roots in evalxc ' xco = -1.d0/xt/k1*(- 2.d0*k2x + 3.d0*xk2x + xt + t32)/4.d0 if (xco.le.0.0) then write (*,*) ' both roots negative in evalxc' stop end if end if end subroutine fo2buf (fo2) c---------------------------------------------------------------------- c this routine returns the fo2 of buffers as a function of P-T. c corrections for graphite activity and delta fo2 are made from c cst100. c ibuf - 1 - fit to aQFM from Holland and Powell in the range c 298-1200 K, probably not valid above aQ/bQ trans. c should be refit to revised hp data. c ibuf - 2 - fit of fo2 to obtain maximum H2O content. c ibuf - 3 - fo2 input by user c ibuf - 4 - fit of rutile-titantite-a-quartz-calcite-graphite c holland & powell 1993. c expansivity and compressibility ignored. c ibuf - 5 - user buffer function c---------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) common / cst5 /p,t,xo,uc,u2,tr,pr,r,ps * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx * / cst112 /ab,bb,cb,db,eb save ao,bo,co,do,eo,fo,go,ho,io,jo,ko,lo,mo,no,oo,po,qo,ro,so,to data to,ao,bo,co,do,eo,fo,go,ho,io,jo,ko,lo,mo,no,oo,po,qo,ro,so/ *-804.2316 , -.1652445 , -.5376252d-02, -4037433.d0 , *-.2091203d-06, -.4638105d-08, 0.3753368d-04, -.3853404d-02, *-.5442896d-08, 0.6484263d-13, -121.6754 , 127.5998 , *-.1486220 , -164866.6 , -.1863209d-05, 0.9622612 , * 2.097447 , -.9838123d-03, 0.3077560d-02, 0.7829503d-03/ t2 = t*t t3 = t2 * t p2 = p*p if (ibuf.eq.1) then fo2 = 13.50290120 + (-46704.69695 +.2190281453*p)/t * -6145687.892/t2 + 754294046.5/t3 else if (ibuf.eq.2) then lp = dlog(p) lt = dlog(t) fo2 = to + t*(ao + do*p + t*(fo + ho*t) + (po+qo*t)/p + * ro*lp) * + p*(bo + p*(eo + io*p) + so*lt) * + p/t*(jo/t + no*p + oo) * + ko*lt + lo*lp + co/t2 + go*dsqrt(p*t) + mo/p2 else if (ibuf.eq.3) then fo2 = dexp(dlnfo2) goto 99 else if (ibuf.eq.4) then fo2 = 16.8582 + (0.2131248*p - 53946.36)/t - 767509.6/t2 * + 0.9371923/t3 else if (ibuf.eq.5) then fo2 = ab + (bb + cb*p)/t + db/t2 + eb/t3 else call error (28,rj,ibuf,'COHHYB') end if fo2 = dexp (fo2 + dlnfo2) 99 end subroutine cohgra (fo2) c---------------------------------------------------------------------- c subroutine to compute H2O and CO2 fugacities in a COH fluid c consistent with graphite saturatuion and an oxygen fugacity c buffer specified by ibuf in routine fo2buf. c---------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10) integer ins(nsp) common / cst11 /fh2o,fco2/ cst5 /p,t,xc,u1,u2,tr,pr,rj,ps * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx * / cstcoh /xh2o,xco2,xco,xch4,xh2,x(5), * gh2o,gco2,gco,gch4,gh2,g(5),v(nsp) save tol,ins data tol,ins/.001, 1,2,3,4,5,5*0/ t2 = t * t t3 = t2 * t nit = 0 call fo2buf (fo2) c evaluate lnk CO2, CO: kco2 = dexp(.04078341613 + (47681.676177+.06372383931*p)/t * + elag - 134662.1904/t2 + 17015794.31/t3)*fo2/p kco = dexp(10.32730663 + (14062.7396777+.06372383931*p)/t * + elag - 371237.1571/t2 + 53515365.95/t3)*dsqrt(fo2)/p c get pure species fugacities call mrkpur (ins, 5) c check for graphite saturation: xco2 = kco2/gco2 xco = kco/gco if (xco2+xco.ge.1.d0) then write (*,1000) fo2,p,t fco2 = dlog(p*gco2) xco2 = 1.d0 xco = 0. goto 99 end if c evaluate other lnk's: kh2o = dexp(-7.028214449d0 + 30607.34044d0/t - 475034.4632d0/t2 * + 50879842.55d0/t3)*dsqrt(fo2) kch4 = dexp(-13.86241656d0 + (12309.03706d0+.06372383931d0*p)/t * + elag - 879314.7005d0/t2 + .7754138439d8/t3)*p c solve for x's oh2o = 2. 10 xco2 = kco2/gco2 xco = kco/gco c make quadratic 0 = qa * xh2**2 c + qb * xh2 + qc qb = kh2o * gh2/gh2o + 1.d0 qa = kch4 * gh2**2/gch4 xh2 = (-qb + dsqrt(qb**2 - 4.d0*qa*(xco + xco2 -1.d0)))/2.d0/qa xch4 = kch4 * gh2**2 * xh2**2/gch4 xh2o = kh2o * gh2 * xh2/gh2o nit = nit + 1 if (nit.gt.100) then write (*,*) ' not converging in cohgra ',nit if (xco2+xco.gt.0.99) then xco2 = 1.d0 xh2o = 1.d-20 call mrkpur (ins, 5) goto 90 else write (*,*) ' icky,icky' call error (999,xh2o,nit,'COHGRA') end if end if if (dabs(xh2o-oh2o).lt.tol*xh2o) goto 90 oh2o = xh2o call mrkmix (ins, 5) goto 10 90 xc = xco2 goto (91),hu fh2o = dlog(gh2o*p*xh2o) fco2 = dlog(gco2*p*xco2) fo2 = dlog(fo2) goto 99 91 fh2o = dlog(gh2*p*xh2) fco2 = dlog(fo2) 1000 format (' **warning ver222** routine cohgra, specified fO2 (', * g12.6,')',/,' is inconsistent with graphite saturation', * ' at P(bar)=',g12.6,' T(K)=',g12.6,/,' XCO2=1 assumed.') 99 end subroutine cohhyb (fo2) c---------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10) integer ins(nsp) common / cst11 /fh2o,fco2/ cst5 /p,t,xc,u1,u2,tr,pr,rj,ps * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx * / cstcoh /xh2o,xco2,xco,xch4,xh2,x(5), * gh2o,gco2,gco,gch4,gh2,g(5),v(nsp) * / cstchx /gmh2o,gmco2,gmch4,vm(3)/ cst26 /vol save tol,ins data tol,ins/.001, 1,2,3,4,5,5*0/ nit = 0 call fo2buf (fo2) call setup (ghh2o,ghco2,ghch4,kh2o,kco2,kco,kch4) c evaluate lnk CO2, CO: kco2 = dexp(kco2)*fo2/p kco = dexp(kco)*dsqrt(fo2)/p c check for graphite saturation: xco2 = kco2/gco2 xco = kco/gco if (xco2+xco.ge.1.d0) then write (*,1000) fo2,p,t fco2 = dlog(p*gco2) xco2 = 1.d0 xco = 0.d0 goto 99 end if c evaluate other lnk's: kh2o = dexp(kh2o)*dsqrt(fo2) kch4 = dexp(kch4)*p c solve for x's oh2o = 2.d0 10 xco2 = kco2/gco2 xco = kco/gco c make quadratic 0 = qa * xh2**2 c + qb * xh2 + qc qb = kh2o * gh2/gh2o + 1.d0 qa = kch4 * gh2**2/gch4 xh2 = (-qb + dsqrt(qb**2 - 4.d0*qa*(xco + xco2 -1.d0)))/2.d0/qa xch4 = kch4 * gh2**2 * xh2**2/gch4 xh2o = kh2o * gh2 * xh2/gh2o nit = nit + 1 if (nit.gt.100) then write (*,*) ' not converging in cohhyb ',nit if (xco2+xco.gt.0.99) then xco2 = 1.d0 xh2o = 1.d-20 call mrkpur (ins, 5) goto 90 else write (*,*) ' icky,icky' call error (999,xh2o,nit,'COHHYB') end if end if if (dabs(xh2o-oh2o).lt.tol*xh2o) goto 90 oh2o = xh2o call mrkmix (ins, 5) gh2o = ghh2o * gh2o gco2 = ghco2 * gco2 gch4 = ghch4 * gch4 goto 10 90 vol = vol + xh2o * vm(1) * + xco2 * vm(2) * + xch4 * vm(3) xc = xco2 goto (91),hu fh2o = dlog(gh2o*p*xh2o) fco2 = dlog(gco2*p*xco2) fo2 = dlog(fo2) goto 99 91 fh2o = dlog(gh2*p*xh2) fco2 = dlog(fo2) 1000 format (' **warning ver222** routine cohhyb, specified fO2 (', * g12.6,')',/,' is inconsistent with graphite saturation', * ' at P(bar)=',g12.6,' T(K)=',g12.6,/,' XCO2=1 assumed.') 99 end subroutine hsmrkp (ins, isp, jns, jsp) c--------------------------------------------------------------------- c subprogram to get volume and fugacity of pure H2O, CO2, and CH4 c fluids from HSMRK EOS of Kerrick and Jacobs (1981) and Jacobs and c Kerrick (1981). c--------------------------------------------------------------------- implicit double precision (a-h,p,r-y) parameter (nsp=10) integer ins(nsp), jns(3), isp, jsp, ij(3) double precision fg(3) common / cst85 /p,t,t12,r/ cst5 /pbar,tk,xc,u1,u2,tr,pr,rcal,ps * / cst11 /f(2)/ cstcoh /x(nsp),g(nsp),v(nsp) * / cstchx /gm(3),vm(3) save bw, bc, bm data bw, bc, bm, ij /29., 58., 60., 1, 2, 4/ t = tk p = pbar t15 = dsqrt(t**3) t12 = dsqrt(t) rtt = r*t15 t2 = t*t call mrkpur (ins,isp) do 10 k = 1, jsp i = jns(k) j = ij(i) if (i.eq.1) then b = bw c = 290.78d6-0.30276d6*t+0.00014774d6*t2 d = -8374.d6+19.437d6*t-0.008148d6*t2 e = 76600.d6-133.9d6*t+0.1071d6*t2 else if (i.eq.2) then b = bc c = 28.31d6+0.10721d6*t-0.00000881d6*t2 d = 9380.d6-8.53d6*t+0.001189d6*t2 e = -368654.d6+715.9d6*t+0.1534d6*t2 else b = bm c = 13.403d6 + 9.28d4 * t + 2.7 * t2 d = 5.216d9 - 6.8d6 * t + 3.28d3 * t2 e = -2.3322d11 + 6.738d8 * t + 3.179d5 * t2 end if vm(i) = -v(j) gm(i) = g(j) call nurap (b,c,d,e,yz,v(j)) vm(i) = vm(i) + v(j) fg(i) = dlog(p) + fugp (rtt,b,yz,c,d,e,v(j)) g(j) = dexp(fg(i))/p 10 if (i.lt.3) f(i) = fg(i) end subroutine nurap (b,c,d,e,yz,vi) c---------------------------------------------------------- c newton-raphson iteration to solve for hsmrk volume, my c idea here was to evaluate all the constants outside of c the iteration loop, which might make sense if a lot of c iterations were necessary, as it is i'm unsure. someone c should do a comparison with my old nurap routine (xnurap). c jadc, july 96. c----------------------------------------------------------- implicit double precision (a-h,p-y) common /cst85 /p,t,t12,r s1 = r*t*t12 s2 = b*s1 s3 = t12*p*b p0 = -256.d0*s1 b2 = b*b q0 = 256.d0*t12*p q1 = 256.d0*(s3 - s1) q2 = (-160.d0*s3 - 512*s1)*b + 256.d0*c q3 = (-80.d0*s3 + p0)*b2 + 256.d0*d q4 = ((65*s3 + 8*s1)*b - 160.d0*c)*b2 + 256.d0*e q5 = -b2*(((14*s3 - 15*s1)*b - 80.d0*c)*b + 160.d0*d) q6 = b2*((((s3 + 6.d0*s1)*b - 15*c)*b + 80.d0*d)*b - 160.d0*e) q7 = b**3*(((-s1*b + c)*b - 15*d)*b + 80.d0*e) q8 = b**4*(-15.d0*e + d*b) q9 = e*b**5 p1 = 512.d0*c-768.d0*s2 p2 = (-832.d0*s2-256.d0*c)*b+768.d0*d p3 = ((-368.d0*s2 - 64.d0*c)*b - 256.d0*d)*b + 1024.d0*e p4 = -b*(((33.d0*s2 - 64.d0*c)*b + 224.d0*d)*b + 256.d0*e) p5 = 2.d0*b2*(b*((s2-c)*7.d0*b + 72.d0*d) - 192.d0*e) p6 = -b**3*(b*((s2-c)*b + 29.d0*d) - 224.d0*e) p7 = 2.d0*b**4*(d*b-22.d0*e) p8 = 3.d0*q9 do 10 k = 1, 50 cor = ((((((((((q0*vi+q1)*vi+q2)*vi+q3)*vi+q4)*vi+q5)*vi * +q6)*vi+q7)*vi+q8)*vi+q9)*vi)/((((((((p0*vi+p1)*vi * +p2)*vi+p3)*vi+p4)*vi+p5)*vi+p6)*vi+p7)*vi+p8) vi = vi + cor 10 if (dabs(cor).lt.0.01) goto 99 99 yz = vi * p/r/t end subroutine xnurap (b,c,d,e,yz,vi) implicit double precision (a-h,p,r-y) common /cst85 /p,t,t12,r do 10 k = 1, 50 y = b/4.d0/vi x = 1.d0 - y bi = vi + b bi2 = bi*bi vi2 = vi*vi vi3 = vi2*vi x3 = x**3 y3 = y**3 pn = 1.d0 + y + y*y - y3 pa1 = c + d/vi + e/vi2 d1 = -.75d0* b/vi3/x/x3 - 1.d0/vi2/x3 d3 = 1.d0/vi/x3 * (-b/4.d0/vi2 - 2.d0*b*b/16.d0/ vi3 * + 0.046875d0 * b**3/vi/vi3 ) df = ( pn * d1 + d3 )*r*t * - (pa1 * (-1.d0/vi/bi2 - 1.d0/vi2/bi) * + (1.d0 /vi/ bi) * (-d/vi2 - 2.d0* e/vi3))/t12 v = vi - (pn/vi/x3 * r*t - pa1/t12/vi/bi - p)/df diff = dabs(v-vi) vi = v 10 if (diff.lt.0.01) goto 99 99 yz = vi * p/83.1441d0/t end subroutine mrkmix (ins, isp) c----------------------------------------------------------------------- c subroutine to calculate the log fugacities and volume of mixed c species fluids using the RK/MRK EoS. c input: c ins(i) - pointers indicating the species are to be calculated. c isp - the number of species to be calculated. c p,t - input from cst5, p(bars), t(K) c output (to common cstcoh): c g(i) - fugacity coefficient of ith species c v(i) - volume of the ith species c species indices: c 1 = H2O c 2 = CO2 c 3 = CO c 4 = CH4 c 5 = H2 c 6 = H2S c 7 = O2 c 8 = SO2 c 9 = COS c 10 = C2H6 c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (nsp=10) integer ins(nsp) double precision b(nsp),a(nsp),f(nsp),aj2(nsp),ev(3) common/ cst5 /p,t,xco2,u1,u2,tr,pr,rcal,ps * / cstcoh /x(nsp),g(nsp),v(nsp)/ cst26 /vol/ cst11 /fg(2) save a, b, r data r / 83.1441 /, * a /0., 0., 16.98d6, 32.154811d6, 2.821d5, 89.d6, * 17.4d6,133.1d6, 130.d6, 90.d6/, * b /14.6, 29.7, 27.38, 29.681, 15.7699, 29.94, * 22.1, 37.4, 43., 20./ t2 = t*t t3 = t2*t t4 = t3*t dsqrtt = dsqrt(t) rt = r*t a(1) = .1452535403d8 +306893.3587*t * -307.9995871*t2 +.09226256008*t3-.2930106337d-5*t4 c compute dispersion term for co2 a(2) = 92935540.d0 - 82130.73d0*t + 21.29d0*t2 ch = dexp(-11.218d0 + 6032./t - 2782000./t2 + 4.708d8/t3) * * 6912.824964 * t2 * dsqrtt + 79267647.d0 c composition dependent mrk-terms bx = 0.d0 aij = 0.d0 do 400 k = 1, isp i = ins(k) aj2(i) = 0. 400 bx = bx + b(i)*x(i) do 200 k = 1, isp i = ins(k) do 300 l = 1, isp j = ins(l) if (i.eq.1.and.j.eq.2.or.i.eq.2.and.j.eq.1) then aij = aij + x(i)*x(j)*ch/2.d0 aj2(i) = aj2(i) + x(j)*ch else aij = aij + x(i)*x(j)*dsqrt(a(i)*a(j)) aj2(i) = aj2(i) + 2.d0*x(j)*dsqrt(a(i)*a(j)) end if 300 continue 200 continue c solve for molar volume of the mix c1 = -rt/p c3 = -aij*bx/p/dsqrtt c2 = c1*bx + aij/dsqrtt/p - bx*bx call roots3 (c1,c2,c3,ev,vmin,vmax,iroots) if (iroots.eq.3) then vol = vmax else vol = ev(1) end if c segment to compute the c fugacities. d2 = dlog((vol+bx)/vol) d1 = rt*dsqrtt*bx d3 = d2 -bx/(bx+vol) d4 = vol -bx d5 = aij*d3/d1/bx d6 = dlog(rt/d4) do 70 i = 1, isp l = ins(i) if (x(l).gt.0.d0) goto 80 g(l) = 1.d0 goto 75 80 f(l) = dlog(x(l))+b(l)/d4-aj2(l)/d1*d2+b(l)*d5+d6 g(l) = dexp(f(l))/p/x(l) 75 if (l.lt.3) fg(l) = f(l) 70 continue end subroutine crkh2o (pbar,t,vcrk,fh2o) c----------------------------------------------------------------------- c compute ln(f[h2o], bar) and volume h2o (kj/bar) from CORK EoS Holland c & Powell CMP 109:265-273. Input pbar - pressure (bars); tk - temp (K). c J.A.D. Connolly, 1992. c----------------------------------------------------------------------- implicit double precision (a-h,o-z) dimension x(3) data b,r,p0 /1.465d0,8.314d-3,2.d0/ p = pbar/1.d3 rt = r*t rtp = rt/p t12 = dsqrt(t) if (t.lt.695.) then psat = -13.627d-3 + 7.29395d-7*t**2 - 2.34622d-9*t**3 * + 4.83607d-15*t**5 if (p.lt.psat.and.t.lt.673.d0) then a = 1113.4d0 + 5.8487d0*(673.-t) - 2.1370d-2*(673.-t)**2 * + 6.8133d-5*(673.-t)**3 else if (t.lt.673.) then a = 1113.4d0 - 0.88517*(673.-t) + 4.53d-3*(673.-t)**2 * - 1.3183d-5*(673.-t)**3 else a = 1113.4d0 - .22291d0*(t - 673.) * - 3.8022d-4*(t-673.)**2 + 1.7791d-7*(t-673.)**3 end if end if else if (t.lt.673.) then a = 1113.4d0 - 0.88517*(673.-t) + 4.53d-3*(673.-t)**2 * - 1.3183d-5*(673.-t)**3 else a = 1113.4d0 - .22291d0*(t - 673.) * - 3.8022d-4*(t-673.)**2 + 1.7791d-7*(t-673.)**3 end if end if a1 = -rtp a2 = a/t12/p - b*(rtp+b) a3 = -a*b/t12/p call roots3 (a1,a2,a3,x,xmin,xmax,iroots) if (iroots.eq.1) then vol = x(1) else if (p.lt.psat) then vol = xmax else if (t.lt.700.) then vol = xmin else c cork has 3 roots sometimes at c high P, this is a fitting artifact c and in my experience the only real c root is positive. do 1 i = 1, 3 if (x(i).gt.0.d0) then vol = x(i) goto 2 end if 1 continue end if end if 2 cc = a/b/rt/t12 gam = vol/rtp -1.d0 - dlog((vol-b)/rtp) - cc*dlog(1.d0+b/vol) if (p.gt.p0) then c cork classic c c = -3.02565d-2 - 5.343144d-6*t c d = -3.2297554d-3 + 2.215221d-6*t c vcrk = vol + c*dsqrt(dp) + d*dp c gam = gam + (dp**1.5d0*2.d0*c/3.d0 + dp**2*d/2.d0)/rt c new cork c = 2.1111d-3 - 5.20907d-7*t d = -9.17623d-2 + 1.12719d-5*t dp = p - p0 vcrk = vol + c*dp + d*dsqrt(dp) + dp**0.25d0*7.40395d-2 gam = gam + (dp**2*c/2.d0 * + dp**1.5d0*d*0.66666666666666667 * + dp**1.25d0*0.0592316d0)/rt end if if (t.lt.695.d0.and.p.gt.psat) then p = psat a1 = -rtp a2 = -(b*(rt+b*p)-a/t12)/p a3 = -(a*b/t12)/p call roots3 (a1,a2,a3,x,vmin,vmax,iroots) if (iroots.eq.1) goto 99 gam1 = vmax/rtp-1.d0-dlog((vmax-b)/rtp)-cc*dlog(1.d0+b/vmax) gam2 = vmin/rtp-1.d0-dlog((vmin-b)/rtp)-cc*dlog(1.d0+b/vmin) gam = gam1 - gam2 + gam end if 99 fh2o = gam + dlog(pbar) end subroutine crkco2 (pbar,t,fco2) c----------------------------------------------------------------------- c compute ln(f[co2], bar) from Eq 8 of Holland c & Powell CMP 109:265-273. Input pbar - pressure (bars); tk - temp (K). c J.A.D. Connolly, 1992. c----------------------------------------------------------------------- implicit double precision (a-h,o-z) data b,r,tc,pc /3.7852d0,8.314d-3,304.2d0,0.0738d0/ p = pbar/1.d3 rt = r*t bp = b*p a = (5.45963d-5*tc - 8.6392d-6*t)*tc**(1.5d0)/pc c = (-3.30558d-5*tc + 2.30524d-6*t)/pc**(1.5d0) d = (6.93054d-7*tc - 8.38293d-8*t)/pc/pc fco2 = dlog(pbar) + (bp + a/b/dsqrt(t)*(dlog(rt+bp) * - dlog(rt+2.d0*bp)) + 2.d0*c*p**(1.5d0)/3.d0 * + d/2.d0*p*p)/rt end subroutine roots3 (a1,a2,a3,x,vmin,vmax,iroots) c--------------------------------------------------------- c returns real roots (in x) of: x**3 + a1*x**2 + a2*x + a3 c--------------------------------------------------------- implicit double precision (a-h,o-z) dimension x(3) qq = (a1**2-3.d0*a2)/9.d0 rr = (2.d0*a1**3-9.d0*a1*a2+27.d0*a3)/54.d0 a5 = a1/3.d0 dif = qq**3-rr**2 if (dif.ge.0.d0) then phi = dacos(rr/dsqrt(qq**3)) a4 = -2.d0*dsqrt(qq) a6 = phi/3.d0 dphi = 0.d0 vmin = 1.d9 vmax = -1.d9 do 10 i = 1, 3 v = a4*dcos(a6+dphi) - a5 if (v.gt.vmax) vmax = v if (v.lt.vmin) vmin = v x(i) = v 10 dphi = dphi + 2.094395102497915d0 iroots = 3 else a7 = (dsqrt(-dif) + dabs(rr))**(1.d0/3.d0) x(1) = -rr/dabs(rr)*(a7 + qq/a7) - a5 iroots = 1 end if end subroutine lomrk (ins, isp) c----------------------------------------------------------------------- c a semi-fudged (see comments below) MRK to calculate the fugacities c of mixed species fluids in the vicinity of the water critical point. c input: c ins(i) - pointers indicating the species are to be calculated. c isp - the number of species to be calculated. c p,t - input from cst5, p(bars), t(K) c output (to common cstcoh): c g(i) - fugacity coefficient of ith species c v(i) - volume of the ith species c species indices: c 1 = H2O c 2 = CO2 c 3 = CO c 4 = CH4 c 5 = H2 c 6 = H2S c 7 = O2 c 8 = SO2 c 9 = COS c 10 = C2H6 c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (nsp=10) integer ins(nsp) double precision b(nsp),a(nsp),f(nsp),aj2(nsp),ev(3) common/ cst5 /p,t,xco2,u1,u2,tr,pr,rcal,ps * / cstcoh /x(nsp),g(nsp),v(nsp) * / cst26 /vol/ cst11 /fg(2) save a, b, r data r / 83.1441d0 /, * a /0., 0., 16.98d6, 32.154811d6, 2.821d5, 89.d6, * 17.4d6,133.1d6, 130.d6, 90.d6/, * b /14.6, 29.7, 27.38, 29.681, 15.7699, 29.94, * 22.1, 37.4, 43., 20./ t2 = t*t dsqrtt = dsqrt(t) rt = r*t c this is a cheap trick the a-fit c was derived using b=15, then b c was adjusted to lower the H2-H2O c critical T, however the fugacities c seem reasonable so who knows? a(1) = .3930568949d9 - .1273025840d7 * t + 2049.978752 * t2 * -1.122350458 * t2 * t c compute dispersion term for co2 a(2) = 0.9293554d8 - 0.8213073d5*t + 0.2129d2*t2 ch = dexp(-11.218d0 + 6032.d0/t - 2782.d3/t2 + 4.708d8/t2/t) * * 6912.824964d0 * t2 * dsqrtt + 79267647.d0 c composition dependent mrk-terms bx = 0.d0 aij = 0.d0 do 400 k = 1, isp i = ins(k) aj2(i) = 0.d0 400 bx = bx + b(i)*x(i) do 200 k = 1, isp i = ins(k) do 300 l = 1, isp j = ins(l) if (i.eq.1.and.j.eq.2.or.i.eq.2.and.j.eq.1) then aij = aij + x(i)*x(j)*ch/2.d0 aj2(i) = aj2(i) + x(j)*ch else aij = aij + x(i)*x(j)*sqrt(a(i)*a(j)) aj2(i) = aj2(i) + 2.d0*x(j)*sqrt(a(i)*a(j)) end if 300 continue 200 continue c solve for molar volume of mixture c1 = -rt/p c3 = -aij*bx/p/dsqrtt c2 = c1*bx + aij/dsqrtt/p - bx*bx call roots3 (c1,c2,c3,ev,vmin,vmax,iroots) if (iroots.eq.3) then vol = vmax else vol = ev(1) end if c segment to compute the c fugacities. d2 = dlog((vol+bx)/vol) d1 = rt*dsqrtt*bx d3 = d2 -bx/(bx+vol) d4 = vol -bx d5 = aij*d3/d1/bx d6 = dlog(rt/d4) do 70 i = 1, isp l = ins(i) if (x(l).gt.0.) goto 80 f(l) = 0. g(l) = 1. goto 75 80 f(l) = dlog(x(l))+b(l)/d4-aj2(l)/d1*d2+b(l)*d5+d6 g(l) = dexp(f(l))/p/x(l) 75 if (l.lt.3) fg(l) = f(l) 70 continue end subroutine mrkpur (ins, isp) c----------------------------------------------------------------------- c subroutine to calculate the log fugacities and volume of single c species fluids using the RK/MRK EoS. c input: c ins(i) - pointers indicating the species are to be calculated. c isp - the number of species to be calculated. c p,t - input from cst5, p(bars), t(K) c output (to common cstcoh): c g(i) - fugacity coefficient of ith species c v(i) - volume of the ith species c species indices: c 1 = H2O c 2 = CO2 c 3 = CO c 4 = CH4 c 5 = H2 c 6 = H2S c 7 = O2 c 8 = SO2 c 9 = COS c 10 = C2H6 c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (nsp=10) double precision b(nsp),a(nsp),f(nsp),ev(3) integer ins(nsp), isp common/ cst5 /p,t,xco2,u1,u2,tr,pr,rc,ps/ cst26 /vol * / cstcoh /x(nsp),g(nsp),v(nsp)/ cst11 /fg(2) save a, b, r data r / 83.1441d0 /, * a /0., 0., 16.98d6, 32.154811d6, 2.821d5, 89.d6, * 17.4d6,133.1d6, 130.d6, 0.d6/, * b /14.6, 29.7, 27.38, 29.681, 15.7699, 29.94, * 22.1, 37.4, 43., 0./ t2 = t*t t3 = t2*t t4 = t3*t dsqrtt = dsqrt(t) rt = r*t do 70 k = 1, isp i = ins(k) if (i.eq.1) then a(1) = .1452535403d8 +306893.3587*t -307.9995871d0*t2 * +.09226256008d0*t3-.2930106337d-5*t4 else if (i.eq.2) then c compute dispersion term for co2 a(2) = 92935540.d0 - 82130.73d0*t + 21.29d0*t2 end if aij = a(i) bx = b(i) c1 = -rt/p c3 = -aij*bx/p/dsqrtt c2 = c1*bx + aij/dsqrtt/p - bx*bx call roots3 (c1,c2,c3,ev,vmin,vmax,iroots) if (iroots.eq.3) then vol = vmax else vol = ev(1) end if c compute fugacities. d2 = dlog((vol + bx)/vol) d1 = rt*dsqrtt*bx d3 = d2 - bx/(bx + vol) d4 = vol - bx v(i) = vol f(i) = bx/d4 - 2.d0*aij/d1*d2 + bx*aij*d3/d1/bx + dlog(rt/d4) if (i.lt.3) fg(i) = f(i) 70 g(i) = dexp(f(i))/p end subroutine setfs2 (fs2,kh2s,kso2,kcos) c----------------------------------------------------------------------- c fs2, kh2s, kso2, kcos calculation. c----------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) common / cst5 /p,t,xo,u1,u2,tr,pr,rj,ps * / cst100 /rat,elag,gz,gy,gx,ibuf,hu,hv,hw,hx c get sulfur fugacity according to c the value of ibuf: if (ibuf.eq.1) then c get po-py 1/2 ln sulfur fugacity: c from simon poulson: fs2 = .005388049*t + 10.24535 - 15035.91/t + 0.03453878/t*p else if (ibuf.eq.2) then c get suflur fugacity from fe/s ratio c of po, expression derived from c toulmin & barton 1964, with p c correction from Craig & Scott 1982. c (this only takes into account V(S2), c hope it's right. xf = rat/( rat + 1.d0) fs2 = 197.6309*xf + 45.2458*dsqrt(1.d0- 1.9962*xf) - 94.33691 * + (80624.79 + 0.2273782*p - 197630.9*xf)/t else fs2 = rat/2.d0 end if c k's from ohmoto and kerrick kh2s = 10115.3d0/t - 0.791d0 * dlog (t) + 0.30164d0 kcos = 10893.52964d0/t - 9.988613730d0 kso2 = 43585.63147d0/t - 8.710679055d0 end subroutine hosmrk (fo2, fs2) c----------------------------------------------------------------------- c program to calculate H-O-S speciation as a function of XO using c an MRK/HSMRK hybrid. Specifically for high (po+py) fs2. c Species are H2 H2O H2S SO2. c----------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10) integer ins(nsp), jns(3) common / cst11 /fh2o,fco2/ cst5 /p,t,xo,u1,u2,tr,pr,rj,ps * / cstcoh /xh2o,xco2,xco,xch4,xh2,xh2s,xo2,xso2,x(2), * gh2o,gco2,gco,gch4,gh2,gh2s,go2,gso2,g(2),v(nsp) * / cstchx /gmh2o,vm(5)/ cst26 /vol * / cst100 /rat,elag,gz,gy,gx,ibuf,hu,hv,hw,hx save ins,jns data ins, jns/ 1,5,6,8,6*0,1,2*0/ c check if xo is <1, >0, c reset if necessary if (xo.gt.0.9999999999) then xo = 0.9999999999 else if (xo.lt.0.0000000001) then xo = 0.0000000001 end if c this fs2 = 1/2 ln (fs2), c k's are ln(k) call setfs2 (fs2,kh2s,kso2,kcos) kh2o = -7.028214449 + 30607.34044/t - 475034.4632/t/t * + 50879842.55/t/t/t c1 = dexp (kh2s + fs2) c3 = dexp (kso2 - 2.d0*kh2o + fs2)/p xom = xo - 1.d0 xop = xo + 1.d0 xos = xo * xo a = 8.*xo*xom**2 b = 4.d0*xom*(3.d0*xos + 1.d0) c0 = 2.d0*xo*(3.d0*xos - 1.d0) + 4.d0 d = xom*xop**2 c get first gammas call hsmrkp (ins, 1, jns, 1) c hybrid gamma factor for water: ghh2o = gh2o/gmh2o call mrkpur (ins, 4) gh2o = gh2o * ghh2o c get first guess: turd = 1.d0/3.d0 if (xo.lt.turd) then xl = 2.d0 * xo/(1.d0-xo) else if (xo.gt.turd) then xl = 2.d0 * (1.d0-xo)/xop end if c outer iteration loop: do 30 j = 1, 20 c2 = gh2/gh2s c4 = gh2o**2/gso2/gh2**2 c = c0 - 8.*c3*c4*(1.d0+c1*c2)**2 xh2o = xl xi = xl do 10 i = 1, 30 c inner iteration loop: f = a + b * xh2o + c * xh2o**2 + d * xh2o**3 df = b + 2.d0 * c * xh2o + 3.d0 * d * xh2o ** 2 xh2o = xi - f/df xh2 = (-xo*xh2o-xh2o-2.d0*xop)/(1.d0+c1*c2)/2.d0 xh2s = c1*c2*xh2 xso2 = c3*c4*xh2o**2/xh2**2 if (dabs((xi-xh2o)/xh2o).lt.1.d-6) goto 20 c write (*,*) xh2o,xi if (xh2o.ge.1.d0) xh2o = xi + (1.d0-xi)/2.d0 10 xi = xh2o goto 50 20 xh2 = -0.5d0*(xo*xh2o+xh2o+2.d0*xo-2.d0)/(1.d0+c1*c2) if (xh2.le.1.d-5.and.xo.gt.turd) goto 50 xh2s = c1*c2*xh2 xso2 = c3*c4*xh2o**2/xh2**2 if (j.gt.1.and.dabs((xl-xh2o)/xh2o).lt.1.d-6) goto 40 call mrkmix (ins, 4) gh2o = gh2o * ghh2o xl = xh2o xi = xh2o 30 continue write (*,*) ' no convergence in outer loop of hosmrk' stop 40 fh2o = dlog(gh2o*p*xh2o) fo2 = 2.d0 * (fh2o - dlog(gh2 * p * xh2) - kh2o) vol = vol + xh2o * vm(3) goto 99 c if xo.gt.turd, at high fs2 c the fluid will be binary: 50 xh2o = -2.d0*xom/xop xso2 = 1.d0 - xh2o xh2s = 0. xh2 = 0. call mrkmix (ins, 4) vol = vol + xh2o * vm(3) gh2o = gh2o * ghh2o fh2o = dlog(gh2o*p*xh2o) fo2 = dlog(gso2*p*xso2) - kso2 - fs2 99 end subroutine hosrk5 (fo2) c----------------------------------------------------------------------- c program to calculate H-O-S speciation as a function of XO using c an MRK/HSMRK hybrid. Species are H2 O2 H2O H2S SO2. c----------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10) integer ins(nsp), jns(3) common / cst11 /fh2o,fco2/ cst5 /p,t,xo,u1,u2,tr,pr,rj,ps * / cstcoh /xh2o,xco2,xco,xch4,xh2,xh2s,xo2,xso2,x(2), * gh2o,gco2,gco,gch4,gh2,gh2s,go2,gso2,g(2),v(nsp) * / cstchx /gmh2o,vm(5)/ cst26 /vol * / cst100 /rat,elag,gz,gy,gx,ibuf,hu,hv,hw,hx save ins,jns data ins, jns/ 1,5,6,7,8,5*0,1,2*0/ c check if xo is <1, >0, c reset if necessary if (xo.gt.0.9999999999) then xo = 0.9999999999 else if (xo.lt.0.0000000001) then xo = 0.0000000001 end if c this fs2 = 1/2 ln (fs2), c k's are ln(k) call setfs2 (fs2,k1,k2,kcos) c kh2o from robie: k3 = dexp(-7.028214449d0+30607.34044d0/t-475034.4632d0/t/t * +50879842.55d0/t/t/t) xom = xo - 1.d0 xop = xo + 1.d0 xos = xo * xo c1 = dexp(k1 + fs2) c3 = dexp(k2 + fs2) c5 = 1.d0/p/k3/k3 a = -8.*xo*xom**3 b = -4.d0*(3.d0*xos+1.d0)*xom**2 c0 = 2.d0*xom*(-xop * (3.d0*xo*xom + 2.d0)) c7 = 8.*xom*c5 d = -xom**2 * xop**2 c get first gammas call hsmrkp (ins, 1, jns, 1) c hybrid gamma factor for water: ghh2o = gh2o/gmh2o call mrkpur (ins, 5) gh2o = gh2o * ghh2o c get first guess: turd = 1.d0/3.d0 if (xo.lt.turd) then if (xo.gt.0.3333333333333333) xo = 0.3333333333333333 xl = 2.d0 * xo/(1.d0-xo) else if (xo.gt.turd) then if (xo.lt.0.3333334) xo = 0.3333334 xl = 2.d0 * (1.d0-xo)/(1.d0 + xo) end if c outer iteration loop: do 30 j = 1, 20 c2 = gh2/gh2s c4 = go2/gso2 c6 = gh2o**2/gh2**2/go2 c = c0 + c7*c6*(1.d0+c1*c2)**2*(1.d0+c3*c4) xh2o = xl xi = xl do 10 i = 1, 30 c inner iteration loop: f = a + b * xh2o + c * xh2o**2 + d * xh2o**3 df = b + 2.d0 * c * xh2o + 3.d0 * d * xh2o ** 2 xh2o = xi - f/df xh2 = -0.5d0*(xo*xh2o+xh2o+2.d0*xo-2.d0)/(1.d0+c1*c2) xh2s = c1*c2*xh2 xo2 = c5*c6*(xh2o*xh2o)/(xh2*xh2) xso2 = c3*c4*xo2 if (dabs((xi-xh2o)/xh2o).lt.1.d-6) goto 20 c write (*,*) xh2o,xi if (xh2o.ge.1.d0) xh2o = xi + (1.d0-xi)/2.d0 10 xi = xh2o write (*,*) ' no convergence in inner loop of homrk' write (*,*) p,t,xo,xh2o stop 20 xh2 = -0.5d0*(xo*xh2o+xh2o+2.d0*xo-2.d0)/(1.d0+c1*c2) xh2s = c1*c2*xh2 xo2 = c5*c6*(xh2o*xh2o)/(xh2*xh2) xso2 = c3*c4*xo2 if (j.gt.1.and.dabs((xl-xh2o)/xh2o).lt.1.d-6) goto 40 call mrkmix (ins, 5) gh2o = gh2o * ghh2o xl = xh2o xi = xh2o 30 continue write (*,*) ' no convergence in outer loop of hosmrk' stop 40 fh2o = dlog(gh2o*p*xh2o) vol = vol + xh2o * vm(3) if (xo2.lt.xh2) then fo2 = 2.d0 * (fh2o - dlog(gh2 * p * xh2) - dlog(k3)) else fo2 = dlog (go2 * p * xo2) end if end subroutine homrk (fo2) c----------------------------------------------------------------------- c program to calculate H-O speciation as a function of XO using c an MRK/HSMRK hybrid. c----------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10) integer ins(nsp), jns(3) common / cst11 /fh2o,fco2/ cst5 /p,t,xo,u1,u2,tr,pr,rj,ps * / cstcoh /xh2o,xco2,xco,xch4,xh2,xh2s,xo2,x(3), * gh2o,gco2,gco,gch4,gh2,gh2s,go2,g(3),v(nsp) * / cstchx /gmh2o,vm(5)/ cst26 /vol save ins,jns data ins, jns/ 1,5,7,7*0,1,2*0/ c check if xo is <1, >0, c reset if necessary if (xo.gt.0.9999999999) then xo = 0.9999999999 else if (xo.lt.0.0000000001) then xo = 0.0000000001 end if k = dexp(-7.028214449d0 + 30607.34044d0/t - 475034.4632d0/t/t * + 50879842.55d0/t/t/t) turd = 1.d0/3.d0 c1 = 1.d0/dsqrt(p)/k c10 = (xo - 1.d0)/2.d0 c11 = 1.d0 - xo c12 = 1.d0 + c10 c get first gammas call hsmrkp (ins, 1, jns, 1) c hybrid gamma factor for water: ghh2o = gh2o/gmh2o call mrkpur (ins, 3) gh2o = gh2o * ghh2o if (xo.lt.turd) then if (xo.gt.0.3333333333333333) xo = 0.3333333333333333 xl = 2.d0 * xo/c11 else if (xo.gt.turd) then if (xo.lt.0.3333334) xo = 0.3333334 xl = 2.d0 * c11/(1.d0 + xo) end if c outer iteration loop: do 30 j = 1, 10 c c13 = c1 * gh2o/gh2/dsqrt(go2) c14 = c10 * c13/2.d0 xi = xl do 10 i = 1, 10 c inner iteration loop: xo2 = xo + c10 * xh2o if (xo2.gt.1.d-9) then x12 = dsqrt (xo2) xh2o = xi + * (c11 - c12 * xh2o - c13 * xh2o/x12) / * (c12 + c13 * x12 + c14 * xh2o/x12) else xh2o = 2.d0 * xo/c11 end if if (dabs((xi-xh2o)/xh2o).lt.1.d-6) goto 20 if (xh2o.ge.1.d0) xh2o = xi + (1.d0-xi)/2.d0 10 xi = xh2o write (*,*) ' no convergence in inner loop of homrk' write (*,*) p,t,xo,xh2o stop 20 if (xo2.lt.0.) xo2 = 0. xh2 = 1.d0 - xo2 - xh2o if (j.gt.1.and.dabs((xl-xh2o)/xh2o).lt.1.d-6) goto 40 call mrkmix (ins, 3) gh2o = gh2o * ghh2o xl = xh2o xi = xh2o 30 continue write (*,*) ' no convergence in outer loop of homrk' stop 40 fh2o = dlog(gh2o*p*xh2o) vol = vol + xh2o * vm(3) if (xo2.lt.xh2) then fo2 = 2.d0 * (fh2o - dlog(gh2 * p * xh2) - dlog(k)) else fo2 = dlog (go2 * p * xo2) end if end function fugp (rtt,b,yz,c,d,e,v) c--------------------------------------------------------------------- c streamlined JADC July 96. c--------------------------------------------------------------------- implicit double precision (a-h,p,r-y) y = b/4.d0/v vpb = v + b dvbv = dlog(vpb/v) y1 = 1.d0 - y fugp = ((4.d0-3.d0*y)*y + (2.d0-y)*2.d0*y/y1)/y1/y1 * + (d*(dvbv/b + (4.d0*y + 2.d0)/vpb - 3.d0/v) * - c*(dvbv + b/vpb) * + e*((4.d0/b-2.d0/v)/v - dvbv/b/b * + ((2.d0*y - 1.5d0)/v - 3.d0/b)/vpb) * )/rtt/b - dlog(yz) end function fug (rtt,cij,dij,eij,x1,x2,b,yz,c,d,e,b1,c1,d1,e1) c--------------------------------------------------------------------- c streamlined JADC July 96. c--------------------------------------------------------------------- implicit double precision (a-h,p,r-y) common / cst26 /v y = b/4.d0/v vpb = v + b dvbv = dlog(vpb/v) y1 = 1.d0 - y dvbvb = dvbv/b vi = 1.d0/v vi2 = .5d0/v/v fug = ((4.d0-3.d0*y)*y + b1/b*(2.d0-y)*y*2.d0/y1)/y1/y1 * + (c*b1*(dvbvb - 1.d0/vpb) * - (c1*x1+cij*x2)*2.d0*dvbv * + (2.d0*(x1*d1+dij*x2)+d)*(dvbvb - vi) * + d*b1*((vi+2.d0/b)/vpb - 2.d0*dvbvb/b) * + 2.d0*(e1*x1+eij*x2+e)*((vi - dvbvb)/b - vi2) * + e*b1*((vi2 - (1.5d0/v+3.d0/b)/b)/vpb + 3.d0*dvbvb/b/b) * )/rtt/b-dlog(yz) end subroutine brmrk c----------------------------------------------------------------------- c brmrk routine to calculate molar volume and fugacity of co2 from c bottinga and richet 1981, ajs v 281. c j. connolly 1988. c c for the unitiated, input and output is done through common blocks c cst5, cst11, and cst26, the relevant variables are: c c for input: c p = pressure in bars c t = temperature in kelvins c r = gas constant c xco2 = mole fraction of co2 in fluid phase c for output: c v = molar volume (cm3/mole) at p and t c fco2 = the natural log of the co2 fugacity c----------------------------------------------------------------------- implicit double precision (a-g,o-y) double precision vdpdv common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps * / cst11 /fh2o,fco2/ cst26 /v c external function to calculate c the product v(dp/dv) external vdpdv if (xco2.ne.1.d0) then xco2 = 1 call warn (35,xco2,0,'BRMRK ') end if c call brvol to calculate volume c of co2. v1 = brvol (1.d0,t) c v2 = brvol (p,t) c call qromb to evaluate integral c of (dp/dv)dv if (v2.ge.180.) then c if b&r equation is continuous c through limits v1 to v do c integration in one chunk. call qromb (vdpdv,v1,v2,fco2) c else if (v2.gt.47.22) then c else do the integration in parts call qromb (vdpdv,v1,180.d0,f1) call qromb (vdpdv,180.d0,v2,f2) fco2 = f1 + f2 c else c call qromb (vdpdv,v1,180.d0,f1) call qromb (vdpdv,180.d0,47.22d0,f2) call qromb (vdpdv,47.22d0,v2,fco2) fco2 = fco2 + f1 + f2 end if v = v2 fco2 = fco2/(10.*r*t) end function brvol (p,t) c----------------------------------------------------------------------- c brvol function to calculate molar volume of co2 by newton-raphson c iteration of eq 3 of bottinga and richet 1981, ajs v 281. c input: p = pressure in bars, t = temperature in kelvins, v = molar c volume of co2 (cm3) estimated by the mrk equation of state. c output: brvol = molar volume of co2 in (cm3/mole) c j. connolly 1988. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (nsp=10) integer ins(nsp) common/ cst26 /v save rbar, dv, a1, a2, a3, ins data rbar,dv/83.143d0,0.001/,a1,a2,a3/6.566d7,7.276d7,37.3/, * ins/2,9*0/ c--------------------------------------------------------------------- c get initial v estimate. call mrkpur (ins, 1) rt = rbar * t t12 = dsqrt(t) itic = 0 c brvol appoximates the derivative dp/dv c by finite difference with dv (in cm3) dv = 0.00005 c choose parameters for b&r's equation 10 if (v.le.47.22) then b1 = 1.856669 b2 = 0.0637935 else if (v.lt.180.) then b1 = 11.707864 b2 = 0.363955 else b1 = 7.352629 b2 = 0.241413 end if c at v: b = (dlog(v/a3)+b1)/b2 a = a1*( (a3/v)**3 - (a3/v)**6 ) + a2 c and at vp: vp = v + dv c bp = (dlog(vp/a3)+b1)/b2 ap = a1*( (a3/vp)**3 - (a3/vp)**6 ) + a2 c fv = rt/(v-b) - a/(v*(v+b)*t12) - p dpdv = (fv - (rt/(vp-bp) - ap/(vp*(vp+bp)*t12) - p))/dv c refine estimate: v = v + fv/dpdv c if the change from the last estimate c is less than 0.001 cm3 leave. if (dabs(fv/dpdv).lt.0.001) goto 99 itic = itic + 1 if (itic.gt.50) goto 999 goto 10 99 brvol = v return 999 write (*,*) '**error ver403** too many iterations in brvol' stop end function vdpdv (v) c----------------------------------------------------------------------- c dpdv function to calculate the finite difference approximation of c the product v(dp/dt) for the equation of bottinga and richet 1981, c ajs v 281, eq 3. c input: v - molar volume of co2 in j/bar c t - temperature in k (common block cst5) c output:vdpdv - the approximation of v(dp/dv) c j. connolly 1988. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps save rbar, dv, a1, a2, a3 data rbar,dv/83.143,0.001/a1,a2,a3/6.566d7,7.276d7,37.3/ c--------------------------------------------------------------------- rt = rbar * t t12 = dsqrt(t) if (v.le.47.22) then b1 = 1.856669 b2 = 0.0637935 else if (v.lt.180.) then b1 = 11.707864 b2 = 0.363955 else b1 = 7.352629 b2 = 0.241413 end if b = (dlog(v/a3)+b1)/b2 a = a1*( (a3/v)**3 - (a3/v)**6 ) + a2 vp = v + dv bp = (dlog(vp/a3)+b1)/b2 ap = a1*( (a3/vp)**3 - (a3/vp)**6 ) + a2 pp = rt/(v-b) - a/(v*(v+b)*t12) vdpdv = -v*(pp - (rt/(vp-bp) - ap/(vp*(vp+bp)*t12) ))/dv end subroutine simps (func,a,b,dx,ss) c------------------------------------------------------------------ c subroutine to numerically integrate the function func over the c limits a to b, in n increments by simpsons method, the integral c is returned as ss. from conte and deboor 1981. c jadc, aug 1990. c------------------------------------------------------------------ implicit double precision (a-z) integer i, n double precision func external func n = idint (dabs(b-a)/dx) if (n.lt.100) n = 100 h = (b-a)/float(n) hover2 = h/2.d0 ss = 0. half = func (a + hover2) do 1 i = 1, n-1 x = a + h * float(i) ss = ss + func(x) 1 half = half + func (x + hover2) ss = (h/6.) * (func(a) + 4.d0*half + 2.d0*ss + func(b)) end subroutine qromb (func,a,b,ss) c----------------------------------------------------------------------- c subroutine to numerically integrate a function func over the limits c a to b by romberg's method, the value of the integral is returned as c ss. c eps is the relative error acceptable for the estimate. c jmax is the number of iterations permitted. c k is the minimum number of iterations(?) c c from numerical recipes by press et al. 1988. c jadc c----------------------------------------------------------------------- implicit double precision (a-h,o-y),integer (i-m,z) parameter (eps=1.d-8, jmax = 20, jmaxp = jmax+1, k = 5) double precision s(jmaxp), h(jmaxp), func external func h(1) = 1.d0 do 11 j = 1, jmax call trapzd (func,a,b,s(j),j) if (j.ge.k) then c ???? call polint(h,s,j,0.d0,ss,dss) if (dabs(dss).lt.eps*dabs(ss)) return end if s(j+1)= s(j) 11 h(j+1)= .25 * h(j) write (*,*) '**error ver410** didnt converge in qromb' stop end c subroutine polint (h,s,k,x,ss,dss) c----------------------------------------------------------------------- c subroutine to extrapolate a polynomial of degree k - 1 to obtain its c value ss at x, where dss is an estimate of the error in ss. c the array s contains the values of the function at the corresponding c values of the array h. c c from numerical recipes by press et al. 1988. c jadc c----------------------------------------------------------------------- implicit double precision (a-h,o-y),integer (i-n,z) double precision s(k),h(k),c(40),d(40) if (k.gt.40) then write (*,*) '**error ver409** ugabugga polint k=',k stop end if ns = 1 dif = dabs(x-h(1)) do 11 i = 1, k dift = dabs(x-h(i)) if (dift.lt.dif) then ns = i dif = dift end if c(i)=s(i) 11 d(i)=s(i) ss = s(ns) ns = ns - 1 do 13 m = 1, k-1 do 12 i = 1, k-m ho = h(i) - x hp = h(i+m) - x w = c(i+1) - d(i) den = ho - hp if (den.eq.0.0) then write (*,*) '**error ver498** polint' stop end if den = w/den d(i) = hp * den 12 c(i) = ho * den if (2*ns.lt.k-m) then dss = c(ns+1) else dss = d(ns) ns = ns - 1 end if ss = ss + dss 13 continue end subroutine trapzd (func,a,b,s,j) c----------------------------------------------------------------------- c subroutine to compute the j'th stage of refinement of an extended c trapezoidal rule, returned as s, of the function func between a and b c from numerical recipes by press et al. 1988. c jadc c----------------------------------------------------------------------- implicit double precision (a-h,o-y),integer (i-m,z) double precision func external func if (j.eq.1) then s = 0.5 * (b-a)*(func(a)+func(b)) it = 1 else tnm = it del = (b-a)/tnm x = a + 0.5 * del sum = 0.0 do 11 i = 1, it sum = sum + func(x) 11 x = x + del s = 0.5*(s+(b-a)*sum/tnm) it = it * 2 end if end subroutine hprk c----------------------------------------------------------------------- c hprk routine to compute h2o-co2 fugacities from the holland and c powell (cmp 1991) c j. connolly 1992. c for the unitiated input and output is done through common blocks c cst5, cst11, and cst26, the relevant variables are: c for input: c pbars= pressure in bars c t = temperature in kelvins c r = gas constant c xco2 = mole fraction of co2 in fluid phase c for output: c v = molar volume (cm3/mole) at p and t c fco2 = the natural log of the co2 fugacity c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) common/ cst5 /pbars,t,xco2,u1,u2,tr,pr,r,ps * / cst11 /fh2o,fco2 rt = r*t/1000.d0 p = pbars/1000.d0 if (xco2.eq.1d0) then call crkco2 (pbars,t,fco2) else if (xco2.eq.0.d0) then call crkh2o (pbars,t,vol,fh2o) else call crkco2 (pbars,t,fco2) call crkh2o (pbars,t,vol,fh2o) xh2o = 1.d0 - xco2 wco2 = (13.2 - .290 * t**0.5d0) * p**0.25d0 wh2o = (7.0 - 0.15 * t**0.5d0) * p**0.25d0 gco2 = xh2o*xh2o*(wco2+2.d0*xco2*(wh2o-wco2))/rt gh2o = xco2*xco2*(wh2o+2.d0*xh2o*(wco2-wh2o))/rt fco2 = fco2 + gco2 + dlog(xco2) fh2o = fh2o + gh2o + dlog(xh2o) end if end subroutine mrk c--------------------------------------------------------------------- c routine to calculate properties of H2O-CO2 mixtures using the c MRK EoS. c--------------------------------------------------------------------- implicit double precision (a-h,p,r-y) parameter (nsp=10) integer ins(nsp), jns(nsp) common / cst5 /p,t,xc,u1,u2,tr,pr,r,ps * / cstcoh /x(nsp),g(nsp),v(nsp) save jns data jns/ 1, 2, 8*0/ if (xc.eq.1.d0) then ins(1) = 2 call mrkpur (ins, 1) goto 99 else if (xc.eq.0.d0) then ins(1) = 1 call mrkpur (ins, 1) goto 99 end if x(2) = xc x(1) = 1.d0 - xc call mrkmix (jns, 2) 99 end subroutine hsmrk c--------------------------------------------------------------------- c main program and subprograms written by gary k. jacobs june 1980 c modified for therm april 1982 jadc c--------------------------------------------------------------------- implicit double precision (a-h,p,r-y) parameter (nsp=10) integer jns(3), ins(nsp) common / cst85 /p,t,t12,r/ cst5 /pbar,tk,xc,u1,u2,tr,pr,rj,ps * / cst11 /fh2o,fco2 save bw, bc data bw, bc/ 29., 58./ if (xc.eq.1.d0) then ins(1) = 2 jns(1) = 2 call hsmrkp (ins, 1, jns, 1) goto 99 else if (xc.eq.0.) then ins(1) = 1 jns(1) = 1 call hsmrkp (ins, 1, jns, 1) goto 99 end if xw = 1.d0 - xc p = pbar t = tk t15 = dsqrt(t**3) t12 = dsqrt(t) rtt = r*t15 t2 = t*t cc = 28.31d6+0.10721d6*t-0.00000881d6*t2 dc = 9380.d6-8.53d6*t+0.001189d6*t2 ec = -368654.d6+715.9d6*t+0.1534d6*t2 cw = 290.78d6-0.30276d6*t+0.00014774d6*t2 dw = -8374.d6+19.437d6*t-0.008148d6*t2 ew = 76600.d6-133.9d6*t+0.1071d6*t2 bm = bc*xc+bw*xw cij = cc * cw dij = dc * dw eij = ec * ew if (dij.lt.0.0.or.eij.lt.0.0.or.cij.lt.0.0) then write (*,1000) p,t eij = 0.d0 dij = 0.d0 cij = 0.d0 else cij = dsqrt(cij) dij = dsqrt(dij) eij = dsqrt(eij) end if xc2 = xc*xc xw2 = xw*xw xwc2 = 2.d0*xc*xw cm = cc*xc2+cw*xw2+xwc2*cij dm = dc*xc2+dw*xw2+xwc2*dij em = ec*xc2+ew*xw2+xwc2*eij c get initial volume estimate for c newrap from mrk, implicit in newrap c solve for hsmrk volume call newrap (bm,cm,dm,em,yzm) c calculate hsmrk (log) fugacities: fco2 = dlog(xc*p)+ fug(rtt,cij,dij,eij,xc,xw,bm,yzm,cm, * dm,em,bc,cc,dc,ec) fh2o = dlog(xw*p)+ fug(rtt,cij,dij,eij,xw,xc,bm,yzm,cm, * dm,em,bw,cw,dw,ew) 1000 format (' **warning ver678** p,t (',g9.3,1x,g9.3, * ') conditions are out of range for HSMRK',/, * ' your results may be incorrect.') 99 end subroutine newrap (b,c,d,e,yz) implicit double precision (a-h,p,r-y) common / cst85 /p,t,t12,r/ cst26 /vi c call mrk to get initial volume call mrk do 10 k = 1, 50 y = b/4.d0/vi x = 1.d0 - y bi = vi + b bi2 = bi * bi vi2 = vi * vi vi3 = vi2 * vi x3 = x**3 y3 = y**3 pn = 1.d0 + y + y*y - y3 pa1 = c + d/vi + e/vi2 d1 = -.75d0 * b/vi3/x/x3 - 1.d0/vi2/x3 d3 = (-b/4.d0/vi2 - 2.d0 * b * b/16./vi3 * + 0.046875d0 * b**3/vi/vi3 )/ vi/x3 df = ( pn * d1 + d3 )*r*t * - (pa1 * (-1.d0/vi/bi2 - 1.d0/vi2/bi) * + ( -d/vi2 - 2.d0 * e/vi3)/ vi/bi )/t12 v = vi - (pn/vi/x3 * r * t - pa1/t12/vi/bi - p)/df diff = dabs(v-vi) vi = v 10 if (diff.lt.0.01) goto 99 99 yz = vi * p/83.14/t end subroutine qrkmrk c--------------------------------------------------------------------- c qrkmrk calculates the log of water and co2 fugacities using the hsmrk c for pure fluids, and the mrk for the acivities of h2o and co2 in c mixtures c----------------------------------------------------------------------- implicit double precision (a-h,p,r-y),integer (i-m,z) common / cst85 /p,t,t12,r * / cst5 /pbar,tk,xc,u1,u2,tr,pr,rcal,ps/ cst11 /fh2o,fco2 save bw, bc data bw, bc/29., 58./ c----------------------------------------------------------------------- t = tk if (xc.ge.1.d0) then call hsmrk goto 99 else if (xc.le.0.) then call hsmrk goto 99 end if xw = 1.d0 - xc p = pbar t15 = dsqrt(t**3) t12 = dsqrt(t) t2 = t*t dlp = dlog(p) rtt = r*t15 cc = 28.31d6+0.10721d6*t-0.00000881d6*t2 dc = 9380.d6-8.53d6*t+0.001189d6*t2 ec = -368654.d6+715.9d6*t+0.1534d6*t2 cw = 290.78d6-0.30276d6*t+0.00014774d6*t2 dw = -8374.d6+19.437d6*t-0.008148d6*t2 ew = 76600.d6-133.9d6*t+0.1071d6*t2 c mixed volatiles c call mrk to get ln fugacities of c h2o and co2 in binary mixes, c save as fmh2o and fmco2: call mrk fmh2o = fh2o fmco2 = fco2 xmc = xc c call mrk (this call is now c implicit in newrap) to get pure c f(h2o) & estimate V for hsmrk xw = 1.d0 xc = 0.d0 c call newrap for hsmrk h2o volume: call newrap (bw,cw,dw,ew,yzm) c calculate water fugacity save as c th2o. tfh2o = dlp + fmh2o - fh2o * + fug(rtt,cw,dw,ew,xw,xc,bw,yzm,cw,dw,ew,bw,cw,dw,ew) c co2 xc = 1.d0 xw = 0. call newrap (bc,cc,dc,ec,yzm) fco2 = dlp + fmco2 - fco2 * + fug(rtt,cc,dc,ec,xc,xw,bc,yzm,cc,dc,ec,bc,cc,dc,ec) fh2o = tfh2o xc = xmc 99 end subroutine haar c----------------------------------------------------------------------- c this code is modified from a source file obtained from c Bern Univerisity and written by Christian DeCapitani at ubc. c for a documented and more general program see c post 1984 steam tables. c haar routine to calculate molar volume and fugacity of h2o from the c equation of haar et al (steam tables). c jadc 1989. c for the unitiated input and output is done through common blocks c cst5, cst11, and cst26, the relevant variables are: c for input: c p = pressure in bars c t = temperature in kelvins c rref = gas constant c xco2 = mole fraction of co2 in fluid phase c for output: c vh2o = molar volume (cm3/mole) at p and t c fh2o = the natural log of the h2o fugacity c----------------------------------------------------------------------- implicit double precision (a-h,o-z) common/ cst26 /vh2o/ cst11 /fh2o,fco2 common/ cst5 /p,t,v(3),tref,pref,rref,v4 dimension taui(0:6), ermi(0:9), gi(40), ki(40), li(40) dimension rhoi(37:40), ttti(37:40), alpi(37:40), beti(37:40) c -----gi are in (bar cc / g) = 10 * (j / g) save ki, li, gi, rhoi, ttti, alpi, beti, r, t0, amh2o data gi /-.53062968529023d4, .22744901424408d5, .78779333020687d4, 1 -.69830527374994d3, .17863832875422d6, -.39514731563338d6, 2 .33803884280753d6, -.13855050202703d6, -.25637436613260d7, 3 .48212575981415d7, -.34183016969660d7, .12223156417448d7, 4 .11797433655832d8, -.21734810110373d8, .10829952168620d8, 5 -.25441998064049d7, -.31377774947767d8, .52911910757704d8, 6 -.13802577177877d8, -.25109914369001d7, .46561826115608d8, 7 -.72752773275387d8, .41774246148294d7, .14016358244614d8, 8 -.31555231392127d8, .47929666384584d8, .40912664781209d7, 9 -.13626369388386d8, .69625220862664d7, -.10834900096447d8, * -.22722827401688d7, .38365486000660d7, .68833257944332d5, 1 .21757245522644d6, -.26627944829770d5, -.70730418082074d6, 2 -.225d1, -1.68d1, .055d1, -93.0d1/ data ki /4*1, 4*2, 4*3, 4*4, 4*5, 4*6, 4*7, 4*9, 2*3, 1, 5, 3*2, 1 4/ data li /1, 2, 4, 6, 1, 2, 4, 6, 1, 2, 4, 6, 1, 2, 4, 6, 1, 2, 4, 1 6, 1, 2, 4, 6, 1, 2, 4, 6, 1, 2, 4, 6, 0, 3*3, 0, 2, 0, 0/ data rhoi /0.319d0, 0.310d0, 0.310d0, 1.55d0/ data ttti /640.0d0, 640.0d0, 641.6d0, 270.0d0/ data alpi /34.0d0, 40.0d0, 30.0d0, 1050.0d0/ data beti /2.0d4, 2.0d4, 4.0d4, 25.0d0/ data r, t0, amh2o/4.6152d0,647.073d0,18.0152d0/ rt = r * t nlow = 40 nhigh = 20 if (t .lt. 449.35) nhigh = 40 c -----the values (t/t0)**i are stored in the array taui(i) taui(0) = 1.d0 taui(1) = t/t0 do 10 i = 2, 6 10 taui(i) = taui(i - 1) * taui(1) b = -0.3540782 * dlog(taui(1)) + 0.7478629+0.007159876 / 1 taui(3) - 0.003528426/taui(5) bb = 1.1278334 - 0.5944001/taui(1) - 5.010996/taui(2) + 1 0.63684256/taui(4) c ps = 220.55 if (t .le. 647.25) call psat2(t, ps) c get initial guess for volume from c mrk: call mrk rhn = 1.d0/vh2o * amh2o c -----find the true(?) rh(t,p) c -----note: pr = pressure corresponding to guessed rh c dpr = (dp/drh) c the values (1-exp(-rh))**i are stored in the array ermi(i) do 50 loo = 1, 100 rh = rhn if (rh .le. 0.) rh = 1.d-8 if (rh .gt. 1.9) rh = 1.9 rh2 = rh** 2 y = rh * b/4. er = dexp(-rh) y3 = (1.d0 - y)** 3 aly = 11.d0 * y bety = 44.33333333333333 * y * y f1 = (1.d0 + aly + bety)/y3 f2 = 4.d0 * y * (bb/b - 3.5) ermi(0) = 1. ermi(1) = 1.d0 - er do 20 i = 2, 9 20 ermi(i) = ermi(i - 1) * ermi(1) pr = 0. dpr = 0. do 30 i = 1, 36 s = gi(i)/taui(li(i)) * ermi(ki(i) - 1) pr = pr + s dpr = dpr + (2.d0 + rh*(ki(i)*er - 1.d0)/ermi(1)) * s 30 continue do 40 i = nlow, nhigh del = rh/rhoi(i) - 1.d0 rhoi2 = rhoi(i) * rhoi(i) tau = t/ttti(i) - 1.d0 abc = -alpi(i) * del** ki(i) - beti(i) * tau** 2 if (abc .gt. - 100.) then q10 = gi(i) * del** li(i) * dexp(abc) else q10 = 0. end if qm = li(i)/del - ki(i) * alpi(i) * del** (ki(i) - 1) s = q10 * qm * rh2/rhoi(i) pr = pr + s dpr = dpr + s * (2.d0/rh + qm/rhoi(i)) - rh2/rhoi2 * q10 * 1 (li(i)/del/del + ki(i)*(ki(i) - 1)*alpi(i)*del**(ki(i) - 2)) 40 continue pr = rh * (rh*er*pr + rt*(f1 + f2)) dpr = rh * er * dpr + rt * ((1.d0 + 2.d0*aly + 3.d0*bety)/y3 * + 3.d0*y*f1/(1.d0 - y) + 2.d0*f2) if (dpr .le. 0.0d0) then if (p .le. ps) then rhn = rhn * 0.95d0 else rhn = rhn * 1.05d0 end if else if (dpr .lt. 0.01d0) dpr = 0.01d0 x = (p - pr)/dpr if (dabs(x) .gt. 0.1d0) x = 0.1d0 * x/dabs(x) rhn = rh + x end if dp = dabs(1.d0 - pr/p) dr = dabs(1.d0 - rhn/rh) if (dp .lt. 5.d-2 .and. dr .lt. 5.d-2) go to 60 50 continue 60 rh = rhn y = rh * b/4. x = 1.d0 - y er = dexp(-rh) ermi(0) = 1.d0 ermi(1) = 1.d0 - er do 70 i = 2, 9 70 ermi(i) = ermi(i - 1) * ermi(1) c -----calculate base function aa = rt * (-dlog(x) - 43.33333333333333/x + 28.16666666666667/ 1 x/x + 4.d0*y*(bb/b - 3.5) + 15.16666666666667+dlog(rh*rt/ 2 1.01325d0)) c -----calculate residual function do 80 i = 1, 36 80 aa = aa + gi(i)/ki(i)/taui(li(i)) * ermi(ki(i)) do 90 i = nlow, nhigh del = rh/rhoi(i) - 1.d0 tau = t/ttti(i) - 1.d0 abc = -alpi(i) * del** ki(i) - beti(i) * tau** 2 if (abc .gt. - 100.) then aa = aa + gi(i) * del** li(i) * dexp(abc) else aa = aa end if 90 continue call aideal (t/100.,rt,aid) aa = aa + aid gid = (aid * amh2o * 0.1) + rref * t gh2o = (aa + p/rh) * amh2o * 0.1 fh2o = (gh2o - gid)/rref/t c to correct to 1 bar 298 gf + c sliding scale in t ideal gas c gref should be 182051 j / mol c vh2o = amh2o/rh end c---------------------------------------------------------------- subroutine psat2(t, ps) double precision t, ps, a(8), w, wsq, v, ff save a data a /-7.8889166, 2.5514255, -6.716169, 33.239495, 1 -105.38479, 174.35319, -148.39348, 48.631602/ if (t .le. 314.00d0) then ps = dexp(6.3573118 - 8858.843/t + 607.56335/(t**0.6)) else v = t/647.25 w = dabs(1.d0 - v) wsq = dsqrt(w) ff = 0.d0 do 10 i = 1, 8 ff = ff + a(i) * w 10 w = w * wsq ps = 220.93 * dexp(ff/v) end if end c---------------------------------------------------------------- subroutine aideal (tr,rt,aid) c subroutine to compute the helmoltz free energy of ideal steam. implicit double precision (a-h, o-z) dimension ci(18) save ci data ci /.19730271018d2, .209662681977d2, -.483429455355d0, 1 .605743189245d1, 22.56023885d0, -9.87532442d0, 2 -.43135538513d1, .458155781d0, -.47754901883d-1, 3 .41238460633d-2, -.27929052852d-3, .14481695261d-4, 4 -.56473658748d-6, .16200446d-7, -.3303822796d-9, 5 .451916067368d-11,-.370734122708d-13,.137546068238d-15/ w = tr** (-3) aid = 1.d0 + (ci(1)/tr + ci(2)) * dlog(tr) do 1 i = 3, 18 aid = aid + ci(i) * w 1 w = w * tr aid = -rt * aid end subroutine trkmrk c--------------------------------------------------------------------- c trkmrk calculates the log of water and co2 fugacities using the haar c equation of state for pure water and the hsmrk for co2 and c the mrk for activities of h2o and co2 in mixtures c----------------------------------------------------------------------- implicit double precision (a-h,p,r-y),integer (i-m,z) parameter (nsp=10) integer ins(nsp) common / cst11 /fh2o, fco2 / cst111 / vhp, vcp, vmix, fhp, fcp * / cst26 / v / cst5 /p,t,xc,vv(6) data ins/ 1, 2, 8*0/ if (xc.le.0.) then c for pure h2o: xc = 0. call haar goto 99 else if (xc.ge.1.d0) then c for pure co2: call hprk goto 99 end if c mixed volatiles c call mrk to get ln fugacities of c h2o and co2 in binary mixtures, c save as fmh2o and fmco2: xmc = xc call mrk fmh2o=fh2o fmco2=fco2 c now get activities fugacities: call mrkpur (ins, 2) ah2o = fmh2o - fh2o aco2 = fmco2 - fco2 c calculate pure water fugacity call haar tfh2o = fh2o + ah2o c calculate pure CO2 fugacity xc = 1.d0 call hprk fco2 = fco2 + aco2 fh2o = tfh2o xc = xmc 99 end subroutine delany c--------------------------------------------------------------------- c delany calculates ln(fh2o) at P>10 KB using the delaney and c helgeson fit (AJS ca 81?), routine by S. Poli c----------------------------------------------------------------------- implicit double precision (a-h,p,r-y) dimension b(5,5), c(5,5), gh2o(2) common / cst5 /p,tk,xc,u1,u2,tr,pr,rj,ps/ cst11 /fh2o,fco2 save p0, b, c c data for delany & helgeson 1974 eos for h2o above 10 kbar data b/-5.6130073,-1.5285559, -2.6092451, 1.7140501, * -6.0126987, 3.8101798, 1.3752390, 3.5988857, -1.6860893, * 0.,-2.1167697,-1.5586868, -2.7916588, 0., 0., 2.0266445, * 6.6329577, 0., 0., 0.,-8.3225572,0., 0., 0., 0./ c data c/4.,1.,-2.,-5.,-9.,-1.,-4.,-8.,-11.,0.,-6.,-9., * -14.,0.,0.,-11.,-15.,0.,0.,0.,-17.,0.,0.,0.,0./ if (p.le.10000.) then call hsmrk goto 99 end if p0 = p p = 10000. call hsmrk p = p0 p0 = 10000. tc = tk - 273.15 do 90 i=1,2 90 gh2o(i) = 0. do 100 j = 1, 5 do 200 l = 1, 6 - j gh2o(1) = gh2o(1)+(b(j,l)*10**c(j,l))*(tc**(j-1))*(p**(l-1)) 200 gh2o(2) = gh2o(2)+(b(j,l)*10**c(j,l))*(tc**(j-1))*(p0**(l-1)) 100 continue fh2o = ((gh2o(1)-gh2o(2))/(1.9872d0*tk)) + fh2o 99 end subroutine dhhsrk c--------------------------------------------------------------------- c dhhsrk calculates the log of water and co2 fugacities using the haar c equation of state for pure water and the hsmrk for co2 and c the mrk for activities of h2o and co2 in mixtures c----------------------------------------------------------------------- implicit double precision (a-h,p,r-y),integer (i-m,z) parameter (nsp=10) integer ins(nsp) common / cst11 /fh2o, fco2 / cst111 / vhp, vcp, vmix, fhp, fcp * / cst26 / v / cst5 /p,t,xc,vv(6) data ins/ 1, 2, 8*0/ if (xc.le.0.) then c for pure h2o: xc = 0. call delany goto 99 else if (xc.ge.1.d0) then c for pure co2: call hsmrk goto 99 end if c mixed volatiles c call mrk to get ln fugacity c of h2o and co2 in binary mixes, c save as fmh2o and fmco2: xmc = xc call mrk fmh2o=fh2o fmco2=fco2 c now get activities fugacities: call mrkpur (ins, 2) ah2o = fmh2o - fh2o aco2 = fmco2 - fco2 c calculate pure water fugacity call delany tfh2o = fh2o + ah2o c calculate pure CO2 fugacity xc = 1.d0 call hsmrk fco2 = fco2 + aco2 fh2o = tfh2o xc = xmc 99 end subroutine dhcork c--------------------------------------------------------------------- c dhhsrk calculates the log of water and co2 fugacities using the haar c equation of state for pure water and the cork for co2 and c the mrk for activities of h2o and co2 in mixtures c----------------------------------------------------------------------- implicit double precision (a-h,p,r-y),integer (i-m,z) parameter (nsp=10) integer ins(nsp) common / cst11 /fh2o, fco2 / cst111 / vhp, vcp, vmix, fhp, fcp * / cst26 / v / cst5 /p,t,xc,vv(6) data ins/ 1, 2, 8*0/ if (xc.le.0.) then c for pure h2o: xc = 0. call delany goto 99 else if (xc.ge.1.d0) then c for pure co2: call hprk goto 99 end if c mixed volatiles c call mrk to get ln fugacity c of h2o and co2 in binary mixes, c save as fmh2o and fmco2: xmc = xc call mrk fmh2o=fh2o fmco2=fco2 c now get activities fugacities: call mrkpur (ins, 2) ah2o = fmh2o - fh2o aco2 = fmco2 - fco2 c calculate pure water fugacity call delany tfh2o = fh2o + ah2o c calculate pure CO2 fugacity xc = 1.d0 call hprk fco2 = fco2 + aco2 fh2o = tfh2o xc = xmc 99 end subroutine lohork (fo2) c---------------------------------------------------------------------- implicit double precision (a-g,k,l,o-z), integer (h,i,j,m,n) parameter (nsp=10) integer ins(nsp), jns(3) c program to calculate fh2o, fo2 for H2-H2O c mixtures. common / cst11 /fh2o,fh2/ cst26 /vol * / cst5 /p,t,xv,u1,u2,tr,pr,rj,ps * / cstcoh /xh2o,xco2,xco,xch4,xh2,x(5), * gh2o,gco2,gco,gch4,gh2,g(5),v(nsp) * / cstchx /gmh2o,gmco2,gmch4,vm(3) save ins data ins, jns / 1, 5, 8*0, 1, 2*0/ xh2 = xv xh2o = 1.d0 - xh2 c check if xo2 is <1, >0, c reset if necessary. if (xh2.gt.0.9999999999) then xh2 = 0.9999999999 else if (xh2.lt.0.0000000001) then xh2 = 0.0000000001 end if c get pure species fugacities call hsmrkp (ins, 2, jns, 1) ghh2o = gh2o/gmh2o c get mrk fugacities: call lomrk (ins, 2) c evaluate lnk's kh2o = -7.028214449 + 30607.34044/t - 475034.4632/t/t * + 50879842.55/t/t/t gh2o = ghh2o * gh2o fh2o = dlog(gh2o*p*xh2o) fh2 = dlog(gh2*p*xh2) fo2 = 2.d0 * (fh2o - fh2 - kh2o) vol = vol + xh2o * vm(1) end subroutine saxfei c--------------------------------------------------------------------- c subroutine to calculate compressibilities and fugacity coefs c according to saxena & fei (1987a, b, 1988) for species in the c c-o-h system c literature: saxena s.k. & fei y. (1987a): c fluids at crustal pressures and temperatures. contr c mineral. petrol., vol. 95, pp 370 - 375. c c saxena s.k. & fei y. (1987b): c high pressure and hig temperature fluid fugacities. c geochim. cosmochim,. acta, vol. 51, pp 783 - 791. c c saxena s.k. & fei y. (1988): c the pressure-volume-temperature equation of hydrogen c geochim. cosmochim. acta, vol. 52, pp 1195 - 1196. c--------------------------------------------------------------------- c written by peter ulmer, geophysical lab, ciw, june, 1988 c modified for vertex/frendly by j. connolly, may, 1989. c never tested. c--------------------------------------------------------------------- c z = compressibility = pv/rt c species are: 1 = co2 c 2 = co c 3 = ch4 c 4 = o2 c 5 = h2o c 6 = h2 c--------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) dimension tcr(6),pcr(6),fcoeff(6),a(3),b(2),c(3),d(2),aw(4), * bw(3),cw(4),dw(4),ah(4),bh(4),ch(4),dh(4),at(6), * bt(6),ct(6),dt(6),pr1(6),pr2(6) common/ cst5 /p,t,xco2,u1,u2,trr,prr,r,ps/ cst11 /f(2) save tcr, pcr, a, b, c, d, aw, bw, cw, dw, ah, bh, ch, dh data tcr / 304.2, 133.15, 190.7, 154.8, 647.3, 33.1/ data pcr / 73.9, 35.0, 46.4, 50.8, 221.3, 13.0/ data a / 2.0614, -2.2351, -0.39411/ data b / 5.5125e-2, 3.9344e-2/ data c / -1.8935e-6, -1.1092e-5, -2.1892e-5/ data d / 5.0527e-11, -6.3033e-21/ data aw / 1.4937, -1.8626, 0.80003, -0.3941/ data bw / 4.2410e-2, 2.4097e-2, -8.9634e-3/ data cw / -9.016e-7, -6.1345e-5, 2.2380e-5, 5.2335e-7/ data dw / -7.6707e-9, 4.1108e-8, -1.4798e-8, -6.3033e-21/ data ah / 1.6688, -2.0759, -9.6173, -0.1694/ data bh / -2.0410e-3, 7.9230e-2, 5.4295e-2, 4.0887e-4/ data ch / -2.1693e-7, 1.7406e-6, -2.1885e-4, 5.0897e-8/ data dh / -7.1635e-12, 1.6197e-10, -4.8181e-9, 0.00000/ c--------------------------------------------------------------------- c calculation of the a, b, c, d terms for the different species c calculation for co2, co, ch4, and o2 ** change to index i c the indices on lhs of equalities. if (p.lt.5000.) then write (*,*) 'this version of saxfei is invalid at ', * 'pressures less than 5000 bars.' end if c do 1010 i = 1, 4 tr = t/tcr(1) at(1) = a(1) + a(2)*tr**(-2) + a(3)*log(tr) bt(1) = b(1)/tr + b(2)*tr**(-2) ct(1) = c(1)/tr + c(2)*tr**(-2) + c(3)*tr**(-3) dt(1) = d(1)/tr + d(2)*tr**3 c1010 continue c--------------------------------------------------------------------- c calculation of a, b, c, d for h2o tr = t/tcr(5) at(5) = aw(1) + aw(2)*tr**(-2) + aw(3)*tr**(-3) + aw(4)*log(tr) bt(5) = bw(1)/tr + bw(2)*tr**(-2) + bw(3)*tr**(-3) ct(5) = cw(1)/tr + cw(2)*tr**(-2) + cw(3)*tr**(-3) + . cw(4)*log(tr) dt(5) = dw(1)/tr + dw(2)*tr**(-2) + dw(3)*tr**(-3) + . dw(4)*tr**3 c calculation of the fugacity coefficients according to c saxena & fei (1987b) for co2 pr1(1) = p/pcr(1) pr2(1) = 5000./pcr(1) fcoeff(1) = r*t * (at(1)*log(pr1(1)/pr2(1)) + . bt(1)*(pr1(1)-pr2(1)) + . ct(1)*(pr1(1)**2 - pr2(1)**2)/2.d0 + dt(1)* . (pr1(1)**3 - pr2(1)**3)/3) - r*t*log(p) . + 1000*(-40.468 + 0.06702*t + 8.3481*log(t)) fcoeff(1) = fcoeff(1)/(r*t) fcoeff(1) = exp(fcoeff(1)) if (xco2.ne.0.) then f(2) = dlog(xco2 * fcoeff(1) * p) else f(2) = 0.0 end if c calculation of the fugacity coefficients according to c saxena & fei (1987b) for h2o pr1(5) = p/pcr(5) pr2(5) = 5000./pcr(5) fcoeff(5) = r*t * (at(5)*log(pr1(5)/pr2(5)) + . bt(5)*(pr1(5)-pr2(5)) + . ct(5)*(pr1(5)**2 - pr2(5)**2)/2 + dt(5)* . (pr1(5)**3 - pr2(5)**3)/3) - r*t*log(p) . + 1000*(-130.517 + 0.06497*t + 19.4831*log(t)) fcoeff(5) = fcoeff(5)/(r*t) fcoeff(5) = exp(fcoeff(5)) if (xco2.ne.1.) then f(1) = dlog( (1.-xco2) * fcoeff(5) * p) else f(1) = 0.0 end if end subroutine browoo c--------------------------------------------------------------------- c browoo calculates the log of water and co2 fugacities using the c Brodholt & Wood, 1993 equation of state for pure water and c the hsmrk for co2 and the mrk for activities of h2o and co2 in mixtures c this and ancillary subroutines were written by c Stefano Poli, Milan, 1996. modified, but untested, July 96, JADC. c----------------------------------------------------------------------- implicit double precision (a-h,p,r-y),integer (i-m,z) parameter (nsp=10) integer ins(nsp) common / cst11 /fh2o, fco2 / cst111 / vhp, vcp, vmix, fhp, fcp * / cst26 / v / cst5 /p,t,xc,vv(6) save ins data ins/ 1, 2, 8*0/ c----------------------------------------------------------------------- if (xc.le.0.d0) then c for pure h2o: xc = 0. call browop goto 99 else if (xc.ge.1.d0) then c for pure co2: call hsmrk goto 99 end if c mixed volatiles c call mrk to get ln fugacity c of h2o and co2 in binary mixes, c save as fmh2o and fmco2: xmc = xc call mrk fmh2o=fh2o fmco2=fco2 c now get activities fugacities: call mrkpur (ins, 2) ah2o = fmh2o - fh2o aco2 = fmco2 - fco2 c calculate pure water fugacity call browop tfh2o = fh2o + ah2o c calculate pure CO2 fugacity xc = 1.d0 call hsmrk fco2 = fco2 + aco2 fh2o = tfh2o xc = xmc 99 end subroutine browop c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) parameter (nsp=10) integer ins(nsp), isp, jns(3) common / cstcoh /x(nsp),g(nsp),v(nsp)/ cst11 /fh2o,fco2 * / cst5 /p,t,xc,u1,u2,tr,pr,rj,ps save r,a1,a2,a3,a4,c1,c2,c3,c4,b1,b2,ins,jns data r,a1,a2,a3,a4,c1,c2,c3,c4,b1,b2/83.144126d0,-582468.d0, * -3038.79d0,-9.24574d-3,3.02674d9,36490.5d0,-1.02451d7, * -1.79681d8,2.18437d9,-3.90463d-2,-0.991078d0/ data ins,jns/nsp*1,1,1,1/ c----------------------------------------------------------------------- c at p < 10000 kerrick if (p.lt.10000.d0) then call hsmrk goto 9999 end if c get portion of fh2o below 10 kbar from hsmrk p0 = p p = 10000.d0 call hsmrkp (ins, 1, jns, 1) f0 = fh2o c brodholdt & wood part: c call mrk to get initial guess for volume vbw = v(1) t2 = t*t a = a1 + t*(a2 + a3*t) + a4/t2 rt = r*t t12 = dsqrt (t) t15 = t*t12 t25 = t2*t12 do 100 i = 1, 100 b = b1 + b2*vbw e0 = vbw - b e1 = b + vbw e3 = 1.d0 + b2 cor = (rt/e0+(c1-a/t12/e1+(c2+(c3+c4/vbw)/vbw)/vbw)/vbw - p) * / * ((((-4.d0*c4/vbw-3.d0*c3)/vbw-2.d0*c2)/vbw-c1)/vbw**2 * -rt*e3/e0/e0 * +(e3/e1+1.d0/vbw)/vbw/e1*(a4/t25+a2*t12+a1/t12+a3*t15)) vbw = vbw - cor 100 if (dabs(cor).lt.0.001d0) goto 99 c compute integral vdp by parts to get fugacity c vdP = V2*P2 - V1*P1 - (integral of PdV) c do integral of PdV for V1 first. 99 b = b1 + b2*vbw pdv = rt/(1.d0-b2)*dlog(vbw-b) + a/t12/b1*dlog((vbw+b)/vbw) * + c1*dlog(vbw) - (c2 + (c3/2.d0 + c4/vbw/3.d0)/vbw)/vbw c v1 is vbw at 10 kb v1 = vbw c v at p = p bar p = p0 call mrkpur (ins, isp) vbw = v(1) do 110 i = 1, 100 b = b1 + b2*vbw e0 = vbw - b e1 = b + vbw e3 = 1.d0 + b2 cor = (rt/e0+(c1-a/t12/e1+(c2+(c3+c4/vbw)/vbw)/vbw)/vbw - p) * / * ((((-4.d0*c4/vbw-3.d0*c3)/vbw-2.d0*c2)/vbw-c1)/vbw**2 * -rt*e3/e0/e0 * +(e3/e1+1.d0/vbw)/vbw/e1*(a4/t25+a2*t12+a1/t12+a3*t15)) vbw = vbw - cor 110 if (dabs(cor).lt.0.001d0) goto 999 999 b = b1 + b2*vbw fh2o = (vbw*p - v1*10000. - a/t12/b1*dlog((vbw+b)/vbw) * - c1*dlog(vbw) + (c2 + (c3/2.d0 + c4/vbw/3.d0)/vbw)/vbw * + pdv ) / rt + f0 - dlog(vbw-b)/(1.d0-b2) 9999 end subroutine hcrk c--------------------------------------------------------------------- c hcrk calculates the log of water and co2 fugacities using the c Halbach Chatterjee 82 equation of state for pure water and c the hsmrk for co2 and the mrk for activities of h2o and co2 in c mixtures. this and ancillary subroutines were written by c Stefano Poli, Milan, 1996. modified, but untested, July 96, JADC. c----------------------------------------------------------------------- implicit double precision (a-h,p,r-y),integer (i-m,z) parameter (nsp=10) integer ins(nsp) common / cst11 /fh2o, fco2 / cst111 / vhp, vcp, vmix, fhp, fcp * / cst26 / v / cst5 /p,t,xc,vv(6) save ins data ins/ 1, 2, 8*0/ c----------------------------------------------------------------------- if (xc.le.0.) then c for pure h2o: xc = 0. call hcrkp goto 99 else if (xc.ge.1.d0) then c for pure co2: call hsmrk goto 99 end if c mixed volatiles: c call mrk to get ln fugacity c of h2o and co2 in binary mixes, c save as fmh2o and fmco2: xmc = xc call mrk fmh2o=fh2o fmco2=fco2 c now get activities fugacities: call mrkpur (ins, 2) ah2o = fmh2o - fh2o aco2 = fmco2 - fco2 c calculate pure water fugacity call hcrkp tfh2o = fh2o + ah2o c calculate pure CO2 fugacity xc = 1.d0 call hsmrk fco2 = fco2 + aco2 fh2o = tfh2o xc = xmc 99 end subroutine hcrkp c-------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-m,z) common / cst11 /fh2o,fco2/ cst5 /p,t,xc,u1,u2,tr,pr,rj,ps save r, pdisc, nincr data r, pdisc, nincr/ 83.144126d0, 1000., 1000/ c-------------------------------------------------------------- c at p < pdisc kerrick if (p.lt.pdisc) then call hsmrk goto 9999 end if c get portion of fh2o below 10 kbar from MRK p0 = p p = pdisc call hsmrk p = p0 c now Halbach & Chatterjee part a = 2.d0 b = 4.d0 c Newton-Raphson numerical integration pincr = (p - pdisc) / nincr vdp = pincr * (volhc(pdisc,t) - volhc(p,t))/3.d0 aincr = pincr do 44 i = 1, nincr - 1 d = b c = a b = c a = d pres = pdisc + aincr vdp = vdp + a*pincr*volhc(pres,t)/3.d0 44 aincr = aincr + pincr fh2o = vdp/r/t + fh2o 9999 end function volhc (p,t) c-------------------------------------------------------------- implicit double precision (a-h,p,r-y) double precision ev(3) save r,a1,a2,a3,b1,b2,b3,b4,b5,b6 data r,a1,a2,a3,b1,b2,b3,b4,b5,b6/83.144126d0,1.616d8,-4.989d4, * -7.358d9,3.4505d-4,3.898d-9,-2.7756d-15,6.3944d-2,2.3776d-5, * 4.5717d-10/ c-------------------------------------------------------------- t12 = dsqrt (t) ahc = (a1 + a2*t + a3/t)/ p / t12 bhc = (1.d0 + p*(b1 + p*(b2 + b3*p))) * / (b4 + p*(b5 + b6*p)) c1 = -r*t / p c2 = ahc - r*bhc/p - bhc**2 c3 = -ahc*bhc call roots3 (c1,c2,c3,ev,vmin,vmax,iroots) if (iroots.eq.3) then volhc = vmax else volhc = ev(1) end if end