minstab1.inc

# Used by ht1/grid calculations to calculate a mineral stability diagram
SELECTED_OUTPUT
  reset false
  high_precision true
USER_PUNCH
-start
10 acc = 1e-10
20 nout1 = 0
30 nout2 = 0
40 nout3 = 0
50 nout4 = 0
60 nout5 = 0
70 IF(STEP_NO = 0) THEN 670

80 t = SYS("phases", n, n$, t$, six)
90 PRINT "No. of potential phases=", n, " Max. SI=", t
100 PRINT "Phases saturated or supersatd:"
110 DIM mineral$(n), cmin(n)
120 IF(t < -acc) THEN 270
130 j = 0
140 k = 0
150 FOR i = 1 TO n
160   IF(INSTR(n$(i), "Fix")) > 0 THEN 250
170   found$ = ""
180   IF ABS(six(i)) <= acc THEN found$ = "Stable mineral"
190   IF(six(i) < -acc) THEN 270
200   PRINT i, n$(i), t$(i), six(i), found$
210   IF(ABS(six(i)) > acc) THEN 250
220   k = k+1
230   cmin(k) = EQUI(n$(i))
240   mineral$(k) = n$(i)
250 NEXT

260 REM
270 IF(k > 0) THEN 400
271 k = 1
272 mineral$(1) = "No minerals present"
273 nout1 = 1
274 GOTO 400

275 REM use this block if you want to show most abundant aq species
280 totaq = SYS("aq", naq, naq$, taq$, aq)
290 REM for i=1 to naq
300 REM   print i,naq$(i),taq$(i),aq(i)
310 REM next
320 mineral$(1) = naq$(1)
330 cmin(1) = aq(1)
340 k = 1
350 IF naq = 1 THEN 400
360 mineral$(2) = naq$(2)
370 cmin(2) = aq(2)
380 k = 2
390 GOTO 430

400 n = k
410 REM Sort
420 GOSUB 750
430 PRINT "Species sorted by mol present, largest first:"
440 FOR i = k TO 1 STEP -1
450   PRINT i, mineral$(i), cmin(i)
460 NEXT i
470 PUNCH mineral$(k), cmin(k)
480 nout1 = 1
490 IF k < 2 THEN 560
500 PUNCH mineral$(k-1), cmin(k-1)
510 nout1 = 2
520 IF k < 3 THEN 560
530 PUNCH mineral$(k-2), cmin(k-2)
540 nout1 = 3

560 REM now add usual system variables
570 IF(SI("H2(g)") > 0) THEN 580 ELSE 600
580 PUNCH"H2(g) > 1 atm", 10^SI("H2(g)")
590 nout3 = nout3+1
600 IF(SI("O2(g)") > -0.677) THEN 610 ELSE 630
610 PUNCH"O2(g) > 0.21 atm", 10^SI("O2(g)")
620 nout3 = nout3+1
630 IF(SI("CH4(g)") > 0) THEN 640 ELSE 660
640 PUNCH"CH4(g) > 1 atm", 10^SI("CH4(g)")
650 nout3 = nout3+1

660 REM
670 po2 = SI("O2(g)")
680 ph2 = SI("H2(g)")
690 ph = -LA("H+")
700 pe = -LA("e-")
710 PUNCH"pH", ph, "pe", pe, "PO2", po2, "PH2", ph2, "temp", tc
720 nout5 = 5
730 PUNCH nout1, nout2, nout3, nout4, nout5
740 END

750 REM insertion sort (from Numerical Recipes) - returned in ascending order
760 FOR j = 2 TO n
770   c = cmin(j)
780   tm$ = mineral$(j)
790   FOR i = j-1 TO 1 STEP -1
800     IF(cmin(i) <= c) THEN GOTO 850
810     mineral$(i+1) = mineral$(i)
820     cmin(i+1) = cmin(i)
830   NEXT i
840   i = 0
850   cmin(i+1) = c
860   mineral$(i+1) = tm$
870 NEXT j
880 RETURN
-end