ht1.inc

# Used by ht1 and grid calculations to calculate a predominance diagram based on the 'predominance' criterion.
PRINT
   -reset false
PHASES
Fix_H+
   H+ = H+
   log_k 0.0
SELECTED_OUTPUT
   -high_precision true
   -reset false
USER_PUNCH
-headings
-start
20 nout1 = 0
40 nout3 = 0
50 main$ = "<mainspecies>"
60 IF STEP_NO = 0 THEN 370 ELSE 70
70 totel = SYS(main$, n, n$, t$, c)
80 h2o = TOT("water")

90 REM Filter out gases
100 FOR i = 1 TO n
110   IF(c(i) <= 0) THEN 190
120   j = INSTR(n$(i), "(g)")
130   IF(j > 0) THEN 180
140   PUNCH n$(i), c(i)/h2o
150 REM   PRINT i, n$(i), c(i)/h2o
160   nout1 = nout1+1
170   IF(nout1 >= 3) THEN 190
180 NEXT i

190 REM no forced mineral species (see ht1stability.inc)
191 nout2 = 0

200 REM add these constraints
201 IF(SI("H2(g)") > 0) THEN 210 ELSE 230
210 PUNCH "H2(g) > 1 atm", 10^SI("H2(g)")
220 nout3 = nout3+1
230 IF(SI("O2(g)") > -0.677) THEN 240 ELSE 260
240 PUNCH "O2(g) > 0.21 atm", 10^SI("O2(g)")
250 nout3 = nout3+1
260 IF(SI("CH4(g)") > 0) THEN 270 ELSE 290
270 PUNCH "CH4(g) > 1 atm", 10^SI("CH4(g)")
280 nout3 = nout3+1

290 REM No 'carry' variables for viewing
291 nout4 = 0

300 REM the 5 system variables (always required)
301 po2 = SI("O2(g)")
310 ph2 = SI("H2(g)")
320 ph = -LA("H+")
330 pe = -LA("e-")
340 PUNCH "pH", ph, "pe", pe, "PO2", po2, "PH2", ph2, "temp", TC

350 REM output the 5 counters
351 nout5 = 5
360 PUNCH nout1, nout2, nout3, nout4, nout5
370 END
-end