program cohscont c---------------------------------------------------------------------- c ************************ c * * c * cohsrk.6.1993 * c * * c ************************ c---------------------------------------------------------------------- c COHSRK is a fortran program to call various subroutines to calculate c C-O-H-S fluid speciation as a function of X(O), X(CO2), X(S), X(C) c a(C), f(S2), f(O2) or O2-S2-buffer assemblage. c The variables X(O), X(S), and X(C) are defined as: c X(O) = n(O)/{n(O)+n(H)} c X(S) = n(S)/{n(S)+n(C)} c X(C) = n(C)/{n(S)+n(C)+n(O)+n(H)} c where n(C), n(S), n(O) and n(H) are the total number of moles of c carbon, sulfur, oxygen and hydrogen (as opposed to the amounts c of the species) in the fluid. c The user may choose from the following routines, identified by c number: c 0 - Modified Redlich-Kwong (MRK) c 1 - Kerrick & Jacobs 1981 hard sphere MRK (HSMRK) c 2 - Hybrid MRK/HSMRK c 3 - Saxena & Fei 1987 pseudo-virial expansion c 4 - Bottinga & Richet 1981 (CO2 RK) c 5 - Holland & Powell 1990 (CORK) c 6 - Hybrid Haar et al 1979/HSMRK (TRKMRK) c 7 - Graphite buffered COH MRK fluid c 8 - Graphite buffered COH hybrid-RK fluid c 9 - Maximum X(H2O) GCOH fluid Cesare & Connolly 1993 c 10 - X(O) GCOH-fluid hybrid-MRK Connolly & Cesare 1993 c 11 - X(O) GCOH-fluid MRK Connolly and Cesare 1993 c 12 - X(O) GCOHS Connolly & Cesare 1993 c 13 - X(H) H2-H2O-hybrid c 14 - hogbrd, don't use this if you dont know what it is. c 15 - X(H) low T H2-H2O-hybrid c 16 - X(O) H-O HSMRK/MRK hybrid c 17 - X(O) H-O-S HSMRK/MRK hybrid c 18 - Delany/HSMRK/MRK hybrid, for P > 10 kb c 19 - X(O)-X(S) GCOHS Connolly & Cesare 1993 c 20 - X(O)-X(C) GCOHS Connolly & Cesare 1993 c Routines 0-3, 5-6, and 18 are for conventional P-T-X(CO2) c calculations, where X(CO2) is the mole fraction of CO2 in c a binary H2O-CO2 mixture. The Saxena & Fei routine was c programmed by someone else, and I suspect it has been c entered incorrectly. c Routine 4 is for pure CO2 fluids. c Routines 7 and 8 are for C-O-H fluids as a function of c P-T and f(O2) or f(O2)-buffer at specified c graphite activity (usually 1). I recommend choice 8. c Routine 9 returns the fugacities of O2, H2O, and CO2 c for a H:O = 2 graphite saturated C-O-H fluid, i.e., c the fugacities computed with routine 10 at X(O) = 1/3 c and a(graphite) = 1. c Routines 10 and 11 are for C-O-H fluids as a function c of X(O) at specified graphite activity. I recommend 10. c Routines 12 and 17 are for C-O-H-S and H-O-S fluids as c a function of P-T-X(O) conditions at specified sulfur c fugacity or buffer, for C-O-H-S fluids graphite activity c must also be specified. c Routines 13 and 15 are for binary H2/H2O mixtures c (effectively H-O fluids at X(O) < 1/3), routine 15 uses c parameters optimized for the solvus region. c Routine 16 is for H-O fluids as a function of P-T-X(O). c Routines 2, 8, 12, 13, 14, 16 and 17 are a hybrids of the c HSMRK EoS (Kerrick & Jacobs) and the MRK (Redlich & Kwong/ c DeSantis et al/Holloway) EoS as described by Connolly & c Cesare, J Met Geol, 11:379-388, and for most purposes I would c recommend these equations. c Routine 19 calculates C-O-H-S graphite-UNDERSATURATED fluid c properties as a function of X(O) and the atomic S/C ratio c (expressed by the variable X(S)) at specified sulfur fugacity. c This routine often may not converge. The EoS is identical to c that used in Routine 12. c Routine 20 calculates C-O-H-S graphite-UNDERSATURATED fluid c properties as a function of X(O) and the atomic fraction of c carbon in the fluid (X(C)) at specified sulfur fugacity. c This routine can be made to calculate simple COH fluid c speciation as a function of bulk composition by setting c sulfur fugacity to such a low value that the concentration of c sulfur becomes negligible. c This routine often may not converge. The EoS is identical to c that used in Routine 12. c COHSRK is primarily intended to provide an example of how the c various speciation routines in PERPLEX can be called. If you c can follow the logic of the program then you shall probably be c able to customize the program to use input, and generate output, c more in line with your specific needs. I do have another program c that generates a plot file with various species, concentrations, c and fugacities plotted against X(O), and I shall be happy to c provide it to interested users. c---------------------------------------------------------------------- c Limitations/warnings: c The C-O-H-S and H-O-S routines do not consider the species C2H6, c O2, and SO3 and pure S species, these will not be significant at c the conditions realized in crustal metamorphic environments. c It is relatively easy to incorporate new species in the routines c so if you suspect the species are important, please contact me c and I will add the species. c The HSMRK equation of state for water seems unreliable at P > c 15 Kbar, and hence all the hybrid EoS speciation routines c as well. This has not bothered me because I am mostly concerned c with crustal metamorphic conditions; however, should you wish c to make calculations at high pressure you need only replace the c calls which get pure fluid fugacities with calls to the high c pressure routine of your choice (e.g., DeSantis-Holloway MRK, c Delany & Helgeson, Bottinga & Richet, etc. etc.). c Functions for equilibrium constants are fit in the range 400-1400K. c No warnings are issued if P-T-X(O) conditions are within a solvus. c All X(O) routines fail at X(O) = 0 or 1, to avoid this problem c the routines reset extreme values of X(O) to 1d-10 or 0.9999999999 c as appropriate. c All X(O) routines may fail at X(O) = 1/3 if the fluid becomes an c essentially mono-species H2O fluid, this is most likely to occur c in simple H-O fluids. c The C-O-H-S and H-O-S routines permit calculations for a fluid c in equilibrium with pyrrhotite with a fixed atomic Fe/S ratio, c this is not the same as N(FeS), the mole fraction of troilite c relative to S2, as often reported for pyrrhotite analyses. c Routines 19 and 20 may not converge, or may converge to c invalid roots, in the vicinity of graphite saturated and c supersaturated conditions; and at very carbon under-saturated c conditions numerical slop may lead to errors on the order of c 0.4 log units in calculated properties such as f(O2). c In speciation routines, if the concentration of a species becomes c negligible (as determined by the numerical precision of the c computer being used), its properties may be undefined or c may have meaningless fluctuations. c---------------------------------------------------------------------- c NOTES ON THIS SOURCE AND COMPILING: c The MAIN program COHSRK and a BLOCK DATA are included in this file, c all the subroutines called by COHSRK are in the file flib.f. c The files may be compiled separately and the objects linked together, c or you may concatenate the two files. c This program was put together from routines used in PERPLEX (see c below) and for this reason it is not as compact or as flexible c as might be desired. The size of the flib.f may cause problems c for DOS compilers, these can be overcome by splitting the source c into two or three blocks that can then be combined during linking. c Alternatively, routines that are not considered may be useful c can be eliminated if the calls to these routines from subprogram c cfluid are also eliminated. c Specifically I would recommend eliminating: c subroutine brmrk c subroutine simps c subroutine qromb c subroutine polint c subroutine trapzd c function brvol c function vdpdv c subroutine hosrk5 c subroutine cohfit c subroutine haar c subroutine psat2 c subroutine aideal c subroutine trkmrk c subroutine saxfei c subroutine hprk c subroutine cohgra c subroutine hh2ork c subroutine lohork c subroutine lomrk c this requires elimination of the calls to: brmrk, cohfit, c trkmrk, saxfei, hprk, cohgra, hh2ork, and lohork from cfluid. c subroutines warn and error can also be eliminated if the c calls to these routines, from subroutines rfluid, brmrk, c cohgra, cohsgr and cohhyb, are replaced by statements that write c an appropriate warning or error message. c---------------------------------------------------------------------- c These routines are incorporated in the PERPLEX programs for calculating c phase equilibria and diagrams, if you would like a copy of these c programs, or if you have any problems with this program, please c contact me at: c James Connolly c IMP-ETHZ c CH-8092 Zuerich c by e-mail at: c jamie@erdw.ethz.ch c or by telephone/fax at: c 0041-1-632-3955/0041-1-632-1088 c If you have funding the cost of the PERPLEX programs is 150 US$, c otherwise I will provide the programs at no cost (although in that c case i'd appreciate it if you send 2.8 Mbytes of 3.5" diskettes with c your request). c----------------------------------------------------------------------- c I/O: most input and output is done through the common blocks below, c the significance of the variables as named in the main program is: c ifug - number indexing the requested EoS. c p - pressure, bars. c t - temperature, kelvins. c xo - X(O) for multispecies routines, and X(CO2) or X(H2) for c binary routines. c vol - molar volume for all multispecies routines, and binary c routines 0, 1, 13, 15. c fh2o - natural log (f(H2O)) for all routines. c fco2 - natural log of the species other than H2O (i.e., CO2 or H2) c in all binary routines. c fo2 - natural log (f(O2)). c fs2 - 1/2 natural log (f(S2)). c the following variables are only used for multispecies calculations: c ins(i) - pointers indicating the species are to be calculated. c isp - number of species to be calculated. c nsp - dimensioning for the maximum number of species. c ibuf - pointer indicating method of calculating f(O2) for c routines 7 and 8, or f(S2) for routines 12 and 17. c dlnfo2 - for routines 7 and 8 the displacement of the f(O2) c relative to a buffer (in log units) or the absolute c ln(f(O2)) (if ibuf = 3). For routines 12 and 17, c if ibuf = 2 then dlnfo2 is the atomic Fe/S of pyrrhotite, c if ibuf = 3 then dlnfo2 is 1/2 the natural log (f(S2)). c elag - natural log of graphite activity. c g(i) - fugacity coefficient of the ith species. c x(i) - mole fraction of the ith species. c the indices of nine species presently defined are: 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 O2 should be replaced by SO3, and ethane should be added. c----------------------------------------------------------------------- implicit double precision (a-g,o-y),integer (h-n) parameter (nsp=10,l2=5,h5=7) character specie(nsp)*3, vname*8, y*1 integer ier, igo, ins(nsp), ifug, irk double precision nc, nh, no, ns common / cst11 /fh2o,fco2/ cst5 /p,t,xo,u(6) * / cstcoh /xs(nsp),g(nsp),v(nsp)/ cst26 /vol * / cst10 /iff(2),idss(h5),ifug,ifyn,isyn * / cst100 /dlnfo2,elag,gz,gy,gx,ibuf,hu,hv,hw,hx * / cst24 /ipot,jv(l2),iv(l2) data tentoe, fo2, fs2, specie /2.302585093d0, 0.d0, 0d0, * 'H2O','CO2','CO ','CH4','H2 ','H2S','O2 ', * 'SO2','COS','UNK'/ c----------------------------------------------------------------------- open (8,file = 'cohscont.out') igo = 0 fs2 = -9999.*tentoe/2. elag = 0. c initialize species fractions/volumes do 60 i = 1, nsp 60 xs(i) = 0. vname = ' X(CO2) ' c get the users choice of EoS: call rfluid (irk, ifug, 6) c for multispecies fluids set c up species indices: if (ifug.gt.6.and.ifug.lt.13.or.ifug.ge.19) then vname = ' X(O) ' if (ifug.eq.7.or.ifug.eq.8) vname = 'log(fo2)' isp = 5 do 1 i = 1, 6 1 ins(i) = i if (ifug.ge.19) then isp = 8 ins(7) = 8 ins(8) = 9 else if (ifug.ge.12) then isp = 9 ins(7) = 7 ins(8) = 8 ins(9) = 9 end if else if (ifug.eq.16) then vname = ' X(O) ' isp = 3 ins(1) = 1 ins(2) = 5 ins(3) = 7 else if (ifug.eq.17) then vname = ' X(O) ' isp = 4 ins(1) = 1 ins(2) = 5 ins(3) = 6 ins(4) = 8 else if (ifug.eq.13.or.ifug.eq.15) then vname = ' X(H2) ' end if 40 write (*,1120) 10 if ((igo.eq.0.and.(ifug.eq.7.or.ifug.eq.8)).or. * (ibuf.ne.3.and.(ifug.eq.7.or.ifug.eq.8)).or. * ifug.eq.4.or.ifug.eq.9) then c get P-T conditions: write (*,1010) xo = 1. igo = 1 read (*,*,iostat=ier) p, t call rerror (ier, *10) else if (ifug.eq.19) then write (*,1340) read (*,*,iostat=ier) p, t, xo, elag call rerror (ier, *10) else if (ifug.eq.20) then write (*,1360) read (*,*,iostat=ier) p, t, xo, elag call rerror (ier, *10) else c or get P-T-X/f conditions: write (*,1000) vname read (*,*,iostat=ier) p, t, xo call rerror (ier, *10) if (ifug.eq.7.or.ifug.eq.8) dlnfo2 = tentoe * xo end if c if sulfur dependent, get c fs2 if user hasn't opted for c a buffer: if ( (ifug.eq.12.and.ibuf.eq.3.and.igo.ne.0).or. * (ifug.eq.17.and.ibuf.eq.3.and.igo.ne.0).or. * (ifug.ge.19.and.ibuf.eq.3.and.igo.ne.0) ) then 30 write (*,1110) read (*,*,iostat=ier) dlnfo2 call rerror (ier,*30) dlnfo2 = tentoe * dlnfo2 end if write (*,*) 'enter pmin, pmax, dp, xo' read (*,*) pmin, pmax, dp, xo write (*,*) 'enter tmin, tmax, dt' read (*,*) tmin, tmax, dt itl = (tmax - tmin) / dt + 1 ipl = (pmax - pmin) / dp + 1 write (8,*) 3+isp write (8,*) dt, dp, tmin, xmin write (8,*) itl, ipl t = tmin - dt do i = 1, itl t = t + dt p = pmin - dp do j = 1, ipl p = p + dp c call fluid routine: call cfluid (fo2, fs2) if (xs(2).eq.1.d0) then do k = 3,isp xs(ins(k)) = 0. end do xs(1) = 0. end if ns = xs(6) + xs(8) + xs(9) no = xs(1) + xs(2)*2.+ xs(3) + xs(7)*2. * + xs(8)*2. + xs(9) nc = xs(2) + xs(3) + xs(4) + xs(9) nh = (xs(1) + xs(5) + xs(6))*2. + xs(4)*4. tot = ns + no + nh + nc xox = no/(no+nh) c write the properties you want out here: write (*,2000) p,t-273.15,fo2,xox,(xs(ins(k)),k=1,isp+3) write (8,2000) p,t-273.15,fo2,xox,(xs(ins(k)),k=1,isp+3) 2000 format (20(g12.6,1x)) c output results: end do end do 1000 format (/,' Enter p(bar), T(K), ',a8': ') 1010 format (/,' Enter p(bar), T(K): ') 1110 format (/,' Enter the log10[f(S2)]:',/) 1120 format (/,' Enter a zero for pressure to quit.',/) 1130 format (/,10x,'f(H2O) = ',g12.5,/,10x,'f(CO2) = ',g12.5,/) 1140 format (/,10x,'f(CO2) = ',g12.5,/) 1150 format (/,10x,'f(H2O) = ',g12.5,/ * ,10x,'f(CO2) = ',g12.5,/ * ,10x,'log[f(O2)] = ',g12.5,/) 1160 format (/,10x,'f(H2O) = ',g12.5, * /,10x,'f(H2 ) = ',g12.5, * /,10x,'log[f(O2)] = ',g12.5,/) 1170 format (/,10x,'log[f(O2)] = ',g12.5, * /,10x,'log[f(S2)] = ',g12.5, * /,10x,' a(gph) = ',g12.5,/) 1180 format (10x,4(5x,a3,5x)) 1190 format (5x,'x',4x,4(g12.5,1x)) 1200 format (5x,'f',4x,4(g12.5,1x),/) 1210 format (/,' Change EoS, buffer, or graphite activity (y/n)?',/) 1220 format (a1) 1230 format (22x,'Speciation/Fugacities',/) 1240 format (/,22x,'Atomic Proportions/Volume',//, * 10x,'C',12x,'H',12x,'O',12x,'S',6x,'V(cm3/mol species)') 1250 format (4x,5(g12.5,1x),/) 1260 format (5x,'v',4x,4(g12.5,1x),/) 1270 format (5x,'Back-calculated X(O) = ',g16.9) 1290 format (5x,'S/H = ',g12.5,/) 1300 format (10x,'V(cm3/mol) = ',g12.5,/) 1310 format (5x,'Back-calculated X(S) = ',g16.9) 1320 format (5x,'a(gph) = ',g12.5,/) 1330 format (/,5x,'COMPOSITION IS SUPERSATURATED WITH ', * 'RESPECT TO GRAPHITE!!',/) 1340 format (/,' Enter p(bar), T(K), X(O), X(S): ') 1350 format (5x,'Back-calculated X(C) = ',g16.9) 1360 format (/,' Enter p(bar), T(K), X(O), X(C): ') 1370 format (/,5x,'Sum of species fractions: ',f14.9,/) 1380 format (/,5x,'INVALID SPECIIATION!!',/) end block data c----------------------------------------------------------------------- implicit double precision (a-g,o-y) common/ cst5 /p,t,xco2,u1,u2,tr,pr,r,ps/ cst85 /pp,tt,yy,rr data r,rr/8.3144126d0,83.14d0/ end