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
|
|