Nucelation RATES M=0

[Isil]

Hello,
I am trying to implement a precipitation rate for aragonite, considering the nucleation. The system fails if M is 0 at the beginning. It seems that moles are calculated but not saved. Thanks to the print function, I saw that moles are calculated, but M is not saved, and no precipitation of Aragonite occurs.
Thank you for any help in advance

Code: [Select]

RATES
    Aragonite
-start
 10 si_ar = SI("Aragonite")
 20 sr_thr = SR("Aragonite")
 30 omega_ar = 10^(si_ar)
 40 moles = 0
 50 sak = 10^(1.05)
 60 threshold = 1e-10
 70 IF(M <= 0 AND sr_thr < 1) THEN GOTO 130
80 IF(M <= 0 AND sr_thr <= 1) THEN GOTO 150
90 IF(M >= threshold) THEN GOTO 130
100 r_ar=1
110 rate = 1e-10
120 GOTO 140
130 rate = -sak*M
140 moles = rate*((omega_ar-1)^(2.71))*TIME/3600/10^6
150 SAVE moles
160 PRINT M rate moles omega_ar
-end
     END


[dlparkhurst]

Here is a script you can work from. However, your rate equation is undefined for dissolution because you are raising a negative number to a fractional power. You will have failures with different time steps because the method may generate some dissolution in intermediate sub time steps.

Code: [Select]

RATES
    Aragonite
-start
 10 si_ar = SI("Aragonite")
 20 sr_thr = SR("Aragonite")
 30 omega_ar = 10^(si_ar)
 40 moles = 0
 50 sak = 10^(1.05)
 60 threshold = 1e-10
 70 IF(M <= 0 AND sr_thr < 1) THEN GOTO 130
80 IF(M <= 0 AND sr_thr <= 1) THEN GOTO 150
90 IF(M >= threshold) THEN GOTO 130
100 r_ar=1
110 rate = 1e-10
120 GOTO 140
130 rate = -sak*M
140 moles = rate*((omega_ar-1)^(2.71))*TIME/3600/10^6
150 SAVE moles
160 PRINT M rate moles omega_ar
-end
END
RATES
    Aragonite
-start
 10 si_ar = SI("Aragonite")
 20 sr_thr = SR("Aragonite")
 30 omega_ar = 10^(si_ar)
 40 moles = 0
 50 sak = 10^(1.05)
 60 threshold = 1e-10
 65 REM 0 moles and SR <= 0
 70 IF(M <= 0 AND sr_thr < 1) THEN GOTO 200
 75 REM 0 moles and supersaturated
 80 IF(M <= 0 AND sr_thr > 1) THEN GOTO 300
 85 REM moles > 0
 90 GOTO 400
200 REM 0 moles and SR <= 0
205 PRINT "0 moles and SR <= 0"
210 rate = 0
220 GOTO 500
300 REM 0 moles and supersaturated
305 PRINT "0 moles and supersaturated"
310 rate = 1e-10
320 GOTO 500
400 REM moles > 0
405 PRINT "moles > 0"
410 rate = sak*M
500 REM calculate moles
510 k = 1/3600/10^6
515 k = 1e-1
520 moles = -k*rate*((omega_ar - 1)^(2.71))*TIME
530 SAVE moles
540 PRINT M, rate, moles, omega_ar
-end
END
SOLUTION 1
-pH 7.5
Ca 20 charge
C(4) 20
END
INCREMENTAL_REACTIONS true
USE solution 1
KINETICS
-step 10*360
Aragonite
-M 0
USER_GRAPH
10 GRAPH_X TOTAL_TIME / 3600
20 GRAPH_Y SI("Aragonite")
END


[Isil]

Thank you for the help. I have tested the rates, and with a proper set of constants (for sak and n), the system has not given an error. But I will keep an eye on it.
Thank you