|
|
# Used by 'hunt and track' method to calculate a 'predominance' diagram based on the 'stability' criterion
PRINT
#### -reset false
PHASES
Fix_H+
H+ = H+
log_k 0.0
SELECTED_OUTPUT
-high_precision true
-reset false
USER_PUNCH
-headings
-start
10 ns = 0
20 h2o = TOT("water")
# NB if <type> is defined by numberOfCharacterTags, specify with nested quotes, eg '"Stability"'
#60 type$ = <type>
70 type$ = "Stability"
80 main$ = "<mainspecies>"
90 REM maximum number of surfaces=100
100 DIM s$(100)
110 IF STEP_NO = 0 THEN 880 ELSE 170
130 REM Print a list of minerals to the output file for pasting in
140 REM if required
150 REM This will require any PRINT -reset false statement to be
160 REM commented out or reset to true
170 totm = SYS("phases", nm, nm$, tm$, cm)
180 FOR i = 1 TO nm
190 k = instr(nm$(i), "(g)")
200 k = k+instr(nm$(i), "Fix")+instr(nm$(i), "FIX")+instr(nm$(i), "fix")
210 IF(k > 0) THEN 230
220 PRINT" ", CHR$(35)+pad(nm$(i), 30), "0", "0"
230 NEXT i
240 nout1=0
250 REM Find top 3 species omitting surface species
260 totel = SYS(main$, n, n$, t$, c)
270 REM Filter out surfaces and gases
280 FOR i = 1 TO n
290 IF(c(i) <= 0) THEN 380
300 IF(t$(i) = "surf") THEN 360
310 j = instr(n$(i), "(g)")
320 IF(j > 0) THEN 360
325 IF (t$(i)="equi" AND INSTR(n$(i),"(g)")<>0) THEN nt$=n$(i)+"(s)" #if aqueous and solid species have same name
330 PUNCH n$(i), c(i)/h2o
340 nout1 = nout1+1
#341 PRINT "Dominant species ",nout1," ",n$(i)," ",c(i)/h2o
350 IF(nout1 >= 3) THEN 380
360 NEXT i
380 REM Now add any SURFACES...
390 totsurf = SYS("surf", ne, ne$, te$, ce)
400 FOR i = 1 TO ne
410 IF(te$(i) <> "surf") THEN 540
420 k = instr(ne$(i), "_")
430 IF(k > 1) THEN ne$(i) = MID$(ne$(i), 1, k-1)
440 PRINT i, ne$(i), te$(i), ce(i)/h2o
450 IF ns <= 0 THEN 490
460 FOR ii = 1 TO ns
470 IF ne$(i) = s$(ii) THEN 540
480 NEXT ii
490 ns = ns+1
500 s$(ns) = ne$(i)
510 s = SURF(main$, s$(ns))
512 if (s<1e-14) then 540 # remove non-existent surface species
520 PUNCH s$(ns)+"-"+trim(main$), s/h2o
530 nout1 = nout1+1
540 NEXT i
560 REM Now MINERAL constraints...
570 IF(type$ <> "Stability") THEN 650
580 FOR i = 1 TO n
590 j = instr(n$(i), "(g)")
600 IF(j > 0) THEN 640
610 IF(t$(i) = "equi" AND c(i) > 0) THEN 620 ELSE 640
620 PUNCH n$(i), c(i)/h2o
630 nout2 = nout2+1
#631 PRINT "Mineral constraints ",nout2," ",n$(i)," ",c(i)/h2o
640 NEXT i
650 nout3=0
660 REM Next any OTHER constraints...
670 IF(SI("H2(g)") > 0) THEN 680 ELSE 700
680 PUNCH"H2(g) > 1 atm", 10^SI("H2(g)")
690 nout3 = nout3+1
700 IF(SI("O2(g)") > -0.677) THEN 710 ELSE 730
710 PUNCH"O2(g) > 0.21 atm", 10^SI("O2(g)")
711 print "SI O2=",SI("O2(g)")
720 nout3 = nout3+1
730 IF(SI("CH4(g)") > 0) THEN 740 ELSE 760
740 PUNCH"CH4(g) > 1 atm", 10^SI("CH4(g)")
750 nout3 = nout3+1
760 REM Add any 'carry' variables here...
761 nout4=0
770 REM Now add 'SYSTEM' variables ...
780 REM Always include these 4...
790 po2 = SI("O2(g)")
800 ph2 = SI("H2(g)")
810 ph = -LA("H+")
820 pe = -LA("e-")
830 PUNCH"pH", ph, "pe", pe, "PO2", po2, "PH2", ph2, "temp", tc
840 nout5 = 5
860 REM Finally send the counters...
870 PUNCH nout1, nout2, nout3, nout4, nout5
880 END
-end
|
|