risk.inc

# calculates a simplified predominance diagram with fields classified into one of:
#  'Mobile' (= dissolved), 'Immobile' (= solid) or "May be mobile" (= adsorbed)
PRINT
   -reset false
PHASES
Fix_H+
   H+ = H+
   log_k 0.0
SELECTED_OUTPUT
    -high_precision true
   -reset false
USER_PUNCH
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

#130 REM Print a list of minerals to the output file for pasting in
#140 REM if required
#150 REM This will require the PRINT -reset false statement to be
#160 REM commented out or reset
#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

260 totel = sys(main$, n, n$, t$, c)
280 FOR i = 1 TO n
290 IF(c(i) <= 0) THEN 670
300 REM Filter out gases
310 j = instr(n$(i), "(g)")
320 IF(j > 0) THEN 360
321 if (t$(i)="aq") then 322 else 324
322 t$(i)="Mobile"
323 goto 330
324 if (t$(i)="equi") then 325 else 327
325 t$(i)="Immobile"
326 goto 330
327 if (t$(i)="surf") then 328 else 330
328 t$(i)="May be mobile"
330 PUNCH t$(i), c(i)
340 nout1 = nout1+1
350 IF(nout1 >= 3) THEN 670
360 NEXT i

670 IF(SI("H2(g)") > 0) THEN 680 ELSE 700
680 PUNCH"Hydrogen", SI("H2(g)")
690 nout3 = nout3+1
700 IF(SI("O2(g)") > -0.677) THEN 710 ELSE 730
710 PUNCH"Oxygen", SI("O2(g)")
720 nout3 = nout3+1
730 IF(SI("CH4(g)") > 0) THEN 740 ELSE 790
740 PUNCH"Methane", SI("CH4(g)")
750 nout3 = nout3+1

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