* ******************************************************************* * * * * * PROGRAM GASWORKS * * * VERSION 1.2 * * * * * ******************************************************************* * This program evolved from program chiller, originally written by * M. Reed, U.C. Berkeley, 1976. * Chiller was further developed by N. Spycher and M. Reed * (1980-1988) then adapted to GASWORKS in late 1988 by M. Reed, * M. Hirschmann, R. Symonds, and N. Spycher * ******************************************************************* * This program needs data file GASTHERM to run * ******************************************************************* * This version designed to run with GASTHERM which contains * regression coefficients as well as log K'S. * individual log K'S are not used in GASWORKS. * ****************************************************************** * GASWORKS is used to calculate properties of overall equilibrium in * solid-liquid-gas systems upon: * 1. Heating or cooling a gas (or fluid-solid-gas mixture) at * constant composition. * 2. Heating or cooling, as in (1.), but with solid fractionatio * (bulk composition change). * 3. Mixing of two gases of different temperature and/or * composition. * 4. Reaction of a gas with a solid (e.g. rock), * causing a chemical reaction. * ****************************************************************** * ****************************************************************** * Revised: * Dec 88 M.R. reorganized real and integer statements to * eliminate a tropical cockroach ghost that munched on FO2. Also * corrected OXY and made numerous changes in main for constant * FO2 calculations * Dec 88 M.R. fixed output formats; tried new oxy form (failed, * so far); squished underflow bugs involving output of ALA. * 22 Feb 89. This version sent to Symonds * 28 Feb 89. LOOC = 3 to stop execution. MR * 25 April 89. Gave the FO2 bees their honey (changed PF to PFLUID) * and changed OXY back to the old format * Jan-Feb 90: (RS) Converted K'S to log K'S for minerals. Changed * arrays to handle 60 components. IPUNCH = 2 to plot log(moles) * for gases, otherwise log(fugacities) are plotted. * Mar-April 90: (RS) updated GASWORKS to handle the new GASTHERM. * Changes size of the mineral arrays and added a new POWERK * subroutine to calculate log K as F(T). * Oct 90 (RS) reset mintry of kicked out mineral to zero (just * after 777 statement in main and 195 statement in MINSAT) * modified 407 statement so LOOC=3 print plotting file * Nov 90 (RS) fixed the "DO 319" loop so that MTRY's of derived * species are set correctly before first iteration * Dec 90 (RS) converted GASWORKS to Fortran 77. Changed * convergence criteria. * March 91 (RS) converted GASWORKS to compile with Lahey Fortran * compiler. * May 91 (JD) adjustable dimensions for some some arrays * May 91 (MR) adjusted output to 80-column width from chiller * August 91 (RS) added GASTOT to matrix and to the WORKRUN file, * shortened the 313 to 358 zone * Oct 91 (RS) modified convergence criteria for FO2 equation * Oct 91 (MR,JD) input formats for SINC and TOTMIX changed to allow * exponential notation. added new switches and separate * routines for Symonds' and Reed's (MINPLOT) plotting programs * (see IPUNCH comment). Bufout now eliminates arrays (ALGG, * ALGM, etc.). BG3 parameter installed to adjust MAXSAQ. * Jan 92 (RS) added BG4 parameter and some documentation on BG * parameters. Corrected error in mole fraction & partial * pressure output * Mar 92 (RS) corrected error in the calculation of MTOT(H2) if FO2 * is input. Also corrected a precision-related error in * calculating AQM(NS) values (in the DO 350 loop). * May 02 (MR) change I/O unit numbers for IOUT1 and IOUT2 * 2022 May 25: Change format statements to allow lines longer * than 72 characters in Gastherm.dat. J.Palandri * ****************************************************************** PROGRAM GASWORKS IMPLICIT DOUBLE PRECISION (A-H,O-Z) *---------------------------------------------------------------------| * BG1 = parameter to scale MAXNM (minerals) | * BG2 = parameter to scale MAXNS (total species) | * BG3 = parameter to scale MAXSAQ (component species) | * MAXIN, AND MAXOX | * BG4 = parameter to scale MAXN = BG3 + BG3 + 1 | * BG5 = parameter to scale the size of the SATMIN array | * see box below for more definitions | *---------------------------------------------------------------------| INTEGER BG1,BG2,BG3,BG4,BG5 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) PARAMETER(BG5=BG3+1) DOUBLE PRECISION LGK,MTRY,MTOT,LGKM,MINTRY,LOGACT INTEGER SAQ, QAS, SAQTOT, SPECOX,OFLAG INTEGER ERLOOP, SPEC, BSTOT, SEQ, BSTP1, SPECM,SOLS,QES,BSTP2 CHARACTER*1 BLNK,SUP CHARACTER*8 AMOLE,BLANK,DUMMY,ETOIL,MIN,NAME,NOMOX,SATMIN,SUPNAM CHARACTER*8 TEXT,TYP,TITLE1,TITLE2 CHARACTER*80 BUFOUT(BG3) DIMENSION NAME(BG2),TITLE1(10),TITLE2(10),IT(BG2),DENS(BG1), 1 W(BG3),X(BG4),SUPNAM(20),FOLD(BG4), 2 FACT(BG4),CONV(BG4), SATMIN(BG5),COX(8),ISPOX(8) ** DIMENSION ACTM(20) COMMON/MAXSIZ/ MAXSAQ,MAXNS,MAXNM,MAXSTO,MAXSSL,MAXN,MAXMEM,MAXIN COMMON/PRECIS/ RMIN, RLOGMN, RLOGMX COMMON /GEORGE/ AZERO(BG2),XI, NS, NST COMMON/MORSEL/ GAMIN(BG1),B(BG1),SUMSOL(15),SOLTRY(15), 1 DGDN(BG3,20),SOLS(15,20),ISOLS(BG1),LIMSOL,JSOLS(BG1),IST(15) COMMON/MORRIS/ MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW 1 (BG1), TEMPK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1) COMMON/MORT/MTRY(BG2),GAMMA(BG2),WATM,FO2,FO2LOG,TERMLN,KK 1,OFLAG COMMON/TRIX/ TMIN(BG3) COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8) COMMON/BULL/ BSTOT,BSTP1,JGSTOT COMMON/TURNIP/ SAQ(BG3), QAS(BG3), SAQTOT COMMON/SLOUGH/ A(BG4,BG4),F(BG4), T(BG4,BG4),IV(BG4), N,ITREF COMMON/FROID/ SINC,BIT,SLIM,AKLOG,APFCL(BG2,5), 1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4), 2 XDH(4),ALK(5),LGK(BG2),TEMP,TEMPC, HEATC,COMTOT(BG3), HEATM, 3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3), 4 CVCCF(7,3), IFRAC COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,PLOG 1,GASTOT,NMTT,LOOC ** COMMON/FUGA/GASLOG(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3, ** 1 NGAS,IMIX,IDEAL,ITGAS COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH COMMON/PP/ TS(200),PS(200),NPP,INCP COMMON/BESS/ WTPC(BG3), COEFOX(BG3,8), XMW(BG3), 1 NOMOX(BG3),ITOX(BG3),NOOX, SPECOX(BG3,8),MINSOL COMMON/VAL/BVAL(20), CVAL(20) DATA ETOIL/'********'/, BLNK/' '/, BLANK/' '/ DATA SUPNAM/20*' '/, SATMIN/BG5*' '/, 1 TYP/'GRAMS '/,AMOLE/'MOLES '/,BUFOUT/BG3*' '/ DATA FACT/BG4*1./, FOLD/BG4*1.D+65/ *=====================================================================| * ASSIGNMENTS OF DEVICE NUMBERS FOR I/O OPERATIONS | * INP1 : WORKRUN (data specific to each run, including new | * data after each loop | * INP2 : GASTHERM (thermochemical data for gases, solids, and | * liquids) | * INP3 : INPUTS PRESSURE CURVE DATA. *AT PRESENT THIS FILE DOES | * NOT EXIST* THE PROGRAM DOES NOT LOOK FOR THIS FILE AND | * IT SHOULD NOT BE ASSIGNED IN THE WORKEX EXEC(IBM) OR | * .COM (VAX) FILE. | * INP4 : oxide or solid data for MINSOL option | * IPUN : plot file | * IOUT1: main output | * IOUT2: terminal remarks | * IOUT3: pickup file written after a mineral is added or | * subtracted from the matrix | *-------------------------------------------------------------------- | INP1 = 4 INP2 = 1 ** INP3 = 7 INP4 = 3 IOUT1 = 9 IOUT2 = 6 IOUT3 = 8 IPUN = 2 *---------------------------------------------------------------------| * OPEN FILES INTERNALLY | *---------------------------------------------------------------------| * OPEN(UNIT=1,FILE='gastherm.dat',STATUS='OLD') OPEN(UNIT=2,FILE='workplot.dat',STATUS='UNKNOWN') OPEN(UNIT=3,FILE='minoxgas.dat',STATUS='OLD') OPEN(UNIT=4,FILE='workrun.dat',STATUS='OLD') OPEN(UNIT=9,FILE='workout.dat',STATUS='UNKNOWN') ***** OPEN(UNIT=7,FILE='curve.dat',STATUS='OLD') OPEN(UNIT=8,FILE='workrun.sav',STATUS='UNKNOWN') *---------------------------------------------------------------------| * ARRAY MAXIMUM SIZES: | * MAXSAQ = component species | * MAXNS = derived gas species | * MAXMN = solids + liquids | * MAXTO = stoichimetric components | * MAXN = matrix size | * MAXSSL = solid solutions | * MAXMEM = solid solution end memembers | * MAXIN = maximum number of phases or phase end members | * (solids + liquids) in matrix | * MAXOX = oxides | *---------------------------------------------------------------------| MAXSAQ = BG3 MAXNS = BG2 MAXNM = BG1 MAXSTO = 8 MAXN = BG4 MAXSSL = 15 MAXMEM = 20 MAXIN = BG3 MAXOX = BG3 *---------------------------------------------------------------------| * COMPUTER PRECISION LIMITS TO AVOID OVERFLOWS, UNDERFLOWS ETC. | * For 386 PCS, smallest double precision number is about | * 10**-300 AMD ABOUT 10**300. For IBM mainframes, these are | * lower. | * RLOGMN = minimum log10 | * RLOGMX = maximum log10 | * RMIN = minimum real | *---------------------------------------------------------------------| RMIN = 1.D-300 RLOGMN = -300. RLOGMX = 300. ERLOOP= 0 MINCT= 0 NNMT = 0 NST = 0 NMT = 0 NDUMP= 0 JSTOP = 0 GRAMSO = 0. IDEAL = 0 OFLAG = 0 FO2 = 0. BIT = 0. BASE = 10. TERMLN = 1./DLOG(BASE) DO 10 I = 1,MAXOX NOMOX(I) = BLANK WTPC(I) = 0. 10 XMW(I) = 0. DO 15 I = 1,MAXSAQ TMIN(I) = 0. W(I) = 0. COMTOT(I) = 0. MTOT(I) = 0. SAQ(I) = 0 AQM(I) = 0. 15 QAS(I) = 0 DO 16 NS = 1,MAXNS DO 6 N = 1,5 6 APFCL(NS,N) = 0. DO 8 I = 1,MAXSTO SPEC(NS,I) = 0 8 COEF(NS,I) = 0. NAME(NS) = BLANK AZERO(NS) = 0. IT(NS) = 0 MTRY(NS) = 0. LGK(NS) = 0. ACT(NS) = 0. 16 GAMMA(NS) = 0. DO 17 I = 1,MAXMEM CVAL(I) = 0. 17 BVAL(I) = 0. DO 27 NM = 1, MAXNM DO 23 N = 1,5 23 BPFCL(NM,N) = 0. DENS(NM) = 0. SCALE(NM) = 0. GAMIN(NM) = 1.0 QES(NM) = 0 LGKM(NM) = 0. SEQ(NM) = 0 B(NM) = 0. ISOLS(NM) = 0 ITM(NM) = 0 XMOLW(NM) = 0. MIN(NM) = BLANK JSOLS(NM) = 0 27 MINTRY(NM) = 0. READ(INP1,870) TITLE1 READ(INP1,870) TITLE2 *---------------------------------------------------------------------| * TITLE IN TWO, 80-CHARACTER LINES | *---------------------------------------------------------------------| READ(INP1,816) READ(INP1,834) ERPC, FO2LOG,PFLUID, TEMP, SINC,SLIM,GASTOT READ(INP1,816) READ(INP1,836) TEMPC,TOTMIX,SOLMIN,RM,GASGRM, SUPRNT IF (FO2LOG.NE.0.) FO2 = 10.**FO2LOG *---------------------------------------------------------------------| * ERPC IS THE FRACTIONAL CONVERGENCE LIMIT ON IONIC STRENGTH. * IF IT IS NOT SPECIFICALLY INPUT A VALUE OF 1D-12 IS USED * IF FO2LOG EQUALS NOT ZERO, FO2 IS FIXED AT SPECIFIED VALUE. USUALL * FO2LOG IS SET TO ZERO, IN WHICH CASE THE SPECIFIED MTOT(H2) * WILL BE USED IN H2 MASS BALANCE EQUATION AND FO2 IS CALCULATED. * THE INITIAL MTOT(H2) IS USUALLY DRAWN FROM AN EARLIER SOLVGAS * RUN ON THE DESIRED WATER COMPOSITION. * PFLUID , IN ATM, IS THE INITIAL PRESSURE OF THE SYSTEM. * IF USING WITH NON IDEAL GASES, CHECK FOR PRESSURE LIMITS * FOR COMPUTATION OF FUGACITY COEFFICIENTS (SUBROUTINE PHI). * SLIM IS THE LAST STEP AFTER WHICH * WE WANT TO STOP THE CALCULATION. DEFAULT IS 0 C * IF COOLING, 1227 C IF HEATING. * TEMP IS THE INITIAL TEMPERATURE FOR THE UNMIXED/UNCOOLED * SOLUTION. INPUT IN DEGREES C. * SINC CAN BE (+) OR (-), AND IS THE STEP INCREMENT (T, P, GRAMS, OR * MIX FRACTION). TYPICALLY SINC = -16 DEG.C IS A WORKABLE * MAXIMUM STEP. * IF USING GASMIX, SINC IS THE FRACTIONAL * TITRATION FACTOR, E.G. .01. THE CHANGE OF MTOTS IS * SINC * COMTOT, AS CALCULATED IN MIXER. * SINC CANNOT BE ZERO! IF ZERO, IT IS DEFAULTED TO -5 FOR * T CHANGE, OR 0.1 IF MIXING * TEMPC IS THE TEMPERATURE OF THE PHASE TO BE TITRATED IN * (USUALLY THE COOLER ONE, BUT NOT NECESSARILY) WHEN USING * GASMIX. IF NOT USING GASMIX, SET TEMPC = 0 !!, THIS IS THE * TEST FOR ACTIVATING THE GASMIX OPTION. * TOTMIX IS THE TOTAL FRACTION OF MIXING GAS OR OTHER SUBSTANCE * (E.G. GAS, ROCK) ADDED WHEN USING MIXING (TEMPC OR MINSOL * NON ZERO). FOR AN INITIAL RUN,TOTMIX=0; THEREAFTER, TOTMIX * IS SIMPLY THE SUM OF ALL PREVIOUS * MIXING INCREMENTS (SINCS). NOTICE THAT TOTMIX (AND SINC, WHEN * USING THE MIXING OPTION, BECAUSE TOTMIX IS A SUMMATION OF * SINCS) IS SIMPLY A DIMENSIONLESS NUMBER THAT DETERMINES THE * TITRATION STEP SIZE FOR GAS-GAS MIXING OR ROCK-GAS MIXING. * THE ACTUAL QUANTITY OF TITRATED MATERIAL (IN MOLES OR GRAMS) IS * FIXED BY THE VALUES OF COMTOT(I) (SEE BELOW) USED. * SOLMIN MINIMUM MINTRY VALUE ALLOWED FOR GASES OR MINERALS * DURING ITERATION. IF COMPUTED MINTRY GETS SMALLER DURING * ITERATIONS, THE MINERAL OR GAS WILL BE KICKED OUT OF THE * MATRIX DURING THE ITERATION PROCESS (SAVES TIME, EASES CONVER- * GENCE). MUST BE 0 OR NEGATIVE. DURING ITERATION, A NEGATIVE * MINTRY GENERALLY MEANS THE PHASE DOES NOT BELONG ANYMORE. THIS I * NOT ALWAYS THE CASE HOWEVER, AND A TOO LARGE SOLMIN MIGHT RESULT * IN SWAPPING A PHASE IN, AND OUT OF THE MATRIX INDEFINITELY. IF * SET TO A LARGE NEGATIVE VALUE, EQUILIBRATION WITH NEGATIVE MINTRY * IS ALLOWED, REQUIRING ONE MORE EQUILIBRATION STEP TO REMOVE * THE PHASE WITH NEGATIVE MINTRY (OLD METHOD). SOLID SOLUTIONS ARE * KICKED OUT IF MINTRY IS LESS THAN ZERO. * SET SOLMIN > 0 TO SKIP THIS OPTION (EQUILIBRATION W/ NEGATIVE * AMOUNTS ALLOWED); SOLMIN > 0 WILL KEEP SOLID SOLUTIONS IN * MATRIX, HOWEVER. IN THE EVENT ONE SOL-SOL SHOULD UNDERSATURATE, * CONVERGENCE WILL FAIL W/ SOLMIN > 0. * RM IS NUMBER BETWEEN 0 AND 1 TO MULTIPLY F'S TO HELP CONVERGENCE. * SET TO 0 UNLESS THERE ARE CONVERGENCE PROBLEMS. * GASGRM GRAMS OF STARTING GAS PHASE, TO COMPUTE WATER/GAS * RATIO IF MINSOL OPTION IS SELECTED * SUPRNT DERIVED SPECIES MINIMUM MOLE NUMBER FOR OUTPUT. ONLY * SPECIES WITH MOLE NUMBER GREATER OR EQUAL TO SUPRNT WILL PRINT * DEFAULT IS 1.D-10 *---------------------------------------------------------------------| IF(SUPRNT.EQ.0.) SUPRNT = 1.D-10 *---------------------------------------------------------------------| * PMAX IS USED IN CHILLER TO SET UPPER PRESSURE LIMITS FOR TREATMENT * OF NON-IDEAL GASES. NOT IN USE HERE SO SET TO 10 KB. *---------------------------------------------------------------------| PMAX = 10000. IF(PFLUID.EQ.0.) PFLUID = 1. IF(ERPC.EQ.0.) ERPC = 1.0D-12 READ(INP1,816) * COLUMN 5 10 15 20 25 READ(INP1,800) BSTOT, IFRAC, IPUNCH, NLOOP, ISTEP 1, LIMSOL, LOOC, ITREF, IDEAL, INCREM, INCP, MINSOL, OFLAG * COL 30 35 40 45 50 55 60 65 IF(IPUNCH.EQ.1.OR.IPUNCH.EQ.2) THEN WRITE(IPUN,870) TITLE1 WRITE(IPUN,870) TITLE2 END IF *---------------------------------------------------------------------| * BSTOT IS THE NUMBER OF COMPONENT SPECIES IN GASTHERM (E.G. 35). * IFRAC = 0 RESULTS IN NO FRACTIONATION DURING COOLING * IFRAC = 1 RESULTS IN FRACTIONATION OF MINERALS. * IPUNCH NOT ZERO RESULTS IN CREATING A PLOT FILE, BELOW. * There are two basic plotting options controlled by the value of * IPUNCH. Each option has two internal variations, thus four * non-zero values of IPUNCH are used to control plotting, * as follows: * IPUNCH Routine Effect * ------ ------- ------ * 1 Symonds PLOT FILE WILL HAVE LOG(FUGACITIES) * FOR GAS SPECIES * 2 Symonds PLOT FILE WILL CONTAIN LOG (MOLES) * 3 Minplot Use 3 only for very first equilibration * of very first run. It puts the header * on the plot file, then IPUNCH is reset * to 4. * 4 Minplot Use for all "pickup" runs, after the * first plot file segment has been * written. * NLOOP = 0 NO CALCULATIONS ARE PERFORMED. PROGRAM STOPS AFTER * ECHOING ALL INPUT DATA, ALL SELECTED INPUT PARAMETERS, AND AN * INDEXED LIST OF ALL GAS SPECIES, AND MINERALS. USEFUL TO * GET MINERAL INDEXES (SEQ), AND CHECK ALL INPUT DATA. * NLOOP > 0 LIMIT ON THE NUMBER OF NEWTON ITERATIONS ALLOWED * BEFORE ARBITRARILY STOPPING. 70 IS A GOOD LIMIT. USUALLY 8-20 * LOOPS ARE USED. SOMETIMES MORE THAN 70 ARE NEEDED, EVEN FOR * A SOLVABLE CASE, E.G., 120 MIGHT SOMETIMES BE NEEDED AND MAY BE * TRIED ON OCCASSION TO WORK THROUGH A COMPUTATIONAL KNOT. * ISTEP NOT ZERO RESULTS IN PRINTING OF EXTRA OUTPUT FOR CLOSER * FOLLOWING OF CALCULATION PROGRESS (USE ONLY FOR DEBUGGING). * LIMSOL NOT ZERO ALLOWS INCLUSION OF SOLID SOLUTIONS (OR MIXED * GASES). IF LIMSOL = 0, THE AQUEOUS PHASE ONLY IS OF VARIABLE * COMPOSITION (THEREFORE ENDMEMBERS ARE TREATED AS REGULAR PURE * MINERALS). * LIMSOL = 2 RESULTS IN KICKING OUT OF THE MATRIX ENDMEMBERS * WITH CALCULATED AMOUNTS LESS THAN 0. THIS IS THE ONLY WAY TO * THROW OUT A SOLSOL THAT BECOMES UNDERSATURATED, BUT MIGHT * RESULT IN KICKING OUT SOLSOLS THAT ARE STILL SATURATED. * LIMSOL = 1 WILL KEEP SOLSOL IN, AND CONVERGENCE WILL FAIL IF * A SOLSOL BECOME UNDERSAT. * LOOC DETERMINES MINERAL SELECTION METHOD IF MORE THAN ONE * MINERAL SUPERSATURATE. LOOC = 0 RESULTS IN SELECTION OF MOST * SUPERSATURATED MINERAL. LOOC NOT 0 AND NOT 3 * RESULTS IN BACK STEPPING * (REVERSE T INCREMENT OR "BACK-MIX") TO THE POINT WHERE ONLY ONE * MINERAL SUPERSATURATES. * LOOC = 3 CAUSES EXECUTION TO STOP AFTER THE FIRST EQUILIBRATION. * ITREF = 0 NO ITERATIVE REFINEMENT IS USED WHEN SOLVING * THE EQUATIONS. THIS SHOULD ALMOST ALWAYS WORK. * = 1 NO ITERATIVE REFINEMENT IS USED, BUT THE * RESIDUAL ERRORS ARE CALCULATED. ISTEP SHOULD * BE SET TO 1 TO PRINT THESE VALUES. * > 1 THE RESIDUAL ERRORS ARE CALCULATED AND MINIMISED * BY ONE OR MORE ITERATIVE REFINEMENTS. THIS * OPTION SHOULD BE USED ONLY IN CASES WHEN CONVER- * GENCE DOES NOT OCCUR DUE TO MACHINE PRECISION. * IDEAL = 0 NON IDEAL GASES AND NON IDEAL MIXING ASSUMED * = 1 NON IDEAL GASES AND IDEAL MIXING ASSUMED * = 2 IDEAL GASES AND IDEAL MIXING ASSUMED * INCREM NOT ZERO PRODUCES IMMEDIATE TITRATION OF SINC-SIZED CHUNK * OF REACTANT, BEFORE ANY EQUILIBRATIONS. RELEVANT ONLY * IF TEMPC IS NON ZERO. USE TO JUMP AHEAD IN TITRATION. * CANNOT BE USED WITH FRACTIONATION, AND WILL RESET IFRAC TO 0 * INCP = 1: P DECREASE AT CONSTANT T. USE SINC AND SLIM WITH * P VALUES INSTEAD OF T. THIS OPTION IS OVERWRITTEN BY OTHER P * VARIATIONS OPTIONS (SEE COOLER). * MINSOL: USED TO ENABLE SOLID TITRATION FOR WATER/ROCK REACTIONS * TEMPC MUST BE SET TO ZERO FOR THIS OTION * = 0 DISABLED * = 1 TITRATE SINC AMOUNT OF SOLIDS IN GRAMS * = 2 SAME AS 1, BUT STOPS TITRATING A MINERAL IF EQUILIBRATED * = 3 TITRATE SINC AMOUNT OF SOLIDS IN MOLES * = 4 SAME AS 3, BUT STOPS TITRATING A MINERAL IF EQUILIBRATED * * OFLAG: DETERMINES WHETHER FO2 IS ALLOWED TO CHANGE AFTER FIRST * EQUILIBRATION. * >0 FO2 DOES NOT CHANGE * =0 FO2 IS CALCULATED AFTER EACH TEMPERATURE CHANGE OR * TITRATION * THERE ARE 4 POSSIBILITIES: * 1. OFLAG> 0 AND FO2 IS GIVEN IN WORKRUN. * 2. OFLAG> 0 AND FO2 IS NOT GIVEN IN WORKRUN. * 3. OFLAG= 0 AND FO2 IS GIVEN IN WORKRUN. * 4. OFLAG= 0 AND FO2 IS NOT GIVEN IN WORKRUN. * THE RESULTS OF WHICH WILL BE: * 1. FO2 IS HELD AT GIVEN VALUE THROUGHOUT RUN. * 2. FO2 IS CALCULATED DURING FIRST EQUILIBRATON AND HELD * AT THAT VALUE FOR THE REST OF THE RUN. * 3. FO2 IS AT GIVEN VALUE DURING FIRST EQUILIBRATION, BUT * ALLOWED TO VARY DURING REST OF RUN. * 4. FO2 IS ALWAYS RECALCULATED. *----------------------------------------------------------------------| IF(SINC.NE.0.) GO TO 22 SINC = -5. IF(TEMPC.NE.0..OR.MINSOL.NE.0) SINC = 0.1 22 IF(SLIM.NE.0.) GO TO 25 IF(MINSOL.EQ.0.AND.TEMPC.EQ.0..AND.SINC.LT.0.) SLIM = 25. IF(MINSOL.EQ.0.AND.TEMPC.EQ.0..AND.SINC.GT.0.) SLIM = 1227. 25 CONTINUE IF(MINSOL.GT.2) TYP = AMOLE IF(INCREM.GT.0) IFRAC = 0 IMINT=BSTOT LIMSS=LIMSOL * TEMPC NON ZERO GASMIX OPTION ON, CANNOT HAVE MINSOL ON TOO IF(TEMPC.NE.0.) MINSOL = 0 24 CONTINUE * NOW WE READ DATA RELATED TO THE CHEMICAL COMPOSITION OF THE SYSTEM READ(INP1,816) I = 1 100 READ(INP1,803) SAQ(I), MTOT(I), MTRY(I), GAMMA(I), COMTOT(I) * SAQ SPECIFIES THE SEQUENCE POSITIONS OF THE COMPONENT SPECIES * TO BE DRAWN FROM GASTHERM FOR THIS CALCULATION. E.G., IF YOU * WANT H2,H2O,CO2,H2S,HCL,HF,HBR,NACL,KCL,FECL3,CD,ALF3,CUCL2,ETC. * PUT IN 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, ETC. * THESE NUMBERS ARE PICKED FROM THE ORDER OF COMPONENT * SPECIES LISTED IN SOLTHERM. H2 AND H2O (1, AND 2) MUST BE * ENTERED FIRST, IN THAT ORDER FOR PRESENT VERSION OF GASTHERM. * OTHER SPECIES MAY BE ENTERED IN ANY ORDER. * MTOT IS THE TOTAL NUMBER OF MOLES OF EACH BASIS SPECIES IN THE * STARTING GAS OR IN STARTING (GAS PLUS SOLIDS). * MTRY (I) CONTAINS TRIAL MOLE NUMBERS FOR COMPONENT SPECIES. USUALL * THESE ARE KNOWN FROM AN EARLIER SOLVGAS RUN. * GAMMA ARE FUGACITY COEFFICIENTS. IF 0 THEY ARE RESET TO 1. * IN SOME CASES WHERE NUMERICAL STABILITY IS A PROBLEM, GAMMA(I) * NOT EQUAL TO 1 MAY BE NEEDED FOR INITIAL VALUES. * COMTOT CONTAINS TOTAL MOLES OF COMPONENT SPECIES IN MIXER * PHASE (THE ONE TO BE TITRATED IN), FOR USE IN GASMIX. ** COMTOT MAY ALTERNATIVELY BE USED TO SPECIFY COMPOSITION OF ANY ** REACTANT TO BE TITRATED IN. E.G. TO REACT SOLUTION WITH ONE ** KILOGRAM OF DOLOMITE, SPECIFY ** COMTOT (1) = -10.8454 (NEGATIVE 10.8454 MOLES OF H+ IN ONE ** KILOGRAM OF DOLOMITE) ** " (2) = 0 (NO H20 IN DOLOMITE) ** " (3) = 0 (NO CL- IN DOLOMITE) ** " (4) = 0 (ETC. FOR ZERO VALULES OF COMTOT(I)) ** " (5) = 10.8454 (POSITIVE 10.8454 MOLES OF BICARBONATE IN ** ONE KILOGRAM OF DOLOMITE) ** " (6) = 0 ** " (7) = 0 ** " (8) = 0 ** " (9) = 5.4227 (POSITIVE 5.4227 MOLES OF CALCIUM IN ONE ** KILOGRAM OF DOLOMITE) ** " (10) = 5.4227 (POSITIVE 5.4227 MOLES OF MAGNESIUM IN ONE ** KILOGRAM OF DOLOMITE) ** " (11 - 29) = 0 ** I.E. 1 KG. OF DOLOMITE CONTAINS -10.8454 MOLES OF H+, +10.454 ** MOLES OF HCO3-, +5.4227 MOLES OF CA++ AND MG++. THE LIMSOL ** OPTION, HOWEVER, IS A LOT EASIER TO USE FOR THIS TYPE OF ** REACTION. * READS UNTIL BLANK LINE IF(SAQ(I).EQ.0) GO TO 144 IF(SAQ(I).GT.MAXSAQ) GO TO 143 J = SAQ(I) QAS(J) = I * SET DEFAULT VALUES IF(GAMMA(I).LE.0.) GAMMA(I) = 1. IF(MTRY(I).LE.0.) MTRY(I) = 1.D-10 I = I + 1 IF(I.LE.MAXSAQ) GO TO 100 READ(INP1,805) IDUM * (IN CASE WE HAVE EXACTLY MAXSAQ COMPONENTS) IF(IDUM.EQ.0) GO TO 144 WRITE(IOUT1,630) MAXSAQ WRITE(IOUT2,630) MAXSAQ 630 FORMAT(//5X,'WARNING! THE NUMBER OF COMPONENTS GIVEN IN WORKRUN 1 EXCEEDS THE'/5X,'CURRENT ARRAY DIMENSION (',I2,'). CHECK', 2 ' WORKRUN DATA FILE.'/) STOP 143 WRITE(IOUT1,631) I, SAQ(I), MAXSAQ WRITE(IOUT2,631) I, SAQ(I), MAXSAQ 631 FORMAT(//5X,'WARNING! SAQ(',I2,') = ',I4,'; MAXIMUM VALUE' 1 ,' ALLOWED IS ', I3/) STOP ************************************** * THE DO 123 LOOP THAT FOLLOWS IS TO ESTIMATE THE INITIAL NUMBER * OF MOLES OF GAS. GASTOT IS THE TOTAL NUMBER OF MOLES OF GAS. 144 SAQTOT = I - 1 IF(GASTOT.NE.0.) GO TO 124 GASTOT=0.0 DO 123 I=1, SAQTOT 123 GASTOT = GASTOT +MTOT(I) 124 GSTOLD = GASTOT ************************************** TEMPK = TEMP + 273.15 * READS EQUILIBRIUM ASSEMBLAGE (MINERALS AND OR GASES) AND AMOUNTS * ACT IS USED HERE AS A DUMMY INPUT VARIALBLE FOR MINTRY, AFTER * STATEMENT 165 MINTRY IS SET = ACT * MINTRY CONTAINS THE NUMBER OF MOLES OF MINERALS CURRENTLY IN * EQUILIBRIUM AT P, T, AND X. IF THIS IS THE INITIAL RUN, NORMALLY * NO MINERALS ARE KNOWN TO BE IN EQUILIBRIUM. IN ANY CASE, EVEN * IF CERTAIN MINERALS ARE WRONGLY INCLUDED, GASWORKS WILL THROW * THEM OUT. READ(INP1,817) I = 1 153 READ(INP1,823) SATMIN(I), ACT(I) IF(ACT(I).EQ.0.) ACT(I) = 1.D-9 IF(SATMIN(I).EQ.BLANK) GO TO 154 I = I + 1 IF(I.LE.MAXIN) GO TO 153 READ(INP1,823) DUMMY IF(DUMMY.EQ.BLANK) GO TO 154 WRITE(IOUT1,652) MAXIN WRITE(IOUT2,652) MAXIN 652 FORMAT(/5X,'TOO MANY SATURATED MINERALS SPECIFIED IN WORKRUN.' $ ,' LIMIT =' ,I3//) STOP 154 NMT= I-1 * READS NAMES AND AMOUNTS OF SOLID REACTANT FOR MINSOL OPTION. * MAXOX SOLIDS MAXIMUM. IF NONE, INPUT 3 BLANK LINES. * WTPC IN WEIGHT % IF MINSOL IS 1 OR 2, IN MOLES % IF MINSOL IS 3,4 I = 1 READ(INP1,817) 133 READ(INP1,823) NOMOX(I), WTPC(I) IF(NOMOX(I).EQ.BLANK) GO TO 132 I = I + 1 IF(I.LE.MAXOX) GO TO 133 READ(INP1,823) DUMMY IF(DUMMY.EQ.BLANK) GO TO 132 WRITE(IOUT1,654) MAXOX WRITE(IOUT2,654) MAXOX 654 FORMAT(/5X,'TOO MANY SOLID REACTANTS SPECIFIED IN WORKRUN.' $ ,' LIMIT =' ,I3//) STOP 132 NOOX = I-1 * NOW WE READ NAMES OF SPECIES, GASES OR MINERALS TO BE SUPPRESSED * READS UNTIL BLANK OR EOF READ(INP1,817) DO 70 NSUP = 1,20 READ(INP1,822,END=72) SUPNAM(NSUP) 70 IF(SUPNAM(NSUP).EQ.BLANK) GO TO 72 72 IF(NSUP.NE.20) NSUP = NSUP-1 * IF THE MINSOL OPTION IS SELECTED, THE OXIDE OR SOLID DATA * FILE IS READ BELOW, AND THE NAMES OF SOLIDS GIVEN ABOVE ARE * MATCHED. THE COMPOSITIONS OF NOMOX (REACTANTS) ARE ARBITRARY. USE * MINERALS, OXIDES, ETC. WILL READ UNTIL ALL NAMES MATCH OR EOF. * *** WARNING: A REACTANT CANNOT BE ENTERED TWICE IN THE OXIDE DATA * FILE. IF IT IS, THIS MIGHT RESULT IN SKIPPING READING OF FURTHER * OXIDE DATA (MESSES UP ICOUNT !). IF(MINSOL.EQ.0) GO TO 161 ICOUNT = 0 152 READ(INP4,854,END=156) DUMMY,SUP,WTX,ITOT, 1 (COX(J),ISPOX(J), J = 1, ITOT) IF(ITOT.LE.MAXSTO) GO TO 150 WRITE(IOUT1,659) DUMMY, ITOT, MAXSTO 659 FORMAT(/5X,'SOLID ',A8,' IN MINOX HAS TOO MANY COMPONENTS' $ /5X,'ITOT FOR THAT SOLID IS ',I3,' ; MAXIMUM IS ',I3//) STOP 150 CONTINUE DO 155 I = 1,NOOX 155 IF(NOMOX(I).EQ.DUMMY) GO TO 158 GO TO 152 * IF EOF REACHED, WE CHECK WHICH SOLIDS WERE NOT FOUND IN MINOX 156 CONTINUE DO 157 I = 1,NOOX IF(XMW(I).NE.0.) GO TO 157 WRITE(IOUT1,658) NOMOX(I) WRITE(IOUT2,658) NOMOX(I) 658 FORMAT(/5X,'SOLID: ',A8,' NOT FOUND IN MINOX DATA FILE'/) 157 CONTINUE STOP * 158 ICOUNT = ICOUNT + 1 XMW(I) = WTX ITOX(I) = ITOT DO 159 J = 1, ITOT COEFOX(I,J) = COX(J) JS = ISPOX(J) IF(QAS(JS).GT.0) GO TO 159 WRITE(IOUT1,661) NOMOX(I), JS WRITE(IOUT2,661) NOMOX(I), JS 661 FORMAT(/5X,'REACTANT ',A8, ' CONTAINS SPECIES NUMBER ',I2, 1 ' WHICH IS NOT',/5X,'INCLUDED IN THE LIST OF COMPONENT SPECIES', 2 ' FOR THIS RUN.'//) STOP 159 SPECOX(I,J) = QAS(JS) IF(ICOUNT.LT.NOOX) GO TO 152 161 NPP=0 ** IF THERE IS A FILE CONTAINING A P-T CURVE, TO BE USED FOR A ** BOILING PATH, PFLUID WILL BE CALCULATED FROM THE ** PRESSURE CURVE GIVEN IN INP3. ** THE CORRESPONDING PRESSURE WILL BE INTERPOLATED BETWEEN THE ** TWO CLOSEST GIVEN T-P POINTS. * ** AT PRESENT THERE IS NO SUCH INPUT FILE FOR GASWORKS ** * THERE IS ANOTHER WAY TO GET A P-T CURVE AFTER STATEMENT 160 ** DO 60 I=1,200 ** READ(INP3,835,ERR=65) TS(I),PS(I) * ********************* ** NPP=NPP+1 ** 60 CONTINUE 65 CONTINUE BSTP1 = BSTOT + 1 * ********************************************************************* WRITE(IOUT1,999) * * STARTS READING GASTHERM DATA FILE. INDIVIDUAL LOG K'S ARE * SKIPPED (SECOND ENTRY FOR GASES AND MINERALS, * FIRST ENTRY FOR OTHERS.) THESE ARE READ AS DUM OR SKIPPED * BY A SLASH IN READ FORMATS. * THE FOLLOWING SKIPS THE COMMENTS AT THE TOP OF THE GASTHERM * DATA FILE UNTIL A LINE OF STARS IS REACHED (8 STARS MINIMUM, * FROM 1ST COL.) 62 READ(INP2,822) TEXT IF(TEXT.NE.ETOIL) GO TO 62 ** DO 160 J = 1, 9 ** IF(J.LT.9) READ(INP2,822) DUM ** READ(INP2,817) (XPFC(J,I),I=1,5) ** DO 162 I = 1, 5 **162 COE(I) = XPFC(J,I) ** CALL POWER(COE,TEMP,AKLOG) * ********** ** XF(J) = AKLOG **160 CONTINUE * THE LINES BELOW FORTRAN LINE 62 ARE PRESENTLY DISABLED, BUT CAN * BE USED IN THE FUTURE TO FILL ARRAY XF FROM GASTHERM. * THE ARRAY XF CONTAINS 11 ENTRIES WHICH ARE USED IN *CHILLER* * PRIMARILY FOR DEBYE-HUCKLE PARAMETERS AND TO SPECIFY GAMMAS FOR * DISSOLVED CO2. THESE ARE UNUSED IN THE PRESENT VERSION OF * *GASWORKS*, BUT MAY IN THE FUTURE BE USED FOR NON-IDEAL TREATMENT * OF GASES. XF(9) IS DESIGNATED FOR A POWER FUNCTION GIVING * P(FLUID) AS A FUNCTION OF T, THOUGH NONE IS PRESENTLY INSTALLED. * ALL DATA DESTINED FOR XF SHOULD BE ADDED TO GASTHERM BELOW THE * LINE OF STARS BUT BELOW THE LIST OF COMPONENT SPECIES. ** IF(XF(9).NE.0.) PFLUID=XF(9) * IF(NPP.NE.0) CALL PRESS(PFLUID,TEMP) PLOG = DLOG10(PFLUID) DO 110 NJ = 1, BSTOT * * DO 110 LOOP READS COMPONENT SPECIES FROM GASTHERM, SKIPPING THOSE * NOT SPECIFIED BY SAQ, ABOVE. IF(QAS(NJ).NE.0)GO TO 103 READ(INP2,820) GO TO 110 103 CONTINUE NS = QAS(NJ) READ(INP2,822) NAME(NS), W(NS) * W CONTAINS THE MOLECULAR WEIGHT OF THE BASIS SPECIES 110 CONTINUE NS = SAQTOT * READS PROPERTIES OF DERIVED SPECIES FROM GASTHERM. NAME IS * AS ABOVE. 113 CONTINUE NS = NS + 1 IF(NS.LE.MAXNS) GO TO 114 NS = NS-1 WRITE(IOUT1,660) MAXNS, NAME(NS), NS 660 FORMAT(/5X,'NUMBER OF SPECIES READ IN GASTHERM EXCEEDS MAXIMUM' $ ,' OF ',I3, /5X,'LAST NAME READ IS: ',A8,' (NS=',I3,')'//) STOP 114 READ(INP2,820) NAME(NS),SUP,ITOT,(COEF(NS,I), 1 SPEC(NS,I), I = 1,ITOT) * SUP- SEE COMMENTS BELOW, SUP = NOT BLANK SUPPRESSES SPECIES. * ITOT CONTAINS NUMBER OF COMPONENT SPECIES IN THE DERIVED SPECIES, * PLUS ONE. * COEF(NS,I) CONTAINS STOICHIOMETRIC COEFFICIENT FOR SPECIES "SPEC * (NS,I)" IN THE DERIVED SPECIES. * SPEC(NS,I) IDENTIFIES COMPONENT SPECIES IN THE DERIVED SPECIES. * * THE I = 1 POSITION IS LEFT BLANK. * IF H2 IS IN THE DERIVED GAS, IT SHOULD BE IN THE I=2 POSITION. * TEST FOR END OF DERIVED SPECIES LIST BY BLANK LINE. IF(NAME(NS).NE.BLANK) GO TO 115 NST = NS - 1 GO TO 120 115 CONTINUE IT(NS) = ITOT IF(ITOT.LE.MAXSTO) GO TO 112 WRITE(IOUT1,662) NAME(NS), ITOT, MAXSTO 662 FORMAT(/5X,'SPECIES ',A8,' IN GASTHERM HAS TOO MANY COMPONENTS' $ /5X,'ITOT FOR THAT SPECIES IS ',I3,' ; MAXIMUM IS ',I3//) STOP * THE I = 1 POSITION IS LEFT BLANK * IF H2 IS IN THE COMPLEX, IT SHOULD BE IN THE I=2 POSITION 112 READ(INP2,818) (COE(I),I = 1, 5) * COE CONTAINS POWER FUNCTION COEFFICIENTS TO GENERATE LOG K AS F(T) DO 121 I = 1,5 121 APFCL(NS,I) = COE(I) CALL POWERK(COE,TEMPK,AKLOG) * * GASTHERM CONTAINS ONE SET OF COEFFICIENTS TO CALCULATE LOG K AS * F (T) WHERE T IS IN DEGREES KELVIN. THE ABOVE STATEMENTS * STORE THE POWER FUNCTION COEFFICIENTS THAT WILL BE NEEDED LATER * IN ARRAYS APFCL. * 119 LGK(NS) = AKLOG * * ALL GAS SPECIES LISTED WITH A NON BLANK CHARACTER IN COL 16 (SUP) * WILL BE SUPPRESSED FROM THE CALCULATION, DISPLAYED ON THE SCREEN * AND PRINTED ON THE OUTPUT. * IF(SUP.EQ.BLNK) GO TO 117 WRITE(IOUT2,885) NAME(NS) WRITE(IOUT1,885) NAME(NS) GO TO 114 117 CONTINUE * * WE ALSO SUPPRESS THE GASES MENTIONED IN THE INPUT FILE DO 130 I = 1,NSUP IF(SUPNAM(I).NE.NAME(NS)) GO TO 130 WRITE(IOUT2,886) NAME(NS) WRITE(IOUT1,886) NAME(NS) GO TO 114 130 CONTINUE * * SORT OUT GASES CONTAINING COMPONENTS BEING USED * AND REDEFINE GAS COMPONENT INDICIES, * ALSO SUMS UP ALL THE GAS SPECIES COEFFICIENTS (COSUM) * COSUM = 0. DO 118 I = 2, ITOT JN = SPEC(NS,I) IF(QAS(JN).EQ.0) GO TO 114 SPEC(NS,I) = QAS(JN) COSUM = COSUM + COEF(NS,I) 118 CONTINUE COEF(NS,1) = -1. SPEC(NS,1) = NS GAMMA(NS) = 1.0 AZERO(NS) = COSUM - 1. GO TO 113 * 120 CONTINUE read (inp2,*) DO 126 I = 1, SAQTOT AQM(I)=MTOT(I) 126 TMIN(I) = 0. * ** CALL HELHUC * *********** NM = 0 DO 26 I = 1,MAXSSL 26 SUMSOL(I) = 0.0 ** ** HERE WE READ COEFFICIENTS TO COMPUTE CROSS VIRIAL COEFFICIENTS FOR ** H2O-CO2-CH4 GAS MIXTURE. B12 IS FOR B(H2O-CO2), B13 FOR B(H2O-CH4 ** B23 FOR B(CO2-CH4), C112 FOR C(H2O-H2O-CO2) ETC. SEE SUBROUTINE ** PHI. THESE COEFFICIENTS ARE ASSIGNED AS FOLLOW: ** ** CVB(1) : B12 ** CVB(2) : B13 ** CVB(3) : B23 ** CVC(1) : C112 ** CVC(2) : C122 ** CVC(3) : C113 ** CVC(4) : C133 ** CVC(5) : C223 ** CVC(6) : C233 ** CVC(7) : C123 ** ** AS OF FEB. 86 WE DO NOT HAVE DATA FOR C113, C133 AND C123. THESE ** VALUES ARE ASSUMED TO BE ZERO. THIS SEEMS TO BE A BETTER ** ASSUMPTION THAN TO CALCULATE THEM WITH A MIXING RULE (LIKE ** A GEOMETRIC AVERAGE). ** ** DO 28 I=1,3 ** READ(INP2,817) (CVBCF(I,J), J=1,3) ** *************** ** 28 CVB(I) = CVBCF(I,1)/TEMPK/TEMPK + CVBCF(I,2)/TEMPK + CVBCF(I,3) ** DO 29 I=1,7 ** READ(INP2,817) (CVCCF(I,J), J=1,3) ** *************** ** 29 CVC(I) = CVCCF(I,1)/TEMPK/TEMPK + CVCCF(I,2)/TEMPK + CVCCF(I,3) ** ** 33 CONTINUE NM = NM + 1 IF(NM.LE.MAXNM) GO TO 30 NM = NM-1 WRITE(IOUT1,664) MAXNM, MIN(NM) 664 FORMAT(/5X,'NUMBER OF MINERALS READ IN SOLTHERM EXCEEDS MAXIMUM' $ ,' OF ',I3,/5X,'LAST NAME READ IS: ',A8//) STOP * read (inp2,*) 30 READ(INP2,854) MIN(NM),SUP,XMOLW(NM),ITOT,(COEFM(NM,I),SPECM(NM,I) 1 ,I = 1, ITOT), B(NM), JSOLS(NM) * *************** * READ MINERAL CHARACTERISTICS FROM GASTHERM. * IF SUP NOT BLANK OR NOT 'G' = SUPPRESS MINERAL. SUP = 'G' FOR * GASES (WILL CAUSE CHILLER TO COMPUTE AND OUTPUT FUGACITY). * XMOLW = MOLECULAR WEIGHT. * COEFM & SPECM ARE AS FOR DERIVED SPECIES, ABOVE. * B(NM) CONTAINS EXPONENT FOR MINERAL END MEMBER TO BE USED TO * CALCULATE ACTIVITY USING IDEAL MULTI-SITE MIXING. * JSOLS(NM) IDENTIFIES THE PARTICULAR SOLID SOLUTION, IF ANY, THAT * THIS MINERAL BELONGS TO. ITM(NM) = ITOT * BLANK LINE IN GASTHERM ENDS READING OF MINERALS. IF(MIN(NM).EQ.BLANK) GO TO 165 * IF(ITOT.LE.MAXSTO) GO TO 34 WRITE(IOUT1,666) MIN(NM), ITOT, MAXSTO 666 FORMAT(/5X,'MINERAL ',A8,' IN GASTHERM HAS TOO MANY COMPONENTS' $ /5X,'ITOT FOR THAT MINERAL IS ',I3,' ; MAXIMUM IS ',I3//) STOP * * SPECIES LISTED WITH CHARACTER OTHER THAN BLANK IN COL 9 * ARE SUPPRESSED FROM THE CALCULATION. THESE MINERALS ARE DISPLAYED * ON THE SCREEN AND PRINTED ON OUTPUT. * 34 IF(SUP.EQ.BLNK) GO TO 35 WRITE(IOUT2,885) MIN(NM) WRITE(IOUT1,885) MIN(NM) READ(INP2,815) GO TO 30 35 CONTINUE * WE ALSO SUPPRESS THE SPECIES MENTIONED IN THE INPUT FILE DO 135 I = 1,NSUP IF(SUPNAM(I).NE.MIN(NM)) GO TO 135 WRITE(IOUT2,886) MIN(NM) WRITE(IOUT1,886) MIN(NM) READ(INP2,815) GO TO 30 135 CONTINUE * * DETERMINES IF MINERAL WAS GIVEN IN SATURATED MINERAL LIST * AND SET SEQ NUMBERS * SEQ CONTAINS THE SEQUENCE POSITIONS OF MINERALS CURRENTLY * INCLUDED IN EQUILIBRIUM ASSEMBLAGE. THESE INTEGERS ARE NOT * SIMPLY SEQUENCE POSITIONS IN SOLTHERM, BUT SEQUENCE POSITIONS * IN THE RELEVANT MINERAL LIST (AS SORTED BY GASWORKS FROM * GASTHERM ACCORDING TO THE COMPOSITION OF THIS PARTICULAR RUN). IF(NMT.EQ.0) GO TO 38 DO 37 I = 1, NMT IF(MIN(NM).NE.SATMIN(I)) GO TO 37 SEQ(I) = NM GO TO 38 37 CONTINUE * 38 SCALE(NM)=0. DO 40 I = 1, ITOT JS = SPECM(NM,I) * SKIP IF STOICHIOMETRY INCLUDES SPECIES THAT WERE NOT GIVEN IF(QAS(JS).NE.0) GO TO 42 READ(INP2,815) GO TO 30 42 SPECM(NM,I) = QAS(JS) SCALE(NM)=SCALE(NM)+DABS(COEFM(NM,I)) * * SCALE(NM) IS USED IN MINSAT TO SCALE THE LOG Q/K VALUES FOR MINERAL * ACCORDING TO THE NUMBER OF COMPONENT SPECIES IN THE MINERAL. * 40 CONTINUE * READ(INP2,818) (COE(I),I = 1,5) * ********************* * COE CONTAINS POWER FUNCTION COEFFICIENTS FOR GENERATING LOG K=F(T) * DO 41 I = 1,5 41 BPFCL(NM,I)=COE(I) CALL POWERK(COE,TEMPK,AKLOG) * ********** * BPFCL STORES POWER FUNCTION COEFFICIENTS FOR LATER RE-USE * THE ABOVE SEQUENCE OF GO TO STATEMENTS DO FOR MINERALS WHAT LINES * 121 AND 122 DO FOR GASES ABOVE. * 45 CONTINUE IF(LIMSOL.EQ.0) JSOLS(NM) = 0 * TURNS OFF SOLID SOLUTIONS IF ALL ARE TURNED OFF BY LIMSOL = 0. * LGKM(NM) = AKLOG J = 0 GO TO 33 165 CONTINUE IF(IDEAL.GT.2) IDEAL = 2 * NMTT = NM - 1 IF(NMT.EQ.0) GO TO 175 DO 48 I = 1, NMT * * DO 48 SORTS OUT MINERALS CURRENTLY INCLUDED, AS SPECIFIED IN * WORKRUN WHERE MINERAL NAMES ARE INPUT * NM = SEQ(I) IF(NM.NE.0) GO TO 49 WRITE(IOUT1,656) SATMIN(I) WRITE(IOUT2,656) SATMIN(I) 656 FORMAT(/3X,A8,', SPECIFIED IN WORKRUN,' , $ ' WAS NOT FOUND IN GASTHERM. CHECK SPELLING.'//) STOP 49 QES(NM)=I MINTRY(NM)=ACT(I) * * SETS TRIAL MOLES OF MINERALS AS INPUT IN ACT ARRAY IN CHILLRUN. * 48 CONTINUE * 175 CONTINUE * BSTOT = SAQTOT BSTP1 = SAQTOT + 1 JSTOP = 0 * * * ****************************************************************** * END INPUT BEGIN DATA ECHO * ****************************************************************** * ****************************************************************** * * WRITE(IOUT1,970) WRITE(IOUT1,943) WRITE(IOUT1,960) TITLE1 WRITE(IOUT1,960) TITLE2 WRITE(IOUT1,943) WRITE(IOUT1,970) WRITE(IOUT1,923) WRITE(IOUT1,940) WRITE(IOUT1,945) (NAME(I),I,MTOT(I), MTRY(I), GAMMA(I), 1 COMTOT(I), I=1,BSTOT) * 50 WRITE(IOUT1,611) ERPC, FO2LOG, PFLUID, SLIM, TEMP,SINC,TEMPC, 1 TOTMIX, SOLMIN, RM, GASGRM,SUPRNT,GASTOT WRITE(IOUT2,611) ERPC, FO2LOG, PFLUID, SLIM, TEMP,SINC,TEMPC, 1 TOTMIX, SOLMIN, RM,GASGRM,SUPRNT,GASTOT 611 FORMAT(//4X,'THE FOLLOWING PARAMETERS WERE GIVEN:' 1 /6X,'ERPC = ',D10.4,' FO2LOG = ',F10.5,' PFLUID = ',F10.5, 2 /6X,'SLIM = ',G10.4,' TEMP = ',F10.5,' SINC = ',G10.4, 3 /6X,'TEMPC = ',F10.5,' TOTMIX = ',G10.4,' SOLMIN = ',D10.4, 4 /6X,'RM = ',F10.5,' GASGRM = ',F10.5,' SUPRNT = ',D10.4, 5 /6X,'GASTOT= ',F10.5) WRITE(IOUT1,612) IFRAC, IPUNCH, NLOOP, ISTEP 1, LIMSS, LOOC, ITREF, IDEAL,INCREM,INCP,MINSOL,OFLAG WRITE(IOUT2,612) IFRAC, IPUNCH, NLOOP, ISTEP 1, LIMSS, LOOC, ITREF, IDEAL, INCREM,INCP,MINSOL,OFLAG 612 FORMAT(//4X,'THE FOLLOWING OPTIONS WERE SELECTED:' 1 /6X,'IFRAC =',I5,' IPUNCH =',I5,' NLOOP =',I5,' ISTEP =',I5, 2 /6X,'LIMSOL=',I5,' LOOC =',I5,' ITREF =',I5,' IDEAL =',I5, 3 /6X,'INCREM=',I5,' INCP =',I5,' MINSOL =',I5,' OFLAG =',I5) IF(NPP.NE.0) WRITE(IOUT1,614) NPP 614 FORMAT(6X,'GIVEN PRESSURE GRID. ',I3,' T-P POINTS READ'/) WRITE(IOUT1,923) WRITE(IOUT1,617) SUPRNT 617 FORMAT(2X,'OF THE DERIVED SPECIES LISTED BELOW, ONLY THOSE WITH' 1 ,/2X,'MOLE NUMBERS ABOVE ', D10.4,' WILL APPEAR ON THE OUTPUT'/) WRITE(IOUT1,950) * DOUBLE COLUMN ECHO OF DERIVED SPECIES AND MINERALS FROM CHILLER * IMPORTED HERE NSTMBS = NST - BSTOT K = NSTMBS/2 IF(MOD(NSTMBS,2).NE.0) K = K + 1 WRITE(IOUT1,955) (NAME(I), I, LGK(I), NAME(I+K), 1 I+K, LGK(I+K), I = BSTP1,BSTOT+K) IF(NMTT.EQ.0) GO TO 176 WRITE(IOUT1,900) L = NMTT/2 IF(MOD(NMTT,2).NE.0) L = (NMTT+1)/2 WRITE(IOUT1,910) (MIN(I), I, LGKM(I), MINTRY(I), MIN(I+L), 1 I+L, LGKM(I+L), MINTRY(I+L), I = 1,L) * IF(MINSOL.EQ.0) GO TO 176 IF(MINSOL.LE.2) WRITE(IOUT1,974) IF(MINSOL.GT.2) WRITE(IOUT1,976) WRITE(IOUT1,975) (NOMOX(I), WTPC(I), I = 1, NOOX) GRMLG = -100. IF(TOTMIX .GT.0.) GRMLG = DLOG10(TOTMIX) WRITE(IOUT1,971) TOTMIX, TYP, GRMLG WR = 0. IF(TOTMIX.NE.0.) WR = GASGRM/TOTMIX WRITE(IOUT1,973) GASGRM, WR * ****************************************************************** * 176 CONTINUE IF(JSTOP.EQ.1) GO TO 388 * CALL SQUOMP * *********** * IF(NLOOP.NE.0) GO TO 170 WRITE(IOUT2, 616) WRITE(IOUT1, 616) 616 FORMAT(/6X,'END DATA ECHO. SET NLOOP > 0 TO START CALCULATIONS'/) STOP 170 CONTINUE * ** TO EASE CONVERGENCE W/ PICKUP RUN, WE CALCULATE GAMMAS SOLID SOL. ** NM = SOLS(10,1) ** JS = ISOLS(NM) * ******** 189 IF(INCREM.NE.0) GO TO 560 * * STEP STARTS HERE * =============================================================== * 190 CONTINUE * N IS THE TOTAL NUMBER OF UNKNOWNS 191 CONTINUE ERLOOP = 0 NDUMP= 0 MINCT= 0 JSTOP= 0 * GASTOT EQUATION IS IN THE BSTOT + 1 POSITION IN THE MATRIX JGSTOT = BSTOT + 1 N = BSTOT + NMT + 1 IF(N.LE.MAXN) GO TO 192 WRITE(IOUT1,668) N, MAXN 668 FORMAT(/5X,'SIZE OF MATRIX REACHED LIMIT. N = ',I3,' MAXIMUM' $ ,' SIZE IS ',I3//) STOP 192 IF(NMT.LE.MAXIN) GO TO 193 WRITE(IOUT1,670) NMT, MAXIN 670 FORMAT(/5X,'TOO MANY PHASES IN MATRIX (',I3,'); MAXIMUM' $ ,' NUMBER ALLOWED IS ',I3//) STOP * * THE PRECEDING "IF" SETS AN APPROXIMATE FH2, ASSUMING THAT H20 IS * THE OVERWHELMINGLY DOMINANT GAS * HERE, WE CALCULATE CONVERGENCE LIMITS FOR F'S BELOW * ERPC SHOULD BE D-12 OR LARGER (32 BIT MACHINE) 193 DO 194 I = 1,N IF (I.EQ.JGSTOT) THEN CONV(I)=ERPC*GSTOLD ELSE IF (I.GT.JGSTOT) THEN II = I-JGSTOT NM = SEQ(II) CONV(I) = ERPC*DABS(LGKM(NM)) IF(DABS(LGKM(NM)).EQ.0.) CONV(I) = ERPC*10. ELSE IF (I.EQ.1.AND.FO2.NE.0.) THEN CONV(1) = ERPC*DABS(LGK(BSTP1)) ELSE IF (DABS(MTRY(I)).GT.DABS(MTOT(I))) THEN CONV(I) = ERPC*DABS(MTRY(I)) ELSE IF (MTOT(I).EQ.0.) THEN CONV(I) = ERPC ELSE CONV(I) = ERPC*DABS(MTOT(I)) ENDIF 194 CONTINUE * * ****************************************************************** * * HERE WE GO TO 313 IN FIRST LOOP TO COMPUTE MTRYS OF DERIVED * SPECIES BEFORE SETTING THE MATRIX GO TO 313 * * >>>>> ITERATIONS START HERE <<<<< * ================================================== * 200 CONTINUE * 202 ERLOOP = ERLOOP + 1 NSTOP = 0 IF(FO2LOG.EQ.0.) GO TO 210 * ABOVE IF SKIPS SUMMATION OF H2-BEARING SPECIES TO GET MTOT(H2) IF * FO2 IS INPUT IF(ERLOOP.EQ.1) GO TO 225 H2TOT = MTRY(1) * * SUMS H2 SPECIES FOR MTOT OF H2 IF FO2 IS INPUT DO 204 NS = BSTP1, NST IF(SPEC(NS,2).NE.1) GO TO 204 H2TOT = H2TOT + MTRY(NS)*COEF(NS,2) 204 CONTINUE H2TOT = H2TOT + TMIN(1) * MTOT(1) = H2TOT 210 CONTINUE 225 CONTINUE * * ****************************************************************** * DO 250 I = 1,N DO 245 J = 1,N A(I,J) = 0.0 245 CONTINUE F(I) = 0. A(I,I) = 1. 250 CONTINUE * INITIALIZE THE MATRIX VALUES FOR THE GASTOT EQUATION A(JGSTOT,JGSTOT)=-1. F(JGSTOT)=-GASTOT DO 300 KK = 1, BSTOT F(JGSTOT) = F(JGSTOT) + MTRY(KK) A(JGSTOT,KK) = 1. * * * THE FOLLOWING IF'S CONTROL PRODUCTION OF MATRIX ENTRIES * FOR O2 IF(KK.EQ.1.AND.FO2LOG.NE.0.) GO TO 300 * * * F(KK) = -MTOT(KK) + MTRY(KK) DO 290 NS = BSTP1, NST * THE DO 290 LOOP THROUGH ALL DERIVED SPECIES FOR INCLUSION IN * APPROPRIATE MASS BALANCE EQUATIONS. ITOT = IT(NS) DO 265 IC = 2,ITOT IF(SPEC(NS,IC).EQ.KK) GO TO 267 265 CONTINUE GO TO 290 267 CONTINUE CFF = COEF(NS,IC) A(KK,JGSTOT) = A(KK,JGSTOT) - CFF*MTRY(NS)*AZERO(NS)/GASTOT A(JGSTOT,KK) = A(JGSTOT,KK) + CFF*MTRY(NS)/MTRY(KK) DO 288 II = 2, ITOT J = SPEC(NS,II) * DTERM IS THE DERIVATIVE OF MOLALITY EXPRESSION FOR DERIVED SPECIES * WITH RESPECT TO COMPONENT SPECIES J (= SPEC(NS,II) DTERM = COEF(NS,II) * MTRY(NS)/MTRY(J) A(KK,J) = A(KK,J) + CFF*DTERM 288 CONTINUE F(KK) = F(KK) + CFF*MTRY(NS) 290 CONTINUE 300 CONTINUE DO 301 NS = BSTP1,NST F(JGSTOT) = F(JGSTOT) + MTRY(NS) A(JGSTOT,JGSTOT) = A(JGSTOT,JGSTOT) - AZERO(NS)* 1 MTRY(NS)/GASTOT 301 CONTINUE * MINERA WILL SET UP MASS ACTION EQS. FOR MINERALS IF(NMT.NE.0) CALL MINERA * *********** IF(FO2LOG.NE.0.) CALL OXY * ******** * PROCEDURE TO HELP CONVERGENCE (SUPPOSEDLY !!). SKIP IF RM=0 OR 1 * TRY ONLY IF THERE ARE CONVERGENCE PROBLEMS (LAST RESOURCE). * ASSURES THAT ALL F'S DECREASE AFTER EACH ITERATION BY MULTIPLYING * F'S BY POWERS OF RM (GENERALLY RM IS BETWEEN .1 AND .9). IF(RM.EQ.0..OR.RM.GE.1.) GO TO 297 DO 292 I = 1,N FFF = F(I) IF(DABS(F(I)).LT.DABS(FOLD(I))) GO TO 296 FACT(I) = FACT(I) * RM F(I) = FFF * FACT(I) IF(FACT(I).LT.1.D-4) FACT(I) = 1.D-4 GO TO 294 296 FACT(I) = 1.0 294 FOLD(I) = FFF 292 CONTINUE * * 297 IF(ERLOOP.LE.5) GO TO 365 * CONVERGENCE CRITERION. SEES THAT ALL F(I) < CONV(I) DEFINED * AFTER STATEMENT 191 DIFTOT = GASTOT DO 360 NS = 1, NST DIFTOT = DIFTOT - MTRY(NS) 360 CONTINUE IF(DABS(DIFTOT).GT.ERPC) GO TO 365 DO 363 I = 1,N IF(DABS(F(I)).GT.CONV(I)) GO TO 365 363 CONTINUE GO TO 410 * 365 IF(ISTEP.EQ.0) GO TO 308 IF(ERLOOP.NE.2) GO TO 308 306 NSTOP = 1 GO TO 388 307 NSTOP = 0 NDUMP = 1 CALL DUMP * ********* 308 CONTINUE IF(ERLOOP.EQ.1.AND.ISTEP.EQ.1) CALL DUMP CALL SLUD * ********* CALL SIR(X) * ******** IF(ERLOOP.EQ.1.AND.ISTEP.EQ.1) THEN WRITE(IOUT1,977) (X(I), I = 1, N) ENDIF IF(NDUMP.EQ.1.AND.ISTEP.EQ.1) CALL DUMP2 * ********** NDUMP = 0 FMAX = PFLUID * GSTOLD PLOGP1 = DLOG10(FMAX) DO 310 KK = 1, BSTOT COLD = MTRY(KK) MTRY(KK) = MTRY(KK) + X(KK) * LIMIT SIZE OF MTRY VALUES TO AID CONVERGENCE IF(MTRY(KK).GT.FMAX) MTRY(KK) = COLD XX = X(KK) 309 CONTINUE IF(MTRY(KK).GT.0.) GO TO 320 XX = XX/2. MTRY(KK) = COLD + XX GO TO 309 320 CONTINUE IF(MTRY(KK).GT.RMIN) GO TO 310 WRITE(IOUT2,676) NAME(KK), RMIN 676 FORMAT(5X,'** WARNING ! MOLALITY OF ',A8,' LESS THAN ', D10.3) MTRY(KK) = RMIN 310 CONTINUE * BRACKETS GASTOT TO EASE CONVERGENCE GASTOT = GASTOT + X(JGSTOT) GASMAX = GSTOLD*10. GASMIN = GSTOLD/10. IF (GASTOT.GT.GASMAX) GASTOT = GASMAX IF (GASTOT.LT.GASMIN) GASTOT = GASMIN NNMT = NMT IF(NMT.EQ.0) GO TO 313 IF(PFLUID.GT.PMAX) PFLUID = PMAX * PLOG = DLOG10(PFLUID) DO 414 I = 1,MAXSSL 414 SOLTRY(I) = 0.0 IQ = 0 DO 314 I = 1, NMT IQ = IQ + 1 NM = SEQ(IQ) M = BSTOT + I + 1 MINTRY(NM) = MINTRY(NM) + X(M) * IF(ISOLS(NM).EQ.0) GO TO 417 JU = ISOLS(NM) IF(MINTRY(NM).LE.0.) GO TO 311 SOLTRY(JU) = SOLTRY(JU) + MINTRY(NM) GO TO 314 * WE KICK OUT SOLSOL IF LE 0 AND LIMSS = 2 311 IF(ERLOOP.GT.6.AND.LIMSS.GT.1) GO TO 420 * WE RESET NEGATIVE SOL SOL TO POSITIVE IF NO KICK OUT MINTRY(NM) = 1.D-20 SOLTRY(JU) = SOLTRY(JU) + MINTRY(NM) GO TO 314 * * PURE MINERAL: WE KICK OUT MINERAL IF MINTRY LT SOLMIN 417 IF(ERLOOP.LE.4.OR.SOLMIN.GE.0.) GO TO 314 IF(MINTRY(NM).GE.SOLMIN) GO TO 314 * * MINERAL KICK OUT SCHEME BELOW 420 WRITE(IOUT2,777) MIN(NM), MINTRY(NM) WRITE(IOUT1,777) MIN(NM), MINTRY(NM) 777 FORMAT(5X,A8,' MINTRY= ',D10.3,' MINERAL REMOVED FROM MATRIX') SEQ(IQ) = 0 QES(NM) = 0 MINTRY(NM) = 1.D-300 * SORT QES AND SEQ AND CONV WITHOUT CURRENT NM MINERAL NNMT = NNMT - 1 * IF NO MINERALS ARE LEFT IN MATRIX, WE SKIP RESORTING IF(NNMT.EQ.0) GO TO 314 DO 418 L = IQ,NNMT LL = L+1 J = BSTOT + L + 1 JJ = J + 1 SEQ(L) = SEQ(LL) CONV(J) = CONV(JJ) F(J) = F(JJ) 418 QES(SEQ(L)) = L IQ = IQ - 1 IF(ISOLS(NM).EQ.0) GO TO 314 * * IF WE KICKED OUT A SOLID SOLUTION ENDMEMBER ISOLS(NM) = 0 JS = JSOLS(NM) IS = IST(JS) ISTOT = 0 * WE COUNT HOW MANY ENDMEMBERS LEFT IN SOLID SOLUTION (ISTOT) DO 336 IL = 1,IS LS = SOLS(JS,IL) IF(QES(LS).EQ.0) GO TO 336 ISTOT = ISTOT + 1 LSBR = LS 336 CONTINUE * * IF ONLY ONE END MEMBER LEFT, WE BREAK SOLID SOLUTION * BY SETTING ISOLS TO 0 IF(ISTOT.NE.1) GO TO 314 WRITE(IOUT1,778) MIN(LSBR) WRITE(IOUT2,778) MIN(LSBR) 778 FORMAT(5X,A8,' NOT IN SOLID SOLUTION ANYMORE') ISOLS(LSBR) = 0 * 314 CONTINUE * * * NMT AND N HAVE TO BE RESET IN CASE WE KICKED OUT MINERALS NMT = NNMT N = NNMT + BSTOT + 1 313 CONTINUE * ** HERE WE COMPUTE ACTIVITY COEFFICIENTS FOR ** NON IDEAL SOLID SOLUTIONS THAT ARE SATURATED (CURRENTLY IN ** MATRIX). WE HAVE ONLY GOLD SO FAR SOL.SOL ID NO. 10 ** ** NM = SOLS(10,1) ** JS = ISOLS(NM) ** IF(JS.NE.0) CALL ELECTR ** *********** BSTP2 = BSTP1 JK = 1 IF(FO2.EQ.0.) GO TO 315 MTRY(BSTP1) = FO2 * GASTOT/(GAMMA(BSTP1) * PFLUID) ACT(BSTP1) = GAMMA(BSTP1) * MTRY(BSTP1) * PFLUID/GASTOT ACT(1) = MTRY(1)*GAMMA(1)*PFLUID/GASTOT JK = 2 BSTP2 = BSTP1 + 1 315 DO 319 NS=BSTP2,NST CSAVE= -DLOG10(GAMMA(NS)) - LGK(NS) ITOT = IT(NS) DO 321 I = 2, ITOT J = SPEC(NS,I) CSAVE2 = (DLOG10(GAMMA(J) * MTRY(J))) * COEF(NS,I) CSAVE= CSAVE + CSAVE2 321 CONTINUE PTERM=AZERO(NS) * DLOG10(PFLUID/GASTOT) RSL=CSAVE + PTERM ATERM = RSL + DLOG10(GAMMA(NS) * PFLUID/GASTOT) IF(RSL.LT.PLOGP1) GO TO 316 RSL = PLOGP1 ATERM = PLOGP1 316 IF(RSL.GT.RLOGMN) GO TO 317 RSL = RLOGMN ATERM = RLOGMN 317 MTRY(NS)=10.**RSL ACT(NS) = 10.**ATERM 319 CONTINUE IF(ERLOOP.EQ.0) GO TO 200 * * ****************************************************************** WRITE(IOUT2,650) ERLOOP, GASTOT 650 FORMAT(5X,'LOOP NO.',I4,' GASTOT: ',D15.7) IF(ERLOOP.EQ.NLOOP) GO TO 410 * COMPUTE THE FUGACITIES (ACT) OF THE COMPONENT SPECIES 352 DO 358 NS = JK, BSTOT FTERM = DLOG10(MTRY(NS)) ATERM = FTERM + DLOG10(GAMMA(NS) * PFLUID/GASTOT) ACT(NS) = 10.**ATERM 358 CONTINUE IF(ISTEP.EQ.1) WRITE(IOUT1,905) ERLOOP,GASTOT IF(ISTEP.EQ.0) GO TO 366 DO 362 NS = 1, BSTOT 362 WRITE(IOUT1,920) NS, NAME(NS), MTRY(NS) WRITE(IOUT1,977) (X(I), I = 1, N) 366 CONTINUE * WE COMPUTE AQM, THE QUANTITY OF GASEOUS COMPONENTS DO 350 I=1,SAQTOT AQM(I) = MTRY(I) DO 370 NS = BSTP1, NST ITOT = IT(NS) DO 364 IC = 2,ITOT IF(SPEC(NS,IC).EQ.I) GO TO 367 364 CONTINUE GO TO 370 367 CONTINUE CFF = COEF(NS,IC) AQM(I) = AQM(I) + CFF*MTRY(NS) 370 CONTINUE 350 CONTINUE GO TO 200 * * 410 CONTINUE GSTOLD = GASTOT * SUM H2-BEARING SPECIES IF FO2 IS INPUT IF (FO2LOG.EQ.0.) GO TO 412 H2TOT = MTRY(1) ACT(1) = MTRY(1)*GAMMA(1)*PFLUID/GASTOT DO 411 NS=BSTP1,NST IF(SPEC(NS,2).NE.1) GO TO 411 H2TOT = H2TOT + COEF(NS,2)*MTRY(NS) 411 CONTINUE H2TOT = H2TOT + TMIN(1) MTOT(1) = H2TOT 412 JSTOP = 1 388 CONTINUE WRITE(IOUT1,999) WRITE(IOUT1,960) TITLE1 WRITE(IOUT1,960) TITLE2 WRITE(IOUT1,943) WRITE(IOUT1,970) WRITE(IOUT1,946) ERLOOP,NLOOP WRITE(IOUT1,943) WRITE(IOUT1,933) TEMP, PFLUID, TOTMIX WRITE(IOUT1,927) GASTOT WRITE(IOUT1,944) WRITE(IOUT1,914) WRITE(IOUT1,915) GRT=0. DO 400 NS = 1, NST IF(NS.LE.BSTOT) GRT=GRT+ AQM(NS)*W(NS) IF(NS.GT.BSTOT) GO TO 396 ALA = DLOG10(ACT(NS)) ** ALGAM = DLOG10(GAMMA(NS)) * ALM IS THE PARTIAL PRESSURE ALM = ACT(NS)/GAMMA(NS) * ALMF IS THE MOLE FRACTION ALMF = ALM/PFLUID WRITE(IOUT1,920) NS,NAME(NS),MTRY(NS),ALMF,ALM,ACT(NS),ALA, 1 GAMMA(NS) GO TO 400 396 IF(MTRY(NS).LT.SUPRNT.AND.NS.NE.BSTP1) GO TO 400 ALA = DLOG10(ACT(NS)) ** ALGAM = DLOG10(GAMMA(NS)) ALM = ACT(NS)/GAMMA(NS) ALMF = ALM/PFLUID WRITE(IOUT1,920) NS,NAME(NS),MTRY(NS),ALMF,ALM,ACT(NS),ALA, 1 GAMMA(NS) 400 CONTINUE GRAMSO = GRT * 402 WRITE(IOUT1,923) GRT = 0. IF(NMT.EQ.0) GO TO 398 WRITE(IOUT1,957) DO 404 NS = 1, BSTOT TMOLAL = AQM(NS)/GASTOT PPM=1.D06 * AQM(NS) * W(NS)/GRAMSO * GRAMSO CONTAINS GRT, WHICH IS TOTAL MASS OF GAS PHASE, * THUS PPM IS PPM IN GAS PHASE OF GIVEN COMPONENT GAS SPECIES * TMOLAL CONTAINS A BULK MOLE FRACTION (OF SORTS) WRITE(IOUT1,962) NS, NAME(NS), MTRY(NS), TMOLAL, PPM 404 CONTINUE WRITE(IOUT1,965) GRAMSO IF(GASGRM.EQ.0.) GASGRM = GRAMSO WRITE(IOUT1,961) DO 405 NS = 1, BSTOT * TMIN WAS COMPUTED IN MINERAL. IT IS THE TOTAL NUMBER * OF MOLES OF COMPONENTS PRESENT IN SOLIDS WRITE(IOUT1,962) NS, NAME(NS), MTOT(NS), AQM(NS), TMIN(NS) 405 CONTINUE WRITE(IOUT1,923) WRITE(IOUT1,930) WRITE(IOUT2,923) WRITE(IOUT2,930) GRT = 0. DO 427 I = 1, NMT NM = SEQ(I) GRT = GRT + MINTRY(NM) * XMOLW(NM) 427 CONTINUE VOLMIN = 0. DO 395 I = 1, NMT NM = SEQ(I) GRAMS = MINTRY(NM) * XMOLW(NM) ALGG=1. ALGM=1. IF(MINTRY(NM).LT.0.) GO TO 429 ALGM = DLOG10(MINTRY(NM)) ALGG = DLOG10(GRAMS) 429 CONTINUE XMIN = 1.0 JS = ISOLS(NM) IF(JS.NE.0) XMIN = MINTRY(NM)/SOLTRY(JS) ACTMIN = GAMIN(NM)*XMIN IF(JS.NE.0) ACTMIN = ACTMIN**B(NM) LOGACT = DLOG10(ACTMIN) IF(GRT.NE.0.) PCMIN = GRAMS/GRT * 100. ** IF(NM.GT.ITGAS) GO TO 435 ** PI = XMIN*PFLUID ** FI = DLOG10(PI*GAMIN(NM)) ** CALCULATES VOLUME FROM PV/RT = Z, Z CALCULATED FROM BVAL AND ** AND CVAL IN SUBROUTINE PHI ** ZGAS = 1 + BVAL(NM)*PFLUID + CVAL(NM)*PFLUID*PFLUID ** CM3 = TEMPK * 83.1432 * ZGAS * MINTRY(NM) / PFLUID * WRITES SATURATED MINERALS DATA 435 IF(DENS(NM).EQ.0.) DENS(NM) = 2.5 CM3 = GRAMS/DENS(NM) VOLMIN = VOLMIN + CM3 WRITE(IOUT1,935) MIN(NM), MINTRY(NM), ALGM, GRAMS, 1 ALGG,PCMIN,CM3 WRITE(BUFOUT(I),937) MIN(NM),XMIN,ACTMIN,GAMIN(NM),LOGACT * TERMINAL OVERLOOK 395 WRITE(IOUT2,935) MIN(NM), MINTRY(NM), ALGM, GRAMS, 1 ALGG,PCMIN,CM3 WRITE(IOUT1,931) WRITE(IOUT1,'(A80)') (BUFOUT(I),I=1,NMT) 408 WRITE(IOUT1,943) * 398 IF(TEMPC.NE.0.) WRITE(IOUT1,936) TOTMIX WRITE(IOUT1,972) GRT, VOLMIN WRITE(IOUT1,968) IF(TEMPC.EQ.0.) GO TO 382 * COMPUTE WEIGHT OF MIXING SOLUTION TOTWT = 0. DO 383 I = 1, BSTOT 383 TOTWT = TOTWT + TOTMIX * COMTOT(I) * W(I) IF(TOTMIX.NE.0.) WR = GASGRM/TOTWT GO TO 380 * 382 IF(MINSOL.EQ.0) GO TO 380 GRMLG = -100. IF(TOTMIX.GT.0.) GRMLG=DLOG10(TOTMIX) WRITE(IOUT1,971) TOTMIX, TYP, GRMLG IF(TOTMIX.NE.0.) WR = GASGRM/TOTMIX WRITE(IOUT1,973) GASGRM, WR 380 IF(NSTOP.EQ.1) GO TO 307 ** ** NON IDEAL SOLID SOLUTIONS. WE NEED TO COMPUTE ACTIVITIES BEFORE ** CALL TO MINSAT, EVEN IF MINSAT WILL RECOMPUTE THESE ACTIVITIES. ** THIS STEP IS ONLY FOR PHASES THAT ARE NOT YET SATURATED. ** WE NEED ALL END MEMBERS ACTIVITIES TO CHECK IF THE SOLID ** SOLUTION SHOULD BE SUPERSATURATED. ** SO FAR, WE HAVE ONLY ELECTRUM SOL-SOL NO. = 10 ** SILVER SHOULD BE ENTERED AFTER GOLD IN SOLTHERM IS=1 AU, =2 AG ** NM = SOLS(10,1) ** LS = IST(10) ** GOLD NOT IN THE CALCULATION, OR ALREADY SATURATED => SKIP ** IF(NM.EQ.0.OR.QES(NM).NE.0) GO TO 515 ** DO 500 IS = 1,LS ** NM = SOLS(10,IS) ** CSAVE = 0. ** ITOT=ITM(NM) ** DO 505 I=1,ITOT ** CFM = COEFM(NM,I) ** SPM = SPECM(NM,I) ** CSAVE1 = CFM*DLOG10(ACT(SPM)) **505 CSAVE = CSAVE + CSAVE1 ** ACTM(IS) = CSAVE - LGKM(NM) **500 CONTINUE ** CALL ELSAT(ACTM(1),ACTM(2),TEMPK) ** ********** **515 CONTINUE ** NOW WE NEED TO COMPUTE GAS FUGACITIES, SATURATION PRESSURE AND ** FUGACTITY COEFFICIENTS, IF GASES ARE NOT SUPERSATURATED, ** BEFORE WE CALL MINSAT. WE TEST IF GASES ARE SUPERSATURATED ** (THAT MEANS IN THE MATRIX) ON H2O GAS. IF H2O GAS IS NOT IN ** THE MATRIX (QES(1)=0) THEN GASES ARE NOT SUPERSATURATED. ** WE REPEAT SOME OF THE ABOVE W/ SOLID SOL, BUT IT IS CLEANER ** IF (QES(1).NE.0.OR.IPSAT.EQ.0) GO TO 399 ** DO 438 NM=1,NGAS ** CSAVE = 0. ** ITOT=ITM(NM) ** DO 440 I=1,ITOT ** CFM = COEFM(NM,I) ** SPM = SPECM(NM,I) ** CSAVE1 = CFM*DLOG10(ACT(SPM)) **440 CSAVE = CSAVE + CSAVE1 ** GASLOG(NM) = CSAVE - AKGASL(NM) **438 CONTINUE ** CALL PSAT * ********* * IF(LOOC.EQ.3) GO TO 406 CALL MINSAT * *********** WRITE(IOUT1,944) IF(ERLOOP.LT.NLOOP) GO TO 406 WRITE(IOUT1,620) NLOOP WRITE(IOUT2,620) NLOOP 620 FORMAT(//5X,'NUMBER OF ITERATIONS REACHED THE GIVEN LIMIT (',I4, 1 ')',/5X,'NON-CONVERGENCE ASSUMED. EXECUTION ABORTED.'//) STOP * 406 IF(IPUNCH.EQ.0) GO TO 407 IF(BIT.NE.SINC.AND.LOOC.NE.3) GO TO 407 * * WRITE PLOTTING FILE ON OPTION * IF(IPUNCH.EQ.1.OR.IPUNCH.EQ.2) THEN * SYMONDS PLOT FILE WRITE(IPUN,832) WRLOG = 0. IF(TOTMIX.NE.0.) WRLOG = DLOG10(WR) ALA1=DLOG10(ACT(BSTP1)) DO 409 NS = 1,NST ALA = DLOG10(ACT(NS)) * PLOT LOG MOLES IF INPUNCH EQUALS 2 IF(IPUNCH.EQ.2) ALA = DLOG10(MTRY(NS)) WRITE(IPUN,832) NAME(NS),ALA,TEMP,PFLUID,ALA1,WRLOG 409 CONTINUE DO 413 I = 1, NMTT *+ XMIN=1. *+ JS=ISOLS(SEQ(I)) *+ IF(JS.NE.0) XMIN=MINTRY(SEQ(I))/SOLTRY(JS) *+ GRAMS = MINTRY(SEQ(I)) * XMOLW(SEQ(I)) *+ IF(GRT.NE.0.) PCMIN=GRAMS/GRT*100. ALA = MINTRY(I) IF(ALA.EQ.0.) ALA = 1.D-300 ALA = DLOG10(ALA) WRITE(IPUN,832) MIN(I),ALA,TEMP,PFLUID,ALA1,WRLOG 413 CONTINUE * ELSE IF(IPUNCH.EQ.3.OR.IPUNCH.EQ.4) THEN * * New plotting code inserted here, modified from CHILLER. This will * allow GASWORKS to run with MINPLOT plotting program. * * CONTROL WRITING OF TITLE AND NST, BSTOT and space holder (1) AT * TOP OF FIRST PLOT SEGMENT IF(IPUNCH.EQ.4) GO TO 415 IPUNCH = 4 WRITE(IPUN,800) BSTOT,NST,1 WRITE(IPUN,870) TITLE1 WRITE(IPUN,870) TITLE2 415 TOT = TOTMIX WRITE(IPUN,833) TEMP,PFLUID,ACT(BSTP1),TOT WRITE(IPUN,830) (MTRY(I), I=1,NST) WRITE(IPUN,830) (MTOT(I), I=1,BSTOT) WRITE(IPUN,830) (AQM(I), I=1,BSTOT) WRITE(IPUN,838) 1,'bogus', 0. DO 416 I = 1, NMT XMIN=1. JS=ISOLS(SEQ(I)) IF(JS.NE.0) XMIN=MINTRY(SEQ(I))/SOLTRY(JS) GRAMS = MINTRY(SEQ(I)) * XMOLW(SEQ(I)) IF(GRT.NE.0.) PCMIN=GRAMS/GRT*100. WRITE(IPUN,837) SEQ(I), MIN(SEQ(I)), MINTRY(SEQ(I)), 1 GRAMS, PCMIN, XMIN 416 CONTINUE WRITE (IPUN, 943) END IF * 407 IF(LOOC.EQ.3) STOP INP=IOUT3 IF(BIT.EQ.SINC) INP=INP1 * * A PICKUP FILE IS CREATED AFTER EACH EQUILIBRATION AND IS WRITTEN * IN THE FILE IOUT3 (E.G. EACH TIME WE PUT A NEW MINERAL IN, OR * KICK A MINERAL OUT). HOWEVER, WHEN WE FINALLY TRUELY EQUILIBRATE, * (NO MORE MINERALS ARE SUPER- OR UNDERSATURATED, AND WE ARE * READY FOR NEXT COOLING OR MIXING STEP), THE PICKUP FILE IS * WRITEN OVER THE LAST PICKUP FILE FOR THE LAST COMPLETED STEP. * THEREFORE, IOUT3 IS CREATED ONLY WHEN MINERALS ARE ADDED OR * OR KICKED OUT IN A GIVEN STEP. THE FOLLOWING CREATES PICKUP * FILES. REWIND INP WRITE(INP,870) TITLE1 WRITE(INP,870) TITLE2 WRITE(INP,801) WRITE(INP,834) ERPC,FO2LOG,PFLUID,TEMP,SINC,SLIM,GASTOT WRITE(INP,802) WRITE(INP,836) TEMPC,TOTMIX,SOLMIN, RM,GASGRM,SUPRNT WRITE(INP,810) WRITE(INP,800) IMINT,IFRAC,IPUNCH,NLOOP,ISTEP, 1 LIMSS,LOOC,ITREF,IDEAL,INCREM,INCP,MINSOL,OFLAG WRITE(INP,804) DO 520 I = 1,BSTOT 520 WRITE(INP,805) SAQ(I),NAME(I),MTOT(I),MTRY(I),GAMMA(I),COMTOT(I) WRITE(INP,806) IF(NMT.EQ.0) GO TO 525 WRITE(INP,823) (MIN(SEQ(I)), MINTRY(SEQ(I)), I=1,NMT) 525 WRITE(INP,807) IF(NOOX.EQ.0) GO TO 535 DO 530 I = 1, NOOX 530 WRITE(INP,823) NOMOX(I), WTPC(I) 535 WRITE(INP,808) DO 540 I = 1,NSUP 540 WRITE(INP,822) SUPNAM(I) WRITE(IOUT2,952) * IF(BIT.EQ.0.) GO TO 190 * TERMINAL OVERLOOK OF CALCULATIONS *% FO2=ACT(BSTP1) FO2LGP=DLOG10(ACT(BSTP1)) WRITE(IOUT2,841) FO2LGP 560 WRITE(IOUT2,840) TEMP,PFLUID IF(INCREM.NE.0) BIT = SINC * 598 CALL COOLER * *********** * * STOPS ARE BUILT IN COOLER WHENEVER SLIM IS REACHED * MIXER, BULK CALLED FROM COOLER * GO TO 190 * STOP IN COOLER * 800 FORMAT (16I5) 801 FORMAT(/'< ERPC >< LOGFO2 >< PFLUID >< TEMP >< SINC >' 1 ,'< SLIM >< GASTOT >') 802 FORMAT(/'< TEMPC >< TOTMIX >< SOLMIN >< RM >< GASGRM >' 1 ,'< SUPRNT >') 803 FORMAT(I5,1X,8X,4X,2D16.9,F8.4,2X,D10.4) 804 FORMAT(/,' saq> < name > < mtot >< mtry >' 1 ,' < comtot >') 805 FORMAT(I5,1X,A8,4X,2D16.9,F8.4,2X,D10.4) 806 FORMAT(/'< min > < mintry >') 807 FORMAT(/' < wtpc >') 808 FORMAT(/'') 810 FORMAT(/' bsto ifra ipun nloo iste lims looc itre idea' 1 ,' incr incp mins ofla') 815 FORMAT (//) 816 FORMAT (/,A80) 817 FORMAT (8X, 6D12.5) 818 FORMAT (//8X, 5D14.7) 820 FORMAT (A8,7X,A1,I1, 7(F6.3,I2),F5.3,I2) 822 FORMAT (A8, 5X,F9.4) 823 FORMAT (A8, 2X,D10.4) 830 FORMAT (8D10.4) 832 FORMAT(A8,5F10.5) 833 FORMAT(2F10.5,2D10.4,2F10.5) 834 FORMAT(D10.4,3F10.5,2G10.4,2F10.5) 836 FORMAT(F10.5,G10.4,D10.4,2F10.5,D10.4) 835 FORMAT(5D16.9) 837 FORMAT (I3,2X,A8,2X,7D10.5) 838 FORMAT(I3,2X,A8,2X,F10.5) 840 FORMAT(/1X,'CALCULATIONS DONE FOR ',F8.3,' DEG.C.;' 1 ,1X,F8.3,' ATM') 841 FORMAT(/5X,'THE LOG(FO2) IS: ',F7.2) 854 FORMAT (A8,A1,1X,F7.2,I2,1X,7(F6.3,I2),F4.3,I2,F6.3,I2,F6.3) 870 FORMAT (10A8) 885 FORMAT(5X,A8,' IS SUPPRESSED IN GASTHERM DATA FILE') 886 FORMAT(5X,A8,' IS SUPPRESSED IN THIS RUN') 900 FORMAT (/1X, 'MINERAL NO. LOG K TRIAL MOLES', 1' MINERAL NO. LOG K TRIAL MOLES'/) 905 FORMAT (2H ,'NO. OF LOOPS = ',I3,5X,'MOLES OF GAS=',D12.7) 906 FORMAT (2H , 'MOLES OF GAS=',D12.7) 910 FORMAT (2(1X, A8, 1X, I3, 3X, F8.3, 3X, D10.4, 4X)) 914 FORMAT (2H , ' MOLE ' 1,'PARTIAL') 915 FORMAT (2H , 'SPECIES MOLE NUMBER FRACTION ' 1,'PRESSURE FUGACITY LOG(F) GAMMA'/) 920 FORMAT (1X,I3,1X,A8,1X,4(D11.5,1X),F9.4,1X,F8.5) 923 FORMAT (/) 924 FORMAT (1X,I3, 1X, A8, 24X, 2(D11.5,1X,F10.4,1X), D11.5) 927 FORMAT (2H , 7X, 'TOTAL MOLES OF GAS = ',D12.7) 930 FORMAT (2H ,'MINERAL MOLES LOG MOLES GRAMS ', 1'LOG GRAMS WT PERCENT VOL (CM3)',/) 931 FORMAT (//2H ,'MINERAL MOLE FRACTION',6X,'ACTIVITY ' 1 ,3X,'ACTIVITY COEF',2X,'LOG ACTIVITY',/) 933 FORMAT (1H , 'TEMPERATURE = ',F8.3,' C.; P(FLUID) =',F8.3, 1' ATM; MIXER FRACTION = ',G10.4,/) 935 FORMAT (2X , A8 ,3X, 2(D10.4,2X,F7.3, 3X), 2(D10.4,3X)) 936 FORMAT(/2X,'TOTAL FRACTION OF MIXING SOLUTION ADDED: ',F12.6) 937 FORMAT (2X , A8 ,7X,2(D10.4,6X),F7.4,8X,F7.3) 940 FORMAT (2H ,'SPECIES NO. TOT.MOLES MTRY ', 1' GAMMA MIXING SOLUTION',/) 943 FORMAT (/) 944 FORMAT (/) 945 FORMAT (2H , A8, I5, 4X, 2D12.5, F10.4,3X,D10.4) 946 FORMAT(' NUMBER OF LOOPS USED = ',I4,' LOOP LIMIT = ',I4) 950 FORMAT (2H ,' SPECIES NO. LOG K', 1 ' SPECIES NO. LOG K',/) 952 FORMAT(/,5X,'PICKUP FILE WRITTEN ON DISK') 955 FORMAT (4X, A8, 3X, I4, 3X, F10.4, 11X, A8, 3X, I4, 3X, F10.4) 957 FORMAT (1X, ' COMPONENT MOLE NUMBER MOLE FRACTION', 1' PPM',/, 35X,'(TOTAL)'/) 960 FORMAT (2H , 10A8) 961 FORMAT (1X, ' COMPONENT TOTAL MOLES GAS PHASE MOLES', 1' SOLID MOLES'/) 962 FORMAT (1X, I2, 1X, A8, 2X, 3(D15.9, 2X)) 965 FORMAT(/16X,'WEIGHT OF GASEOUS SOLUTION: ',F12.6,' GRAMS'//) 968 FORMAT(2X,'IF UNKNOWN, MINERAL DENSITIES ARE ASSUMED TO BE 2.5'/) 970 FORMAT (2H , 78('+')) 971 FORMAT (2H ,'AMOUNT OF SOLID REACTANT CONSUMED:', F12.6,1X,A8 1 ,' LOG=',F9.4) 972 FORMAT (1H ,'AMOUNT OF SOLID PRODUCTS PRODUCED:', F12.6, 1' GRAMS;',1X,'VOLUME:',F14.8,' CM3') 973 FORMAT (2H ,'AMOUNT OF ORIGINAL SOLUTION: ',F12.6,' GRAMS' 1 /,2X,'ORIGINAL SOLUTION / ROCK ADDED = ',F14.4,' (= GAS/ROCK)'/) 974 FORMAT (/1H , 'INITIAL COMPOSITION OF SOLID PHASE (WT%):') 976 FORMAT (/1H , 'INITIAL COMPOSITION OF SOLID PHASE (MOLE%):') 975 FORMAT (/20(3(2X,A8,' = ',D10.4, ',')/)) 977 FORMAT (2H , 7D11.4) 999 FORMAT (1H1) END SUBROUTINE OXY IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DOUBLE PRECISION MTRY,MTOT,LGK,LGKM INTEGER BSTP1, BSTOT,OFLAG,SPEC COMMON/SLOUGH/ A(BG4,BG4),F(BG4),T(BG4,BG4),IV(BG4), N,ITREF COMMON/FROID/ SINC,BIT,SLIM,AKLOG,APFCL(BG2,5), 1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4), 2 XDH(4),ALK(5),LGK(BG2),TEMP,TEMPC, HEATC,COMTOT(BG3), HEATM, 3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3), 4 CVCCF(7,3), IFRAC COMMON /GEORGE/ AZERO(BG2),XI, NS, NST COMMON/MORT/MTRY(BG2),GAMMA(BG2),WATM,FO2,FO2LOG,TERMLN,KK 1,OFLAG COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,PLOG 1,GASTOT,NMTT,LOOC COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8) COMMON/BULL/BSTOT,BSTP1,JGSTOT COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH IF(ISTEP.EQ.1) WRITE(IOUT1,911) 911 FORMAT (2H , 40X, 'SUBROUTINE OXY WAS USED') * * * * PSEUDO SUBROUTINE FH2CALC. * NO PF/GASTOT TERM IN "F" EQUATION BECAUSE OF FO2LOG CANCELS IT. F(1) = -.5*LGK(BSTP1)+DLOG10(GAMMA(2)*MTRY(2)) 1 -DLOG10(GAMMA(1)*MTRY(1))-0.5*FO2LOG * DERIVATIVES W/R H2 AND H2O A(1,1) = -1./MTRY(1)*TERMLN A(1,2) = 1./MTRY(2)*TERMLN DO 57 J = 3,N 57 A(1,J) = 0.0 DO 100 NS = BSTP1,NST CFF=0. IF(SPEC(NS,2).NE.1) GO TO 100 CFF=COEF(NS,2) A(JGSTOT,1) = A(JGSTOT,1) + CFF*MTRY(NS)/MTRY(1) 100 CONTINUE RETURN END * SUBROUTINE MINERA IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DOUBLE PRECISION MSL ,MTRY, MINTRY,MTOT,LGKM CHARACTER*8 MIN INTEGER SOLS, QES,OFLAG INTEGER SPM, SPECM, BSTOT,BSTP1 INTEGER SPEC, SEQ * THIS SUBROUTINE SETS UP ALL PARTIAL DERIVATIVE MATRIX ELEMENTS FOR * MINERALS AND SOLID SOLUTIONS * EQUATIONS BSTP1 THROUGH BSTP1 + NMT ARE HYDROLYSIS Q EXPRESSIONS * IN LOG FORM E.G. F = -LOG K + LOG Q COMMON/MAXSIZ/ MAXSAQ,MAXNS,MAXNM,MAXSTO,MAXSSL,MAXN,MAXMEM,MAXIN COMMON/MORRIS/ MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW 1 (BG1), TEMPK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1) COMMON/TRIX/ TMIN(BG3) COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,PLOG 1,GASTOT,NMTT,LOOC COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8) COMMON/BULL/BSTOT,BSTP1,JGSTOT COMMON/SLOUGH/ A(BG4,BG4),F(BG4), T(BG4,BG4),IV(BG4), N,ITREF COMMON/MORT/MTRY(BG2),GAMMA(BG2),WATM,FO2,FO2LOG,TERMLN,KK 1,OFLAG COMMON/MORSEL/ GAMIN(BG1),B(BG1),SUMSOL(15),SOLTRY(15), 1 DGDN(BG3,20),SOLS(15,20),ISOLS(BG1),LIMSOL,JSOLS(BG1),IST(15) ** COMMON/FUGA/GASLOG(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3, ** 1 NGAS,IMIX,IDEAL,ITGAS COMMON/PRECIS/ RMIN, RLOGMN, RLOGMX * N = BSTOT + NMT + 1 MINCT = MINCT + 1 DO 76 I = 1, MAXSAQ 76 TMIN(I) = 0.0 DO 101 M = 1, NMT NM = SEQ(M) IF(NM.EQ.0) GO TO 105 NMA = M + BSTOT + 1 ITOT = ITM(NM) MSL = 0. DO 90 I = 1, ITOT SPM = SPECM(NM,I) CFM = COEFM(NM,I) * ADD MOLES OF MINERAL TO MASS BALANCE EQUATION F(SPM) = F(SPM) + CFM * MINTRY(NM) A(SPM,NMA) = CFM * COMPUTE Q FOR MINERAL CHARACTERISTIC REACTION * ACT CONTAINS FUGACITY OF GASES CHK=CFM * DLOG10(MTRY(SPM)*GAMMA(SPM)*PFLUID/GASTOT) MSL = MSL + CHK * FOR MASS ACTION EQS IN LOG FORM, DERIVATIVES ARE SAME FOR GASEOUS * AND AQUEOUS FORMS A(NMA,SPM) = CFM * TERMLN/MTRY(SPM) * DERIVATIVE OF MINERAL MASS ACTION EQUATION W/R TO GASTOT A(NMA,JGSTOT) = A(NMA,JGSTOT) - CFM*TERMLN/GASTOT * TMIN CONTAINS MOLES OF COMPONENTS IN SOLIDS 92 TMIN(SPM) = TMIN(SPM) + CFM * MINTRY(NM) 90 CONTINUE IF (MSL.LT.RLOGMN.OR.MSL.GT.RLOGMX) WRITE(IOUT1,201) MIN(NM),MSL IF(ISOLS(NM).EQ.0) GO TO 95 JU = ISOLS(NM) JS = JSOLS(NM) * COMPUTE Q FOR MINERAL SOLID SOLUTIONS. NEGATIVE MINTRY NOT ALLOWED * SEE MAIN. IF B IS USED, GAMIN MUST BE 1, IF GAMIN IS NOT 1, B MUST * BE 1, THUS WE CAN PUT ALL THE EXPRESSION INTO ONE DLOG10. ACM = DLOG10 (MINTRY(NM)/SOLTRY(JU) *GAMIN(NM))*B(NM) IF (ACM.LT.RLOGMN.OR.ACM.GT.RLOGMX) WRITE(IOUT1,202) MIN(NM),ACM, 1MINTRY(NM),SOLTRY(JU) MSL = MSL - ACM 95 CONTINUE * F(NMA) IS MASS ACTION EQUATION FOR MINERAL (IN LOG FORM) F(NMA) = -LGKM(NM) + MSL A(NMA,NMA) = 0.0 IF(ISOLS(NM).EQ.0) GO TO 101 * BELOW FOR SOLID SOLUTIONS ONLY * NOW WE COMPUTE PARTIAL DERIVATIVES OF MASS ACTION EQ. W/ REPECT * TO MOLES OF MINERAL FOR SOLID SOLUTIONS. A(NMA,NMA) = -B(NM)*TERMLN*(1./MINTRY(NM) - 1./SOLTRY(JU)) IS = IST(JS) DO 110 I = 1, IS LS = SOLS(JS,I) * HERE WE INCLUDE DERIVATIVE OF GAMIN(NM) W/ RESPECT TO MOLES OF * NM. THIS IS TO EASE CONVERGENCE IF NEEDED. DGDN ARE COMPUTED IN * THE ROUTINES THAT COMPUTE GAMIN. THESE CAN BE ZERO IN MOST CASES. IF(LS.NE.NM) GO TO 99 A(NMA,NMA) = A(NMA,NMA) - TERMLN/GAMIN(NM)*DGDN(QES(NM),I) GO TO 110 * 99 CONTINUE DO 111 KI = 1, NMT * WE SKIP IF END MEMBER NOT IN MATRIX IF(SEQ(KI).EQ.0) GO TO 112 IF(LS.EQ.SEQ(KI)) NMD = KI + BSTOT + 1 IF(LS.EQ.SEQ(KI)) A(NMA,NMD) = B(NM)*TERMLN/SOLTRY(JU) . - TERMLN/GAMIN(NM)*DGDN(QES(NM),I) * ABOVE LINE CONTAINS DERIVATIVE OF GAMIN(NM) W/ REPECT TO * MOLES OF ENDMEMBER DIFFERENT THAN NM. SEE COMMENTS ABOVE. 111 CONTINUE 112 CONTINUE 110 CONTINUE 101 CONTINUE 105 CONTINUE * 201 FORMAT (/1X,'LOG Q FOR ',A8,' IS ',F7.3) 202 FORMAT (1X,'LOG A FOR ',A8,'IS',F7.3/1X,'MINTRY=',D12.5,'SOLTRY= 1',D12.5) RETURN END SUBROUTINE MINSAT IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DOUBLE PRECISION MSL,MTRY,MINTRY, MTOT, LGKM ,LGK CHARACTER*8 MIN INTEGER SOLS, QES,OFLAG INTEGER SPM, SPECM, BSTOT,BSTP1,SPEC, SEQ * SUBROUTINE MINSAT DETERMINES WHETHER THERE IS SUPERSATURATION OF * THE PRESENT SOLUTION WITH ANY MINERALS OR SOLID SOLUTIONS. * IF MORE THAN ONE MINERAL IS SUPERSATURATED, IT ADDS THE ONE WITH * THE LOWEST AFFINITY TO THE MATRIX. IT ALSO REMOVES ANY MINERALS * PREVIOUSLY IN MATRIX THAT ARE NOW UNDERSATURATED DIMENSION NMSAV(BG3),XLGQM(BG1), MTOSS(5),XM(15,20) COMMON/MAXSIZ/ MAXSAQ,MAXNS,MAXNM,MAXSTO,MAXSSL,MAXN,MAXMEM,MAXIN COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8) COMMON/BULL/BSTOT,BSTP1,JGSTOT COMMON/SLOUGH/ A(BG4,BG4),F(BG4), T(BG4,BG4),IV(BG4), N,ITREF COMMON/MORT/MTRY(BG2),GAMMA(BG2),WATM,FO2,FO2LOG,TERMLN,KK 1,OFLAG COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,PLOG 1,GASTOT,NMTT,LOOC COMMON/FROID/ SINC,BIT,SLIM,AKLOG,APFCL(BG2,5), 1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4), 2 XDH(4),ALK(5),LGK(BG2),TEMP,TEMPC, HEATC,COMTOT(BG3), HEATM, 3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3), 4 CVCCF(7,3), IFRAC COMMON/MORRIS/ MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW 1 (BG1), TEMPK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1) COMMON/MORSEL/ GAMIN(BG1),B(BG1),SUMSOL(15),SOLTRY(15), 1 DGDN(BG3,20),SOLS(15,20),ISOLS(BG1),LIMSOL,JSOLS(BG1),IST(15) ** COMMON/FUGA/GASLOG(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3, ** 1 NGAS,IMIX,IDEAL,ITGAS COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH COMMON/PRECIS/ RMIN, RLOGMN, RLOGMX RT = 1.98719 * TEMPK AFOLD = 0. LREAP = 0 MCNT = 0 NEWCT = 0 NMT = 0 JACK=0 DO 57 I = 1, MAXIN 57 SEQ(I) = 0 DO 201 NM = 1, NMTT MSL = 0.0 ITOT = ITM(NM) * WE COMPUTE LOG Q DO 93 I = 1,ITOT CFM = COEFM(NM,I) SPM = SPECM(NM,I) CSAVE=CFM*DLOG10(ACT(SPM)) MSL = MSL + CSAVE 93 CONTINUE JS = JSOLS(NM) IF(JS.EQ.0) GO TO 96 GAMLOG = DLOG10(GAMIN(NM)) IF(ISOLS(NM).NE.0) GO TO 94 * * GASES OR MINERALS NOT YET SUPERSATURATED * COMPUTES SUM OF MOLE FRACTIONS FOR SOLID SOLUTIONS THAT ARE NOT * YET SUPERSATURATED. IF SUMM ADDS TO 1 OR GREATER, WE WILL SUPERSAT QTERM = (MSL - LGKM(NM))/B(NM) - GAMLOG IF (QTERM.LE.RLOGMN) QTERM=RLOGMN IF (QTERM.GE.RLOGMX) QTERM=RLOGMX QTERM =10.**QTERM SUMSOL(JS) = SUMSOL(JS) + QTERM IS = IST(JS) JACK=1 DO 73 J=1,IS IF(SOLS(JS,J).EQ.NM) XM(JS,J) = QTERM 73 CONTINUE GO TO 96 * * COMPUTE Q FOR SATURATED SOLID SOLUTIONS 94 JU = ISOLS(NM) ACM = DLOG10(DABS(MINTRY(NM)/SOLTRY(JU)))*B(NM) + GAMLOG MSL = MSL - ACM * 96 XLGQM(NM) = MSL DIFM = XLGQM(NM) - LGKM(NM) * REMOVE UNDERSATURATED MINERALS FROM MATRIX * ================================================ * IF SINGLE MINERAL WAS PREVIOUSLY IN MATRIX BUT GOT UNDERSATURATED * (NEGATIVE AMOUNT OF MOLES WERE COMPUTED) WE REMOVE IT IF(QES(NM).EQ.0) GO TO 97 IF(MINTRY(NM).LT.0.) GO TO 195 * * IN CASE MASS BALANCE CONVERGES, BUT MASS ACTION FAILS (SHOULD NOT * HAPPEN WITH ERPC SUFFICIENTLY LOW), THE AFFINITY (-RTLN Q/K ) * WILL BE TOO LARGE => MINERAL MIGHT NOT BELONG, OR PROBLEM. WE * KICK MINERAL OUT (USEFUL TO DUMP TOO SMALL SOL-SOL ENDMEMBERS) AF = -RT*DIFM*2.303 IF(DABS(AF).LT.1.D-6) GO TO 101 WRITE(IOUT1,800) MIN(NM), AF 800 FORMAT(/5X,A8,' DOES NOT SATISFY ANYMORE THE MASS ACTION EQUATION' . ,/13X,' WE REMOVE IT FROM MATRIX. AFFINITY =', D10.4) GO TO 195 * * CHECK FOR SUPERSATURATION OF SINGLE MINERALS * ============================================ 97 CONTINUE DIFM=DIFM/SCALE(NM) * IF UNDERSATURATED, GO TO 201 IF(DIFM.LT.1.D-10) GO TO 201 * NOW, WE SELECT ONLY ONE MINERAL WITH SMALLEST AFFINITY IF * SEVERAL MINERALS ARE SUPERSATURATED IF(DIFM.LT.AFOLD) GO TO 99 AFOLD = DIFM MSTORE = NM 99 NEWCT = NEWCT + 1 GO TO 201 * * 101-195, RE-COUNT SATURATED MINERALS 101 MCNT = MCNT + 1 NMSAV(MCNT) = NM SEQ(MCNT) = NM QES(NM) = MCNT NMT = MCNT GO TO 201 * * HERE WE KICK MINERALS OUT OF THE MATRIX 195 LREAP = LREAP + 1 MASH = NM MTOSS(LREAP)= NM QES(NM) = 0 MINTRY(NM) = 1.D-300 * 201 CONTINUE * IF MINERALS WERE REMOVED FROM MATRIX IF(LREAP.NE.0) GO TO 107 * IF NO NEW MINERALS SUPERSATURATE, AFOLD = 0 IF(AFOLD.EQ.0.) GO TO 207 * MORE THAN ONE NEW SUPERSAT MINERAL: TAKE THE MOST SUPERASTURATED IF(LOOC.EQ.0) GO TO 105 * MORE THAN ONE NEW SUPERSAT MINERAL: BACK TITRATE IF(NEWCT.GE.2) GO TO 107 * * WE PUT IN THE MOST SUPERSAT. MINERAL (INDEX = MSTORE) 105 MCNT = MCNT + 1 NM = MSTORE SEQ(MCNT) = MSTORE QES(MSTORE) = MCNT NMSAV(MCNT) = MSTORE NMT = MCNT IF(JSOLS(NM).EQ.0) GO TO 107 * IF WE JUST PUT IN A MINERAL WHICH BELONG TO A SOLID SOLUTION, THEN * WE NEED TO PUT IN THE REMAINING END MEMBERS. JS = JSOLS(NM) IS = IST(JS) MINTRY(NM) = 1. D-7 DO 106 J = 1, IS M = SOLS(JS,J) * IF OTHER ENDMEMBERS WERE ALREADY IN MATRIX, WE SKIP IF(QES(M).NE.0) GO TO 106 IF(M.NE.NM) MINTRY(M) = MINTRY(NM)*XM(JS,J)/SUMSOL(JS) IF(M.EQ.NM) GO TO 106 MCNT = MCNT + 1 SEQ(MCNT) = M QES(M) = MCNT NMSAV(MCNT) = M NMT = MCNT 106 CONTINUE GO TO 107 * * NOW IN CASE A SOLID SOLUTION, WITH ALL END MEMBERS * INDEPENDENTLY UNDERSATURATED, SUPERSATURATES, WE PUT IT IN 207 JSAV = 0 IF(JACK.EQ.0) GO TO 107 SOLQ = 1.0 DO 213 NM = 1, NMTT JS = JSOLS(NM) IF(JS.EQ.0) GO TO 213 IS = IST(JS) DO 210 J = 1, IS JK = SOLS(JS,J) IF(QES(JK).NE.0) GO TO 213 210 CONTINUE IF(SUMSOL(JS).LT.SOLQ) GO TO 213 SOLQ = SUMSOL(JS) JSAV = JS 213 CONTINUE IF(JSAV.EQ.0) GO TO 107 NEWCT = NEWCT + 1 IS = IST(JSAV) DO 225 J = 1, IS M = SOLS(JSAV,J) MCNT = MCNT + 1 SEQ(MCNT) = M QES(M) = MCNT NMSAV(MCNT) = M MINTRY(M) = 1.D-8*XM(JSAV,J)/SUMSOL(JSAV) NMT = MCNT 225 CONTINUE * 107 CONTINUE LBOTH = LREAP + NEWCT IF(LBOTH.GE.2) GO TO 233 IF(LREAP.NE.1) GO TO 311 MINTRY(MASH) = 1.D-300 GO TO 311 233 IF(LOOC.EQ.0) GO TO 311 IF(BIT.EQ.0.) BIT = SINC BIT =-(BIT/2.) * IF MORE THAN TWO MINERALS ARE UNDER OR SUPER SATURATED (LBOTH.GE. * 2) THEN MUST REPLACE UNDERSAT ONES AND BACK UP IF(LREAP.EQ.0) GO TO 108 DO 237 I = 1, LREAP NM = MTOSS(I) MCNT = MCNT + 1 SEQ(MCNT) = NM QES(NM) = MCNT NMSAV(MCNT) = NM 237 NMT = MCNT GO TO 108 * 311 IF(LBOTH.NE.0) BIT = 0. IF(LBOTH.EQ.0) BIT = SINC 108 CONTINUE CALL SQUIMP WRITE(IOUT1,923) WRITE(IOUT1,970) WRITE(IOUT1,923) IF(MCNT.EQ.0) GO TO 110 WRITE(IOUT1,975) WRITE(IOUT1,980) DO 109 I = 1, MCNT NM = NMSAV(I) XLK = LGKM(NM) AF = (XLK - XLGQM(NM)) * RT * 2.303 QKLOG = XLGQM(NM) - XLK QKSLOG=QKLOG/SCALE(NM) WRITE(IOUT1,977) NM,MIN(NM),XLK,XLGQM(NM),QKLOG,QKSLOG,AF XLGQM(NM) = 666. 109 CONTINUE WRITE(IOUT1,923) 110 WRITE(IOUT1,979) WRITE(IOUT1,980) DO 112 NM = 1, NMTT IF(XLGQM(NM).EQ.666.) GO TO 112 XLK =LGKM(NM) AF = (XLK - XLGQM(NM)) * RT * 2.303 QKLOG = XLGQM(NM) - XLK QKSLOG=QKLOG/SCALE(NM) ** IF(NM.GT.ITGAS) GO TO 130 ** GO TO 112 130 IF(QKLOG.GE.-5.) 1 WRITE(IOUT1,977) NM,MIN(NM),XLK,XLGQM(NM),QKLOG,QKSLOG,AF 112 CONTINUE * 923 FORMAT (/) 970 FORMAT (1H , 78('+')) 975 FORMAT (2H , 'THE FOLLOWING MINERALS ARE PRESENTLY IN MATRIX OR W 1ERE JUST ADDED',//) 977 FORMAT (2X,I3,2X,A8,2X,3(3X,F8.2),3X,F8.3,4X,D13.5,4X,F8.4) 979 FORMAT (2H , ' THE FOLLOWING MINERALS ARE PRESENTLY EXCLUDED FROM 1 MATRIX',/3X,'MINERALS WITH LOG(Q/K) LESS THAN -5 ARE NOT LISTED' 2 ,' BELOW'//) 980 FORMAT (5X,' MINERAL LOG K LOG Q LOG(Q/K)', 1' LOG(Q/K)/S AFFINITY',/) RETURN END SUBROUTINE SQUOMP * TO SORT AND INITIALIZE SOLID SOLUTION ARRAYS. * NM MINERAL INDEX: CURRENT SEQUENCE POSITION AS READ FROM * SOLTHERM FOR EACH PARTICULAR RUN. * JSOLS(NM) SOLID SOLUTION INDEX: FIXED ID NUMBER (E.G.GASES = 6) * B(NM) COEFFICIENT FOR MULTISITE MIXING * SOLS(JS,L) RETURNS (CONTAINS) NM VALUE, WITH JS = JSOLS(NM) AND * L = POSITION OF END ENDMEMBER FOR THAT JS SOLID SOLUTION * (1ST, 2ND ETC. END MEMBER) * IST(JS) TOTAL NUMBER OF END MEMBERS IN SOL SOL JS = JSOLS(NM) * ISOLS(NM) RETURNS (CONTAINS) THE POSITION, IN ARRAY OF CURRENTLY * SATURATED SOL SOL, OF THE SOLID SOLUTION CONTAINING NM * SOLTRY(J) SUM OF MOLES OF ALL END MEMBERS IN SATURATED SOL SOL * J = ISOLS(NM) * SUMSOL(J) SUM OF Q/K'S OF ALL END MEMBERS IN SATURATED SOL SOL * J = ISOLS(NM) * DGDN(J,I) DERIVATIVE OF GAMIN(NM) W/ RESPECT TO MOLES OF * ENDMEMBER I, WITH J = QES(NM) AND I = L ABOVE * IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DOUBLE PRECISION LGKM,MINTRY CHARACTER*8 MIN INTEGER SPECM, SEQ, QES, SOLS COMMON/MAXSIZ/ MAXSAQ,MAXNS,MAXNM,MAXSTO,MAXSSL,MAXN,MAXMEM,MAXIN COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,PLOG 1,GASTOT,NMTT,LOOC COMMON/MORRIS/ MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW 1 (BG1), TEMPK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1) COMMON/MORSEL/ GAMIN(BG1),B(BG1),SUMSOL(15),SOLTRY(15), 1 DGDN(BG3,20),SOLS(15,20),ISOLS(BG1),LIMSOL,JSOLS(BG1),IST(15) COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH WRITE(IOUT1,500) DO 180 J = 1,MAXIN DO 180 I = 1,MAXMEM 180 DGDN(J,I) = 0. DO 183 J = 1,MAXSSL DO 182 I = 1,MAXMEM 182 SOLS(J,I) = 0 IST(J) = 0 183 CONTINUE DO 186 NM = 1, NMTT JS = JSOLS(NM) IF(JS.EQ.0) GO TO 186 L = 0 IF(JS.LE.MAXSSL) GO TO 185 WRITE(IOUT1,400) MIN(NM), JS, MAXSSL 400 FORMAT(/5X,A8,' HAS A TOO LARGE SOLID SOLUTION INDEX (',I3,')', $ /5X,'MAXIMUM VALUE ALLOWED IS ',I3//) STOP * 185 L = L + 1 IF(L.LE.MAXMEM) GO TO 187 WRITE(IOUT1,410) MIN(NM), L, MAXMEM 410 FORMAT(/5X,A8,' IS PART OF A SOLID SOLUTION WITH TOO MANY' $ ,' END MEMBERS (',I3,')'/5X,'MAXIMUM NUMBER ALLOWED IS ',I3//) STOP * 187 IF(SOLS(JS,L).NE.0) GO TO 185 SOLS(JS,L) = NM IST(JS) = L WRITE(IOUT1,510) MIN(NM), JS, B(NM) 186 CONTINUE 500 FORMAT(/6X,'THE FOLLOWING MINERALS AND/OR GASES WITH IDENTICAL', 1/6X,'INDEX ARE IN SOLID SOLUTION:'//,8X, 2 'MINERAL INDEX B(NM)') 510 FORMAT(8X,A8,3X,I5,3X,F6.3) * ENTRY SQUIMP DO 163 NM = 1, NMTT 163 ISOLS(NM) = 0 L = 0 DO 189 I = 1,MAXSSL * RESET SUMSOL TO ZERO FOR MINSAT SUMSOL(I) = 0. M1 = SOLS(I,1) IF(IST(I).GE.2) GO TO 188 IF(IST(I).EQ.0) GO TO 189 SOLS(I,1) = 0 JSOLS(M1) = 0 IST(I) = 0 GO TO 189 188 CONTINUE IS = IST(I) L = L + 1 SOLTRY(L) = 0.0 DO 196 J = 1, IS LS = SOLS(I,J) * SKIP IF ONE END MEMBER IS NOT IN MATRIX IF(QES(LS).EQ.0) GO TO 196 SOLTRY(L) = SOLTRY(L) + MINTRY(LS) ISOLS(LS) = L 196 CONTINUE 189 CONTINUE LIMSOL = L RETURN END SUBROUTINE COOLER * COOLER IS USED TO GO ON TO THE NEXT STEP, WHICH CAN BE: * T CHANGE, P CHANGE, P-T CHANGE, OR COMPOSITION CHANGE * NOTE, CURRENTLY P CHANGE ONLY AFFECTS GASES * IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DOUBLE PRECISION MINTRY,MTRY,MTOT,LGKM,LGK CHARACTER*8 MIN,NOMOX INTEGER QES,SEQ,QAS,SAQ,SPECM,BSTOT,SPEC,SAQTOT,SPECOX,BSTP1 1,OFLAG COMMON /GEORGE/AZERO(BG2),XI, NS, NST COMMON/TRIX/ TMIN(BG3) COMMON/MORT/MTRY(BG2),GAMMA(BG2),WATM,FO2,FO2LOG,TERMLN,KK 1,OFLAG COMMON/MORRIS/ MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW 1 (BG1), TEMPK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1) COMMON/TURNIP/ SAQ(BG3), QAS(BG3), SAQTOT COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8) COMMON/BULL/BSTOT,BSTP1,JGSTOT COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,PLOG 1,GASTOT,NMTT,LOOC COMMON/FROID/ SINC,BIT,SLIM,AKLOG,APFCL(BG2,5), 1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4), 2 XDH(4),ALK(5),LGK(BG2),TEMP,TEMPC, HEATC,COMTOT(BG3), HEATM, 3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3), 4 CVCCF(7,3), IFRAC ** COMMON/FUGA/GASLOG(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3, ** 1 NGAS,IMIX,IDEAL,ITGAS COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH COMMON/PP/TS(200),PS(200),NPP,INCP COMMON/BESS/ WTPC(BG3), COEFOX(BG3,8), XMW(BG3), 1 NOMOX(BG3),ITOX(BG3),NOOX, SPECOX(BG3,8),MINSOL * * * DECREASE P AT CONSTANT T 65 IF(INCP.EQ.0) GO TO 80 PFLUID = PFLUID + BIT WRITE(IOUT1,904) BIT , PFLUID WRITE(IOUT2,846) BIT PLOG = DLOG10(PFLUID) IF(SINC.GT.0..AND.PFLUID.LE.SLIM) GO TO 160 IF(SINC.LT.0..AND.PFLUID.GE.SLIM) GO TO 160 WRITE(IOUT1,628) SLIM, PFLUID WRITE(IOUT2,628) SLIM, PFLUID STOP * * TITRATE SOLIDS AT CONSTANT T 80 IF(MINSOL.EQ.0) GO TO 85 CALL BULK IF(SINC.GT.0..AND.TOTMIX.LE.SLIM) GO TO 160 IF(SINC.LT.0..AND.TOTMIX.GE.SLIM) GO TO 160 WRITE(IOUT1,627) SLIM, TOTMIX WRITE(IOUT2,627) SLIM, TOTMIX STOP * * MIX THE SOLUTION, WITH OR WITHOUT T CHANGE 85 IF(TEMPC.EQ.0.) GO TO 90 CALL MIXER IF(SINC.GT.0..AND.TOTMIX.GT.SLIM) GO TO 88 IF(SINC.LT.0..AND.TOTMIX.LT.SLIM) GO TO 88 IF(TEMPC.EQ.TEMP) GO TO 160 GO TO 100 88 WRITE(IOUT1,627) SLIM, TOTMIX WRITE(IOUT2,627) SLIM, TOTMIX STOP * * CHANGE TEMPERATURE 90 TEMP = TEMP + BIT WRITE(IOUT1,903) BIT , TEMP WRITE(IOUT2,842) BIT IF(TEMP.LE.1227.) GO TO 92 WRITE(IOUT1,625) TEMP WRITE(IOUT2,625) TEMP STOP 92 IF (OFLAG.NE.0) GO TO 91 FO2 = 0. FO2LOG = 0. 91 IF(SINC.GT.0..AND.TEMP.LE.SLIM) GO TO 100 IF(SINC.LT.0..AND.TEMP.GE.SLIM) GO TO 100 WRITE(IOUT1,626) SLIM,TEMP WRITE(IOUT2,626) SLIM, TEMP STOP * * BELOW WE RECALCULATE ALL THERMO DATA IF T OR P WAS CHANGED 100 CONTINUE TEMPK=TEMP + 273.15 ** DO 50 I=1,3 ** 50 CVB(I) = CVBCF(I,1)/TEMPK/TEMPK + CVBCF(I,2)/TEMPK + CVBCF(I,3) ** DO 52 I=1,7 ** 52 CVC(I) = CVCCF(I,1)/TEMPK/TEMPK + CVCCF(I,2)/TEMPK + CVCCF(I,3) ** DO 54 I=1,NGAS ** VB(I)= VBCOE(I,1)/TEMPK/TEMPK + VBCOE(I,2)/TEMPK + VBCOE(I,3) ** 54 VC(I)= VCCOE(I,1)/TEMPK/TEMPK + VCCOE(I,2)/TEMPK + VCCOE(I,3) ** ** DO 127 J = 1, 9 ** DO 123 I = 1,5 **123 COE(I) = XPFC(J,I) ** CALL POWER(COE,TEMP,AKLOG) ** XF(J) = AKLOG 127 CONTINUE * RECOMPUTE DERIVED SPECIES DATA WITH NEW T NBSP1 = SAQTOT + 1 DO 107 NS = NBSP1, NST DO 109 I = 1, 5 109 COE(I) = APFCL(NS,I) 110 CALL POWERK(COE,TEMPK,AKLOG) LGK(NS) = AKLOG 107 CONTINUE * * THERE ARE TWO WAYS CHANGE P AS F(T): POWER FUNCTION P AS * F(T), WHOSE COEFFICIENTS ARE READ IN GASTHERM; A SET OF GIVEN * P-T POINTS READ FROM AN ACCESSORY FILE (INP3). * IF NEITHER OF THESE OPTIONS ARE USED, P IS CONSTANT. IF(XF(9).NE.0.) PFLUID=XF(9) IF(NPP.NE.0) CALL PRESS(PFLUID,TEMP) PLOG = DLOG10(PFLUID) * HERE WE RECOMPUTE GAS AND MINERAL DATA WITH NEW P AND T DO 117 NM = 1, NMTT DO 111 I = 1, 5 111 COE(I) = BPFCL(NM,I) 114 CALL POWERK(COE,TEMPK,AKLOG) LGKM(NM) = AKLOG 117 CONTINUE * * BIT NOT EQUAL TO SINC ONLY IF BACK TITRATION 160 IF(BIT.NE.SINC) GO TO 210 * FRACTIONATION OPTIONS BELOW * NO FRACTIONATON IF(IFRAC.GT.0) GO TO 202 WRITE(IOUT1,907) GO TO 210 * * FRACTIONATES SOLIDS ONLY 202 DO 170 NS = 1, SAQTOT 170 MTOT(NS) = MTOT(NS) - TMIN(NS) WRITE(IOUT1,908) WRITE(IOUT2,908) GO TO 210 * 210 CONTINUE * NEXT TWO STATEMENTS CAUSE MINSAT TO USE THE LOWEST K/Q FOR * MINERAL SELECTION RATHER THAN T-CHANGE TO FIND SINGLE * SUPERSATURATED MINERAL. IF(DABS(BIT).LT..05.AND.(BIT*SINC).LT.0.) LOOC=0 IF (DABS(BIT/SINC).LT..24) LOOC=0 RETURN 625 FORMAT(///5X,'TEMPERATURE GREATER THAN 1227 C. TEMP =',F8.3, * /5X,'EXECUTION STOPPED.'//) 626 FORMAT(///5X,'TEMPERATURE REACHED GIVEN LIMIT OF ',F8.3, * /5X,'TEMP = ',F8.3,' EXECUTION STOPPED.'//) 627 FORMAT(///5X,'MIXING FRACTION REACHED GIVEN LIMIT OF ',G10.4, * /5X,'TOTMIX = ',G10.4,' EXECUTION STOPPED.'//) 628 FORMAT(///5X,'PRESSURE REACHED GIVEN LIMIT OF ',F8.3, * /5X,'PFLUID = ',F8.3,' EXECUTION STOPPED.'//) 630 FORMAT(///5X,'ENTHALPY REACHED GIVEN LIMIT OF ',F10.4, * /5X,'ENTH =',F10.4,' EXECUTION STOPPED.'//) * 842 FORMAT(5X,' NEXT TEMPERATURE CHANGE OF ',F8.3,'DEG.C.',///) 846 FORMAT(5X,' NEXT PRESSURE CHANGE OF ',F8.3,' ATM',///) 847 FORMAT(5X,' NEXT ENTHALPY CHANGE OF ',F10.4,'KCAL',/,' TOTAL ', 1' ENTHALPY ADDED (INCL. STARTING ENTHALPY)',F10.4,' KCAL'///) 903 FORMAT (' TEMPERATURE CHANGED BY ',F8.3,' DEGREES TO ',F8.3,' DE 1GREES C.',/) 904 FORMAT (' PRESSURE CHANGED BY ',F8.3,' ATM TO ',F8.3,' ATM'/) 905 FORMAT (' ENTHALPY CHANGED BY ',F10.4,' KCAL TO ',F10.4,' KCAL'/, 1' TOTAL ENTHALPY ADDED (INCL. STARTING ENTHALPY)',F10.4,' KCAL'/) 907 FORMAT (5X,'SOLIDS AND GASES CARRIED') 908 FORMAT(5X,'SOLIDS ONLY LEFT BEHIND') 909 FORMAT(5X,'GASES ONLY LEFT BEHIND') 910 FORMAT(5X,'SOLIDS AND GASES LEFT BEHIND') END SUBROUTINE MIXER IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DOUBLE PRECISION MTOT,LGK,LGKM INTEGER BSTOT,SPEC,BSTP1 COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8) COMMON/BULL/BSTOT,BSTP1,JGSTOT COMMON/FROID/ SINC,BIT,SLIM,AKLOG,APFCL(BG2,5), 1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4), 2 XDH(4),ALK(5),LGK(BG2),TEMP,TEMPC, HEATC,COMTOT(BG3), HEATM, 3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3), 4 CVCCF(7,3), IFRAC COMMON/MINNIE/ AQM(BG3),ACT(BG2),LGKM(BG1),SCALE(BG1),PFLUID,PLOG 1,GASTOT,NMTT,LOOC COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH * * TITRATE MIXER (COOL) SOLUTION INTO ORIGINAL (HOT) SOLUTION WRITE(IOUT1,925) TOTMIX WRITE(IOUT2,925) TOTMIX DO 101 I = 1, BSTOT 101 MTOT(I) = MTOT(I) + BIT * COMTOT(I) IF (TEMP.EQ.TEMPC) GO TO 102 TEMP = (TEMP*(1+TOTMIX)+TEMPC*BIT)/(1+TOTMIX+BIT) * TEMPERATURE IS CURRENTLY A SIMPLE WEIGHTED SUM OF THE MIXED * COMONENTS. NO ACCOUNT FOR NON-LINEARITIES IN HEAT CAPACITIES * IS TAKEN. 102 TOTMIX=TOTMIX+BIT * WRITE(IOUT2,920) TOTMIX,TEMP WRITE(IOUT1,920) TOTMIX,TEMP RETURN * 920 FORMAT(5X,'NEW MIXER FRACTION: ',G10.4,/, 1 5X,'RESULTING TEMPERATURE: ',F8.3,' DEG.C.'/) 925 FORMAT(/5X,'MIXER FRACTION: ',G10.4,/) END SUBROUTINE BULK * BULK COMPUTES COMPOSITION OF SOLID PHASE IN TERMS OF COMPONENT * SPECIES, GIVEN WT. PERCENTS OF OXIDES, ELEMENTS, OR MINERALS. * IF MINERAL REACTANTS ARE USED, THE MINERAL NAMES SHOULD BE WRITTEN * TO MATCH THE FORM USED IN GASTHERM IF YOU WANT MINSOL = 2 * OPTION TO WORK * IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DOUBLE PRECISION MTOT,MINTRY,LGK INTEGER SPECM, SEQ, QES, SPEC, SPECOX, BSTOT,BSTP1 CHARACTER*8 AMOLE,MIN,NOMOX,TYP DIMENSION SOLMOL(BG3) COMMON/MORRIS/ MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW 1 (BG1), TEMPK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1) COMMON/OXEN/ MTOT(BG3),COEF(BG2,8),SPEC(BG2,8) COMMON/BULL/BSTOT,BSTP1,JGSTOT COMMON/FROID/ SINC,BIT,SLIM,AKLOG,APFCL(BG2,5), 1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4), 2 XDH(4),ALK(5),LGK(BG2),TEMP,TEMPC, HEATC,COMTOT(BG3), HEATM, 3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3), 4 CVCCF(7,3), IFRAC COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH COMMON/BESS/ WTPC(BG3), COEFOX(BG3,8), XMW(BG3), 1 NOMOX(BG3),ITOX(BG3),NOOX, SPECOX(BG3,8),MINSOL COMMON /GEORGE/ AZERO(BG2),XI, NS, NST * DATA TYP/'GRAMS '/, AMOLE/'MOLES '/ * OLDTOT = TOTMIX TOTMIX = TOTMIX + BIT DO 102 I = 1, BSTOT 102 SOLMOL(I) = 0.0 DO 105 I = 1, NOOX IF(XMW(I).EQ.0.) WRITE(IOUT1,600) NOMOX(I) 600 FORMAT(//5X,'WARNING, MOLECULAR WEIGHT OF REACTANT ',A8,' IS 0.'/ 1 5X,'CHECK MINOX DATA FILE FOR ERRORS OR DOUBLE ENTRIES.'//) OXMOL = BIT * WTPC(I) * .01/XMW(I) IF(MINSOL.GT.2) OXMOL = BIT*WTPC(I)* 0.01 * MINSOL = 1 TITRATION IN GRAMS (SINC IN GRAMS, WTPC IN WT%) * MINSOL = 2 ACTIVATES NON-TITRATION OF EQUILIBRATED REACTANTS * WITH GRAM TITRATION OPTION * MINSOL = 3 TITRATION OF IN MOLES (SINC IN MOLES, WTPC IN MOLE%) * MINSOL = 4 ACTIVATES NON-TITRATION OF EQUIL. REACTANTS W/ MOLE * TITRATION OPTION IF(MINSOL.NE.2.AND.MINSOL.NE.4) GO TO 95 DO 90 J = 1, NMT NM = SEQ(J) IF(NOMOX(I).NE.MIN(NM)) GO TO 90 OXMOL = 0.0 WRITE(IOUT1,900) MIN(NM) WRITE(IOUT2,900) MIN(NM) GO TO 95 90 CONTINUE 900 FORMAT(2H ,'REACTANT ',A8,' NOT TITRATED BECAUSE EQUILIBRATED') * 95 CONTINUE ITOT = ITOX(I) DO 103 IX = 1, ITOT CX = COEFOX(I,IX) JSPX = SPECOX(I,IX) IF(JSPX.EQ.0) GO TO 105 SOLMOL(JSPX) = SOLMOL(JSPX) + CX * OXMOL 103 CONTINUE 105 CONTINUE * WE TITRATE SOLID INTO SOLUTION AND RECALCULATE MASS BALANCE DO 205 I = 1,BSTOT MTOT(I) = MTOT(I) + SOLMOL(I) 205 CONTINUE * IF(MINSOL.GT.2) TYP = AMOLE WRITE(IOUT1,910) BIT, TYP WRITE(IOUT2,910) BIT, TYP WRITE(IOUT2,920) OLDTOT,TYP,TOTMIX,TYP WRITE(IOUT1,920) OLDTOT,TYP,TOTMIX,TYP RETURN * 910 FORMAT(/5X,'NEXT SOLID TITRATION INCREMENT OF:'3x,F15.7,1X,A8) 920 FORMAT(5X,'CURRENT TOTAL AMOUNT OF SOLID ADDED: ',F15.7,1X,A8 1 /5X,'NEW TOTAL AMOUNT OF SOLID ADDED: ',4x,F15.7,1x,A8//) END * SUBROUTINE POWER(COE,TEMP, AKLOG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DIMENSION COE(5) COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH COMMON/PRECIS/ RMIN, RLOGMN, RLOGMX AKLOG = COE(1) DO 107 I = 1, 4 J = I + 1 AKLOG = AKLOG + COE(J) * TEMP**I 107 CONTINUE RETURN END * SUBROUTINE POWERK(COE,TEMPK,AKLOG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DIMENSION COE(5) COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH COMMON/PRECIS/ RMIN, RLOGMN, RLOGMX AKLOG = COE(1) AKLOG = AKLOG + COE(2)/TEMPK AKLOG = AKLOG + COE(3)*TEMPK AKLOG = AKLOG + COE(4)/TEMPK**2 AKLOG = AKLOG + COE(5)*DLOG10(TEMPK) RETURN END SUBROUTINE SLUD IMPLICIT DOUBLE PRECISION (A-G,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) COMMON/SLOUGH/ A(BG4,BG4),F(BG4),T(BG4,BG4),IV(BG4),N,ITREF * DO 8100 J=1,N DO 8100 I=1,N T(I,J) = A(I,J) 8100 CONTINUE 8110 CONTINUE * IV(N) = 1 DO 8260 K=1,N PIV = DABS(T(K,K)) IF(K.GE.N) GO TO 8260 * L = K KP1 = K+1 DO 8200 I=KP1,N IF(PIV.GE.DABS(T(I,K))) GO TO 8200 PIV = DABS(T(I,K)) L = I 8200 CONTINUE IV(K) = L * TMP = T(L,K) IF(K.NE.L) GO TO 8210 IF(PIV.GT.0.0) GO TO 8220 IV(N) = 0 GO TO 8260 8210 CONTINUE IV(N) = -IV(N) T(L,K) = T(K,K) T(K,K) = TMP 8220 CONTINUE * TMP = -TMP DO 8230 I=KP1,N T(I,K) = T(I,K)/TMP 8230 CONTINUE * DO 8250 J=KP1,N TMP = T(L,J) T(L,J) = T(K,J) T(K,J) = TMP DO 8240 I=KP1,N T(I,J) = T(I,J) + T(I,K) * TMP 8240 CONTINUE 8250 CONTINUE 8260 CONTINUE IF(PIV.EQ.0.0) IV(N) = 0 RETURN END SUBROUTINE SIR(X) IMPLICIT DOUBLE PRECISION (A-G,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DIMENSION X(BG4), B(BG4), R(BG4) COMMON/SLOUGH/ A(BG4,BG4),F(BG4),T(BG4,BG4),IV(BG4),N,ITREF COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH * * B IS THE VECTOR CONTAINING THE RIGHT-HAND SIDE OF THE * SYSTEM OF LINEAR EQUATIONS WHICH IS -F. NOTE THE MINUS SIGN WHICH * IS COMPATIBLE WITH XNEW = XOLD + X IN MAIN (NOTE PLUS SIGN) * DO 150 I=1,N B(I)=-F(I) 150 CONTINUE * IER = 0 DO 8100 I=1,N X(I) = B(I) 8100 CONTINUE IF(IV(N).NE.0) GO TO 8500 * WRITE(IOUT1,915) 915 FORMAT(/,5X,'******** WARNING ********',/,2X, 1'COMPUTATIONAL SINGULARITY, IV(N) = 0.0',/) STOP * 8500 CONTINUE QUOT = 1.0 * CALL SBS(X) IF(ITREF.EQ.0) RETURN 8200 CALL SAXMB(X,B,R) IF(ITREF.EQ.1) RETURN CALL SBS(R) IER = IER+1 IF(ISTEP.EQ.1) WRITE(IOUT1,909) IER 909 FORMAT(5X,' RESIDUAL VECTOR (IER =',I2,')') IF(ISTEP.EQ.1) WRITE(IOUT1,973) (R(I), I=1,N) * XN1 = 0.0 RN1 = 0.0 DO 8210 I=1,N XN1 = XN1 + DABS(X(I)) RN1 = RN1 + DABS(R(I)) X(I) = X(I) - R(I) 8210 CONTINUE * IF(XN1.GT.0.0) RN1 = RN1/XN1 IF((RN1 + RN1).LE.QUOT) GO TO 8220 IER = -IER WRITE(IOUT1,908) 908 FORMAT(2X,'*************** WARNING *******************',/ 1,5X,'THE SYSTEM IS ILL-CONDITIONNED') RETURN 8220 CONTINUE QUOT = RN1 IF(1.0 + SNGL(RN1).NE.1.0) GO TO 8200 RETURN ENTRY DUMP WRITE(IOUT1,923) WRITE(IOUT1,974) N DO 507 I = 1, N WRITE(IOUT1,973) (A(I,J), J = 1, N), F(I) 507 CONTINUE RETURN ENTRY DUMP2 WRITE(IOUT1,910) IER WRITE(IOUT1,922) RETURN 910 FORMAT(5X,'NUMBER OF ITERATIVE REFINEMENTS : ',I4) 922 FORMAT(//) 923 FORMAT (///) 973 FORMAT (2H , 7D11.4) 974 FORMAT (2H , 'MATRIX OF A(I,I), F(I) I=',I2//) END SUBROUTINE SBS(B) IMPLICIT DOUBLE PRECISION (A-G,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DIMENSION B(BG4) COMMON/SLOUGH/ A(BG4,BG4),F(BG4),T(BG4,BG4),IV(BG4),N,ITREF * DO 8110 K=1,N IF(K.GE.N) GO TO 8110 L = IV(K) TMP = B(L) B(L) = B(K) B(K) = TMP KP1 = K + 1 DO 8100 I=KP1,N B(I) = B(I) + T(I,K)*TMP 8100 CONTINUE 8110 CONTINUE * K = N 8200 B(K) = B(K)/T(K,K) IF(K.LE.1) RETURN TMP = -B(K) KP1 = K K = K-1 DO 8210 I=1,K B(I) = B(I) + T(I,KP1) *TMP 8210 CONTINUE GO TO 8200 END SUBROUTINE SAXMB(X,B,R) IMPLICIT DOUBLE PRECISION (A-G,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) DIMENSION X(BG4), B(BG4), R(BG4) COMMON/SLOUGH/ A(BG4,BG4),F(BG4),T(BG4,BG4),IV(BG4),N,ITREF * DO 8110 I=1,N SUM = -B(I) DO 8100 J = 1,N SUM = SUM + A(I,J)*X(J) 8100 CONTINUE R(I) = SUM 8110 CONTINUE RETURN END SUBROUTINE PRESS(P,T) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER BG1,BG2,BG3,BG4 PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) * * THIS SUBROUTINE USES THE POINTS OF A GIVEN P-T CURVE, READ * IN INP3, AND INTERPOLATES LINEARLY BETWEEN THE TWO CLOSEST * GIVEN POINTS TO RETURN THE PRESSURE. N.S. ADDED OCT 83 * COMMON/PP/TS(200),PS(200),NPP,INCP COMMON/IO/ INP1,INP2,INP3,INP4,IOUT1,IOUT2,IPUN,ISTEP,IPUNCH * * IF(T.GT.TS(NPP)) RETURN * GET THE APPROXIMATE POSITION OF POINT IN CURVE ARRAY ITEM=IDINT((T-TS(1))*NPP/(TS(NPP)-TS(1)))+1 I=ITEM IF(I.GT.NPP) I=NPP IF(TS(I)-T) 40,30,50 30 P=PS(I) RETURN 40 K=I+1 IF(TS(K).GE.T) GO TO 60 I=K GO TO 40 50 K=I-1 IF(TS(K).LE.T) GO TO 60 I=K GO TO 50 60 CONTINUE P=((T-TS(I))*(PS(K)-PS(I))/(TS(K)-TS(I))) + PS(I) RETURN END * SUBROUTINE PHI(PF) * *****THIS SUBROUTINE IS NOT IN USE AT PRESENT.. IN ORDER TO USE * YOU MUST FIRST ADJUST ALL COMMON BLOCKS AND ARRAY SIZES TO CONFORM * WITH PRESENT GASWORKS VERSION * ** IMPLICIT DOUBLE PRECISION (A-H,O-Z) ** INTEGER BG1,BG2,BG3,BG4 ** PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) ** CHARACTER*8 MIN ** COMMON/MORRIS/ MIN(BG1),COEFM(BG1,8),MINTRY(BG1),XMOLW ** 1 (BG1), TEMPK,NMT,MINCT,SPECM(BG1,8),ITM(BG1),SEQ(BG1),QES(BG1) ** COMMON/MORT/MTRY(310),GAMMA(310),WATM,FO2,FO2LOG,KK ** 1,OFLAG ** COMMON/FUGA/GASLOG(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3, ** 1 NGAS,IMIX,IDEAL,ITGAS ** COMMON/MORSEL/ GAMIN(BG1),B(BG1),SUMSOL(15),SOLTRY(15), ** 1 DGDN(BG3,20),SOLS(15,20),ISOLS(BG1),LIMSOL,JSOLS(BG1),IST(15) ** COMMON/VAL/BVAL(20), CVAL(20), BMIX, CMIX ** INTEGER SOLS, SPECM, SEQ, QES, OFLAG ** DOUBLE PRECISION MINTRY ** FUNC1(A1,A2,A3) = W1*A1 + W2*A2 + W3*A3 ** FUNC2(A1,A2,A3,A4,A5,A6) = W1*W1*A1 + W2*W2*A2 + W3*W3*A3 ** & + 2.*( W1*W2*A4 + W1*W3*A5 + W2*W3*A6 ) ** ** THIS SUBROUTINE COMPUTES FUGACITY COEFFICIENTS FOR GASES. ** WE ASSUME NON IDEAL MIXING FOR MIXTURES OF H2O, CO2 AND CH4 ** AND IDEAL MIXING FOR THE OTHER GASES. ** SINCE H2O, CO2 AND CH4 GENERALLY ACCOUNT FO 99 PERCENT OF THE ** THE GAS PHASE, THE RESULTING ERROR WHEN COMPUTING GAS COMPOSITIONS ** IS ACCEPTABLE. THE FUGACITY COEFFICIENTS, PHI, FOR H2O, CO2 AND CH ** ARE THEREFORE DEPENDENT ON T, P AND MOLE FRACTION Y OF THESE ** THREE COMPONENTS. PHI FOR OTHER GASES ARE ONLY DEPENDENT ON T AND ** P (PHI OF PURE GAS = PHI OF GAS IN MIXTURE) ** ENTRIES FOR H2O, CO2 AND CH4 GAS HAVE TO BE AT THE TOP OF THE GAS ** LIST IN DATA FILE SOLTHERM, H2O FIRST, CO2 SECOND AND CH4 THIRD. ** FUGACITY COEFFICIENTS FOR H2O, CO2 AND CH4 ARE CALCULATED WITH ** THEIR MOLE FRACTIONS NORMALIZED TO 1 (W1 W2 W3). ** IMIX IS A FLAG WHICH INDICATES WHICH OF THESE THREE GASES ARE ** PRESENT (H2O IS ALWAYS PRESENT) : ** IMIX=3 H2O-CO2-CH4 NON IDEAL MIXING ** IMIX=2 H2O-CO2 (NO CH4) NON IDEAL MIXING ** IMIX=1 H2O ONLY (NO CH4,CO2) IDEAL MIXING WITH OTHER GASES ** ** PHI ARE CALCULATED FROM COMPRESSIBILITY DATA, USING A VIRIAL ** EQUATION IN PRESSURE Z = P*V/R/T = 1 + VB*PF + VC*PF**2. ** FOR NON IDEAL MIXING WE USE Z = ZMIX. FOR A MULTICOMPONENT ** MIXTURE ZMIX = 1 + BMIX*P + CMIX*P*P WHERE BMIX AND CMIX ARE ** EXPRESSED IN TERMS OF CROSS VIRIAL COEFFICIENTS BIJ' AND CIJK'. ** VALUES FOR BIJ'S AND CIJK'S WHERE OBTAINED FROM LEAST SQUARE ** FIT OF PUBLISHED PVT OR SOLUBILITY DATA. PHI(I)'S ARE THEN ** OBTAINED BY TAKING THE APPROPRIATE PARTIAL DERIVATIVES. ** SEE SPYCHER AND REED, IN PREPARATION. ** CROSS VIRIAL COEFFICIENTS ARE ASSIGNED AS FOLLOW (1=H2O, 2=CO2, ** 3=CH4): ** CVB(1) : B12 ** CVB(2) : B13 ** CVB(3) : B23 ** CVC(1) : C112 ** CVC(2) : C122 ** CVC(3) : C113 ** CVC(4) : C133 ** CVC(5) : C223 ** CVC(6) : C233 ** CVC(7) : C123 ** ********* THIS ROUTINE NEEDS SUBROUTINE MIXUP FOR NON IDEAL MIXING ** ********* AS OF FEB 86, BECAUSE OF LACK OF EXPERIMENTAL DATA TO ** TO COMPUTE THESE VALUES, C113, C133 AND C123 ARE SET TO ** ZERO. ** DATA FOR H2O-CH4 IS EXTRAPOLATED ABOVE 100 DEG.C. DATA FOR ** CO2-CH4 IS EXTRAPOLATED ABOVE 240 DEG.C. ********* WATCH ! THE CURRENT PRESSURE LIMIT FOR COMPUTING ** FUGACITY COEFFICIENTS WITH ACCURACY AROUND 1-2% OR LESS ** IS ABOUT 500 ATM FROM 150 TO 300 DEG.C, AND 200 ATM FROM ** FROM 80 TO 150 DEG.C. IN THIS LOWER T RANGE, ERRORS CAN BE ** UP TO 10% AT 500 BARS. NOT RECOMMENDED BELOW ** 75 DEG.C. EXCEPT AT LOW P. H2 DATA IS VALID UP TO ** 1300 BARS FROM 0 TO 600 DEG.C. A SET OF COEFFS VALID ABOVE ** 350 C UP TO 1000 BARS IS AVAILBALE FOR H2O-CO2. SEE SPYCHER ** AND REED, IN PREPARARION. ** ** SEE SOLTHERM DATA FILE FOR SOURCE OF CURRENT DATA FOR PURE AND ** MIXED GASES. ** ** WRITTEN BY N.SPYCHER FEB 86 ** **** LAST REVISED JAN 87. ADD DGDN TO EASE CONVERGENCE IN SOME ** RARE CASES. N.S. ** ** N1=1 ** IF(IDEAL.NE.0.OR.IMIX.EQ.0) GO TO 10 ** CALL MIXUP(VB,VC,CVB,CVC) ** ********** ** N1 = IMIX + 1 ** ** IF IMIX IS 2, BVAL(3) AND CVAL(3) ARE OVERWRITTEN BELOW ** 10 CONTINUE ** IDEAL MIXING ** DO 15 I = N1,NGAS ** BVAL(I) = VB(I) ** 15 CVAL(I) = VC(I) ** ** NOW WE COMPUTE FUGACITY COEFFICIENTS ** DO 20 I=1,NGAS ** VV = BVAL(I)*PF + CVAL(I)*PF*PF/2. ** SKIP IF VV IS TOO SMALL OR TOO LARGE AND KEEP OLD GAMIN VALUE ** IF(DABS(VV).GT. 10.) GO TO 20 ** GAMIN(I) = DEXP(VV) ** B(I) = 1. ** 20 CONTINUE ** ** C-OUT BELOW AND FUNCTIONS ABOVE IF CONVERGENCE PROBLEMS ARISE, ** PARTICULARTLY AT HIGH CO2 MOLE FRAC. IN GAS PHASE ** ** NOW WE COMPUTE DERIVATIVES OF GAMIN W/ RESPECT TO Y (NON IDEAL ** MIXING) TO EASE CONVERGENCE IN MINERA IN SOME RARE CASES. ** SKIP IF: GASES ARE NOT SATURATED, OR IDEAL MIXING. ** IF(ISOLS(1).EQ.0) RETURN ** IF(IDEAL.NE.0.OR.IMIX.EQ.0) RETURN ** SUMY = MINTRY(1) + MINTRY(2) ** IF(IMIX.GT.2) SUMY = SUMY + MINTRY(3) ** D1 = FUNC1(VB(1),CVB(1),CVB(2)) ** D2 = FUNC1(CVB(1),VB(2),CVB(3)) ** D3 = FUNC1(CVB(2),CVB(3),VB(3)) ** D12 = ( 2.*CVB(1) + 2.*BMIX - 2.*D1 - 2.*D2) * PF ** D13 = ( 2.*CVB(2) + 2.*BMIX - 2.*D3 - 2.*D1) * PF ** D23 = ( 2.*CVB(2) + 2.*BMIX - 2.*D1 - 2.*D3) * PF ** E1 = FUNC2(VC(1),CVC(2),CVC(4),CVC(1),CVC(3),CVC(7)) ** E2 = FUNC2(CVC(1),VC(2),CVC(6),CVC(2),CVC(7),CVC(5)) ** E3 = FUNC2(CVC(3),CVC(5),VC(3),CVC(7),CVC(4),CVC(6)) ** E12 = (6.*FUNC1(CVC(1),CVC(2),CVC(7)) + 6.*CMIX - 6.*E1 ** & - 6.*E2) * PF*PF/2. ** E13 = (6.*FUNC1(CVC(3),CVC(7),CVC(4)) + 6.*CMIX - 6.*E1 ** & - 6.*E3) * PF*PF/2. ** E23 = (6.*FUNC1(CVC(7),CVC(5),CVC(6)) + 6.*CMIX - 6.*E2 ** & - 6.*E3) * PF*PF/2. ** ** DGDN(1,1) = GAMIN(1)/SUMY * (( 2.*VB(1) + 2.*BMIX - 4.*D1) * PF ** & + (6.*FUNC1(VC(1),CVC(1),CVC(3)) + 6.*CMIX - 12.*E1)*PF*PF/2.) ** DGDN(1,2) = GAMIN(1)/SUMY * (D12 + E12) ** DGDN(1,3) = GAMIN(1)/SUMY * (D13 + E13) ** ** DGDN(2,1) = GAMIN(2)/SUMY * (D12 + E12) ** DGDN(2,2) = GAMIN(2)/SUMY * (( 2.*VB(2) + 2.*BMIX - 4.*D2) * PF ** & + (6.*FUNC1(CVC(2),VC(2),CVC(5)) + 6.*CMIX - 12.*E2)*PF*PF/2.) ** DGDN(2,3) = GAMIN(2)/SUMY * (D23 + E23) ** ** IF(IMIX.LT.3) RETURN ** DGDN(3,1) = GAMIN(3)/SUMY * (D13 + E13) ** DGDN(3,2) = GAMIN(3)/SUMY * (D23 + E23) ** DGDN(3,3) = GAMIN(3)/SUMY * (( 2.*VB(3) + 2.*BMIX - 4.*D3) * PF ** & + (6.*FUNC1(CVC(4),CVC(6),VC(3)) + 6.*CMIX - 12.*E3)*PF*PF/2.) ** RETURN ** END ** SUBROUTINE MIXUP(B,C,CB,CC) ** IMPLICIT DOUBLE PRECISION (A-H,O-Z) ** INTEGER BG1,BG2,BG3,BG4 ** PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) ** COMMON/FUGA/GASLOG(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3, ** 1 NGAS,IMIX,IDEAL,ITGAS ** COMMON/VAL/BVAL(20),CVAL(20),BMIX,CMIX ** DIMENSION B(20),C(20),CB(3),CC(7) ** ** THIS ROUTINE COMPUTES BMIX, CMIX AND THEIR DERIVATIVE WITH ** RESPECT TO WI AT T,P WJ NOT WI CONSTANT. THE RESULTING VALUES ** BVAL AND CVAL ARE THEN USED IN PHI TO CALCULATE FUGACITY COEFFS. ** SEE SUBROUTINE PHI FOR COMMENTS. FOR ENTHALPIES, WE PLUG IN ** VALUES OF DB/DT, DC/DT, ETC. AND THEREFORE WE COMPUTES ** DBMIX/DT, DCMIX/DT ETC. WE THEN USE RESULTING VALUES OF BVAL ** AND CVAL TO COMPUTE ENTHALPIES OF GASES. SEE SUBROUTINE ENTHAL ** B, C ARE 2ND AND 3RD VIRIAL COEFFS; CB, CC ARE 2ND AND 3RD CROSS ** VIRIAL COEFFS AT GIVEN TEMP. (IF MIXUP IS CALLED FROM ENTHALP ** THEN B ETC.. ARE -T**2 * DB/DT ETC...) ** CROSS VIRIAL COEFFICIENTS ARE ASSIGNED AS FOLLOW AS IN SUBROUTINE ** PHI. ** ** WRITTEN BY NICOLAS SPYCHER, U.OF.O JAN 1986 ** LAST REVISED: FEB 86 BY: N.S ** ** FUNC1(A1,A2,A3) = W1*A1 + W2*A2 + W3*A3 ** FUNC2(A1,A2,A3,A4,A5,A6) = W1*W1*A1 + W2*W2*A2 + W3*W3*A3 ** & + 2.*( W1*W2*A4 + W1*W3*A5 + W2*W3*A6 ) ** ** WE HAVE TO NORMALIZE THE SUM OF THE MOLE FRACTIONS OF H2O,CO2,CH4 ** TO 1 BECAUSE WE IGNORE THE MIXING EFFECTS OF ALL OTHER GASES. ** ** WTOT = Y(1) + Y(2) ** W3=0. ** IF(IMIX.LT.3) GO TO 5 ** WTOT=WTOT+Y(3) ** W3=Y(3)/WTOT ** 5 W1=Y(1)/WTOT ** W2=Y(2)/WTOT ** IF(Y(1).LE.1.D-23) W1 = 0. ** IF(Y(2).LE.1.D-23) W2 = 0. ** IF(Y(3).LE.1.D-23) W3 = 0. ** NON IDEAL MIXING ** ** BMIX = FUNC2(B(1),B(2),B(3),CB(1),CB(2),CB(3)) ** CMIX = W1*W1*W1*C(1) + W2*W2*W2*C(2) + W3*W3*W3*C(3) ** & + 3.*( W1*W1*W2*CC(1) ** & + W1*W2*W2*CC(2) + W1*W1*W3*CC(3) + W1*W3*W3*CC(4) ** & + W2*W2*W3*CC(5) + W2*W3*W3*CC(6) + 2.*W1*W2*W3*CC(7) ) ** BVAL AND CVAL ALSO NEEDED IN PSAT TO SET UP PARTIAL DERIVATIVES ** OF GAMIN WITH RESPECT TO PRESSURE (TO SOLVE FOR SAT. PRESSURE) ** ** Z = 1 + BMIX*P + CMIX*P*P ** BVAL(1) = 2.*FUNC1(B(1),CB(1),CB(2)) - BMIX ** CVAL(1) = 3.*FUNC2(C(1),CC(2),CC(4),CC(1),CC(3),CC(7) ) ** & - 2.*CMIX ** ** BVAL(2) = 2.*FUNC1(CB(1),B(2),CB(3)) - BMIX ** CVAL(2) = 3.*FUNC2(CC(1),C(2),CC(6),CC(2),CC(7),CC(5) ) ** & - 2.*CMIX ** ** BVAL(3) = 2.*FUNC1(CB(2),CB(3),B(3)) - BMIX ** CVAL(3) = 3.*FUNC2(CC(3),CC(5),C(3),CC(7),CC(4),CC(6) ) ** & - 2.*CMIX ** IF IMIX IS 2, BVAL(3) AND CVAL(3) WILL BE OVERWRITTEN ** RETURN ** END ** SUBROUTINE ENTHAL(TC) ** ** IMPLICIT DOUBLE PRECISION (A-H,O-Z) ** INTEGER BG1,BG2,BG3,BG4 ** PARAMETER(BG1=450,BG2=650,BG3=42,BG4=85) ** COMMON/FROID/ SINC,BIT,SLIM,AKLOG,APFCL(BG2,5), ** 1 BPFCL(BG1,5),XPFC(11,5), XF(11), COE(5), ALG(4), ** 2 XDH(4),ALK(5),LGK(BG2),TEMP,TEMPC, HEATC,COMTOT(BG3), HEATM, ** 3 TOTMIX,CFH(5),CFT(5),VBCOE(20,3),VCCOE(20,3),CVBCF(3,3), ** 4 CVCCF(7,3), IFRAC ** COMMON/FUGA/GASLOG(20),Y(20),VB(20),VC(20),CVB(3),CVC(7),W1,W2,W3, ** 1 NGAS,IMIX,IDEAL,ITGAS ** COMMON/VAL/BVAL(20),CVAL(20) ** DIMENSION WATCF1(4),WATCF2(5),DBDT(20),DCDT(20),DCBDT(3),DCCDT(7) ** DOUBLE PRECISION LGK ** DATA WATCF1/2.335628,-.00441813,2.30965D-5,-777.40446/, ** 1 WATCF2/240.912149,-.35584278,.00021415,-7.195396D+4,7.895231D-6/ ** DATA A1,B1,C1/7.17,2.56D-3,.08D+5/, ** 1 A2,B2,C2/10.55,2.16D-3,-2.04D+5/, ** 2 A3,B3,C3,D3/2.975,18.329D-3,.346D+5,-4.303D-6/ ** FUNCTION DHT FROM CP = A + B*T + C/T**2 ** DHT(A,B,C,T) = A*(T-298.15) + B/2.*(T*T-298.15*298.15) - ** 1 C*(1./T - 1./298.15) ** ** THIS SUBROUTINE CALCULATES ENTHALPIES OF WATER, H2O GAS, ** CO2 AND CH4. WATER ENTHALPY IS ASSUMED TO BE PRESSURE ** INDEPENDENT. ** THE ENTHALPY OF WATER FOR A GIVEN TEMPERATURE IS ** CALCULATED AT THE CORRESPONDING SATURATION PRESSURE. SINCE ** WATER IS NEARLY INCOMPRESSIBLE, THE DIFFERENCE IN ENTHALPY ** AT PRESSURES HIGHER THAN THE SATURATION PRESSURE IS NEGLIGIBLE. ** THE ENTHALPY IS CALCULATED FROM A POWER FUNCTION OF THE ** TEMPERATURE. ** H(TK) = A + B*TK + C*TK**2 + D/TK + E/TK**2 ** WHERE TK IS THE TEMPERATURE IN DEG. KELVIN. THE COEFFICIENTS ** A,B,C,D,E WERE LEAST-SQUARE FITTED FROM HELGESON AND KIRKHAM ** 1974 PART I (AJS, 274, TABLE 37,P.1184). THE VALUES OF THESE ** COEFFICIENTS ARE SUCH THAT THE CALCULATED ENTHALPY H(TK) IS ** IN KCAL/MOLE AND IS ZERO AT 0 DEG.C. ** ENTHALPY OF H2O, CO2 AND CH4 GAS: ** THE ENTHALPIES OF THESE GASES ARE CALCULATED FROM HEAT ** CAPACITY DATA (BARIN&KNACKE,1973,1977,THERMOCHEMICAL PROPERTIES ** OF INORGANIC SUBSTANCES, SPRINGER-VERLAG) AND FROM VIRIAL ** COEFFICIENTS (SAME DATA THAT WE USE TO COMPUTE FUGACITY ** COEFFICIENTS, SEE SUBROUTINE PHI). THE ENTHALPY IS THEN: ** H GAS MIXTURE = SUM OF (MOLES(I) * (DH/DT (I) + DH/DP (I))) ** IF IDEAL MIXING, THEN DH/DP IS FOR PURE GAS. IF NO-IDEAL MIXING ** THEN DH/DP IS FOR THE GAS IN THE MIXTURE. ** DHT IS THE CHANGE OF ENTHALPY WITH TEMPERATURE, CALCULATED FROM ** A MAIER-KELLY HEAT CAPACITY POWER FUNCTION, AND DHP IS THE ** CHANGE WITH PRESSURE, CALCULATED FROM THE FOLLOWING EQUATION: ** DHP = R*TK**2 * ( DB/DT*P + DC/DT*P**2/2 ) ** WHERE B AND C ARE THE SECOND AND THIRD VIRIAL COEFFICIENTS ** OF THE VIRIAL EQ. OF STATE IN PRESSURE. FOR NON IDEAL MIXTURES ** WE TAKE THE APPROPRIATE PARTIAL DERIVATIVES. ** ****** WARNING: DHP MAY BE UP TO 8% IN ERROR AT HIGHER T,P ** HOWEVER THE ERROR IN THE TOTAL ENTHALPY OF THE GASES DOES NOT ** EXEED A FEW PRECENT. MAKE SURE TO RESPECT T,P LIMITS (SEE ** SUBROUTINE PHI). ** ENTHAPIES OF CO2 AND CH4 ARE ZERO AT 0 DEG.C. AND 0 BAR. ** FOR H2O GAS, AT 0 DEG.C. 0.006 BAR, H = ENTH. OF VAPORIZATION. ** (10.77 KCAL/MOLE) ** WRITTEN BY NICOLAS SPYCHER, U.OF.O FEB 86 ** LAST REVISED: BY: ** ** TK=TC + 273.15 ** GETS HEAT OF WATER IF LESS THAN 100 DEG.C. COEFFICIENTS ** WATCF1 VALID ONLY FROM 0 TO 99 DEG.C. ** IF(TEMP.GE.100.) GO TO 20 ** ENTHW = WATCF1(1) + WATCF1(2)*TK + WATCF1(3)*TK*TK + WATCF1(4)/TK ** GETS HEAT OF WATER IF 100 DEG.C.OR ABOVE. WATCF2 COEFFICIENTS ** VALID FROM 100 TO 350 DEG.C. ** 20 CONTINUE ** ENTHW = WATCF2(1) + WATCF2(2)*TK + WATCF2(3)*TK*TK + WATCF2(4)/TK ** & + WATCF2(5)/TK/TK ** RETURN ** ** ENTRY PFENTH(TC,P) ** *************** ** TK=TC + 273.15 ** GET HEAT OF GASES. 1 H2O; 2 CO2; 3 CH4; ** ENTHM(1)= DHT(A1,B1,C1,TK)/1000. + 10.970 ** ENTHM(2)= DHT(A2,B2,C2,TK)/1000. + .210 ** ENTHM(3)= DHT(A3,B3,C3,TK)/1000. + D3*(TK**3.-298.15**3.)/3000. ** 1 + .210 ** ** IF(IDEAL-1) 53,52,51 ** 51 RETURN ** PREF IS ONE ATM SINCE CP DATA IS AT ONE ATM ** IDEAL MIXING/REAL GAS WE GET DHP FOR PURE GASES THEN ADD TO DHT ** 52 CONTINUE ** DO 50 I=1,IMIX ** DB = 2.*VBCOE(I,1)/TK + VBCOE(I,2) ** DC = 2.*VCCOE(I,1)/TK + VCCOE(I,2) ** DHDP(I) = 1.987 * (DB + DC*P)/1000. ** DHP = 1.987 * (DB*(P-1.) + DC*(P*P-1.)/2.)/1000. ** 50 ENTHM(I) = ENTHM(I) + DHP ** RETURN ** ** 53 CONTINUE ** IF(Y(1).EQ.0.) GO TO 52 ** NON IDEAL MIXING/REAL GAS. WE CALCULATE DHP FOR GAS IN MIXTURE ** THEN ADD TO DHDT. PREF IS ONE ATM SINCE CP DATA IS AT 1 ATM. ** BELOW WE GET VALUES OF T**2 * DB/DT (ALSO FOR C AND CROSS B,C) ** DO 60 I = 1,IMIX ** DBDT(I) = 2.*VBCOE(I,1)/TK + VBCOE(I,2) ** 60 DCDT(I) = 2.*VCCOE(I,1)/TK + VCCOE(I,2) ** DO 62 I = 1,3 ** 62 DCBDT(I) = 2.*CVBCF(I,1)/TK + CVBCF(I,2) ** DO 64 I = 1,7 ** 64 DCCDT(I) = 2.*CVCCF(I,1)/TK + CVCCF(I,2) ** ** CALL MIXUP(DBDT,DCDT,DCBDT,DCCDT) ** ********** ** DO 66 I=1,IMIX ** DHDP(I) = 1.984* (BVAL(I) + CVAL(I)*P)/1000. ** DHP= 1.984* (BVAL(I)*(P-1.) + CVAL(I)*(P*P-1.)/2.)/1000. ** 66 ENTHM(I) = ENTHM(I) + DHP ** RETURN ** END