This is file README.COHSRK SEE ALSO EXAMPLE #15 AT WWW.PERPLEX.ETHZ.CH/PERPLEX_EXAMPLES.HTML 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 1991 (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 The EoS is identical to 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. 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 (except routine 12), and SO3 and pure S species, these will c 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 cohsrk.f, c the subroutines called by COHSRK are in the files flib.f & tlib.f. c The files may be compiled separately and the objects linked together, c or you may concatenate the 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 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-----------------------------------------------------------------------