ht1s.inc

# ht1 but adds "(s)" to all solid phases (to avoid confusion with solution species), eg Fe(OH)(s)
PRINT
   reset false
PHASES
Fix_H+
   H+ = H+
   log_k 0.0
SELECTED_OUTPUT
   -high_precision true
   -reset false
USER_PUNCH
-headings
10 ns = 0
20 nout1 = 0
30 nout2 = 0
40 nout3 = 0
50 nout4 = 0
51 nout5 = 0

80 main$ = "<mainspecies>"
100 DIM s$(100)
110 IF STEP_NO = 0 THEN 880 ELSE 260

260 totel = sys(main$, n, n$, t$, c)
265 h2o = TOT("water")
270 REM Filter out gases
280 FOR i = 1 TO n
290 IF(c(i) <= 0) THEN 670
310 j = instr(n$(i), "(g)")
311 simin(i)=-99
320 IF(j > 0) THEN 360
324 nt$=n$(i)
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 nt$,c(i)/h2o
335 if (t$(i)="equi") then simin(i)=SI(n$(i))
340 nout1 = nout1+1
350 IF(nout1 >= 3) THEN 670
360 NEXT i

670 rem no type2's

671 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)")
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


# 4. carry
# 2. top 3 minerals SI's
760 for i=1 to nout1
761 punch "SI_"+chr$(48+i),simin(i)
762 nout4=nout4+1
763 next i


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

870 PUNCH nout1, nout2, nout3, nout4, nout5
880 END
-end
#include printphases.inc  # PS debug 2 resolution 1