Updated to PHREEQC3 12386

git-svn-id: svn://136.177.114.72/svn_GW/phreeqc3/branches/concrete@12387 1feff8c3-07ed-0310-ac33-dd36852eb9cd

Makefile.am

@@ -42,7 +42,6 @@ phreeqc_SOURCES=\
 	Dictionary.h\
 	dumper.cpp\
 	dumper.h\
-	dw.cpp\
 	Exchange.cxx\
 	Exchange.h\
 	ExchComp.cxx\

Makefile.old

@@ -15,15 +15,16 @@ PROGRAM = phreeqc
 MPI_DIR=/usr/lib64/openmpi
 
 CFG1 :=`uname`
-CFG :=$(shell echo $(CFG1) | sed "s/CYGWIN.*/CYGWIN/")
+CFG :=$(shell echo $(CFG1) | sed "s/Linux.*/Linux/")
-ifeq ($(CFG), CYGWIN)
+ifeq ($(CFG), Linux)
-	SPOOL=>
-	SPOOL2=2>&1
-else
 	SPOOL=>&
 	SPOOL2=
+	CONCAT=;
+else
+	SPOOL=>
+	SPOOL2=2>&1
+	CONCAT=&
 endif
-
 all: release_openmp
 
 Debug: debug
@@ -235,7 +236,6 @@ COMMON_COBJS =  \
 	       cvdense.o \
 	       cvode.o \
 	       dense.o \
-	       dw.o \
 	       gases.o \
 	       input.o \
 	       integrate.o \
@@ -685,14 +685,6 @@ Dictionary.o: ../Dictionary.cpp ../Dictionary.h
 dumper.o: ../dumper.cpp ../dumper.h ../StorageBinList.h \
  ../common/PHRQ_base.h ../common/Parser.h ../common/PHRQ_base.h \
  ../PhreeqcKeywords/Keywords.h ../common/PHRQ_io.h ../common/PHRQ_io.h
-dw.o: ../dw.cpp ../Phreeqc.h ../common/phrqtype.h ../cvdense.h ../cvode.h \
- ../sundialstypes.h ../nvector.h ../dense.h ../smalldense.h ../runner.h \
- ../StorageBinList.h ../common/PHRQ_base.h ../dumper.h \
- ../common/PHRQ_io.h ../PhreeqcKeywords/Keywords.h ../SelectedOutput.h \
- ../NumKeyword.h ../UserPunch.h ../Pressure.h ../cxxMix.h ../Use.h \
- ../Surface.h ../SurfaceComp.h ../NameDouble.h ../common/Parser.h \
- ../common/PHRQ_base.h ../common/PHRQ_io.h ../SurfaceCharge.h \
- ../global_structures.h ../NA.h
 gases.o: ../gases.cpp ../Phreeqc.h ../common/phrqtype.h ../cvdense.h \
  ../cvode.h ../sundialstypes.h ../nvector.h ../dense.h ../smalldense.h \
  ../runner.h ../StorageBinList.h ../common/PHRQ_base.h ../dumper.h \
@@ -1161,8 +1153,9 @@ dependencies:
 	   -I../PhreeqcRM/IPhreeqcPhast/IPhreeqc ../PhreeqcRM/IPhreeqcPhast/IPhreeqc/*.cpp ../PhreeqcRM/IPhreeqcPhast/IPhreeqc/*.c
 
 tester:
-	cd ../mytest; make clean; make -k -j 1 $(SPOOL) make.out $(SPOOL2); make diff $(SPOOL) diff.out $(SPOOL2)
+#	cd ../mytest; make clean; make -k -j 1 $(SPOOL) make.out $(SPOOL2); make diff $(SPOOL) diff.out $(SPOOL2)
-	cd ../examples; make -f Makefile.old clean; make -f Makefile.old -k -j 1 $(SPOOL) make.out $(SPOOL2); make -f Makefile.old diff $(SPOOL) diff.out $(SPOOL2)
+	cd ../mytest $(CONCAT) make clean $(CONCAT) make -k -j 1 $(SPOOL) make.out $(SPOOL2) $(CONCAT) make diff $(SPOOL) diff.out $(SPOOL2)
+	cd ../examples $(CONCAT) make -f Makefile.old clean $(CONCAT) make -f Makefile.old -k -j 1 $(SPOOL) make.out $(SPOOL2) $(CONCAT) make -f Makefile.old diff $(SPOOL) diff.out $(SPOOL2)
 	svn status -q ../mytest 
 	svn status -q ../examples
 

PBasic.cpp

@@ -1529,6 +1529,9 @@ listtokens(FILE * f, tokenrec * l_buf)
 		case tokrho:
 			output_msg("RHO");
 			break;
+		case tokrho_0:
+			output_msg("RHO_0");
+			break;
 /* VP: Density End */
 		case tokcell_volume:
 			output_msg("CELL_VOLUME");
@@ -1548,6 +1551,10 @@ listtokens(FILE * f, tokenrec * l_buf)
 		case tokiso_unit:
 			output_msg("ISO_UNIT");
 			break;
+		case tokkinetics_formula:
+		case tokkinetics_formula_:
+			output_msg("KINETICS_FORMULA$");
+			break;
 		case tokphase_formula:
 		case tokphase_formula_:
 			output_msg("PHASE_FORMULA$");
@@ -1580,6 +1587,9 @@ listtokens(FILE * f, tokenrec * l_buf)
  		case tokvm:
  			output_msg("VM"); // mole volume of aqueous solute
  			break;
+ 		case tokphase_vm:
+ 			output_msg("PHASE_VM"); // mole volume of a phase 
+ 			break;
  		case tokdh_a:
  			output_msg("DH_A"); // Debye-Hueckel A
  			break;
@@ -1629,9 +1639,6 @@ listtokens(FILE * f, tokenrec * l_buf)
 		case tokviscos_0:
 			output_msg("VISCOS_0");
 			break;
-		case tokrho_0:
-			output_msg("RHO_0");
-			break;
 		case tokcurrent_a:
 			output_msg("CURRENT_A");
 			break;
@@ -2946,6 +2953,112 @@ factor(struct LOC_exec * LINK)
 			break;
 		}
 
+	case tokkinetics_formula:
+	case tokkinetics_formula_:
+		{
+			require(toklp, LINK);
+			std::string kinetics_name(stringfactor(STR1, LINK));
+			varrec *elts_varrec = NULL, *coef_varrec = NULL;
+			cxxNameDouble stoichiometry;
+			/*
+			*  Parse arguments
+			*/
+			if (LINK->t != NULL && LINK->t->kind == tokcomma)
+			{
+				/* kinetics_formula("calcite", count, elt, coef) */
+				/* return formula */
+				/*int c; */
+				/*  struct varrec *count_varrec, *names_varrec, *types_varrec, *moles_varrec; */
+				/*  struct varrec *count_varrec, *elt_varrec, *coef_varrec; */
+				/* return number of species */
+				LINK->t = LINK->t->next;
+				count_varrec = LINK->t->UU.vp;
+				if (LINK->t->kind != tokvar || count_varrec->stringvar != 0)
+					snerr(": Cannot find count variable");
+
+				/* return number of names of elements */
+				LINK->t = LINK->t->next;
+				require(tokcomma, LINK);
+				elts_varrec = LINK->t->UU.vp;
+				if (LINK->t->kind != tokvar || elts_varrec->stringvar != 1)
+					snerr(": Cannot find element string variable");
+
+				/* return coefficients of species */
+				LINK->t = LINK->t->next;
+				require(tokcomma, LINK);
+				coef_varrec = LINK->t->UU.vp;
+				if (LINK->t->kind != tokvar || coef_varrec->stringvar != 0)
+					snerr(": Cannot find coefficient variable");
+				LINK->t = LINK->t->next;
+				arg_num = 4;
+			}
+			else
+			{
+				arg_num = 1;
+			}
+			require(tokrp, LINK);
+
+			if (arg_num > 1)
+			{
+				free_dim_stringvar(elts_varrec);
+				PhreeqcPtr->free_check_null(coef_varrec->UU.U0.arr);
+				coef_varrec->UU.U0.arr = NULL;
+			}
+			/*
+			*  Call subroutine
+			*/
+			std::string form = PhreeqcPtr->kinetics_formula(kinetics_name, stoichiometry);
+
+			// put formula as return value
+			n.stringval = true;
+			n.UU.sval = PhreeqcPtr->string_duplicate(form.c_str());
+
+			/*
+			*  fill in varrec structure
+			*/
+
+			if (arg_num > 1)
+			{
+				size_t count = stoichiometry.size();
+				*count_varrec->UU.U0.val = (LDBLE) count;
+				/*
+				* malloc space
+				*/
+				elts_varrec->UU.U1.sarr = (char **) PhreeqcPtr->PHRQ_malloc((count + 1) * sizeof(char *));
+				if (elts_varrec->UU.U1.sarr == NULL)
+					PhreeqcPtr->malloc_error();
+				coef_varrec->UU.U0.arr = (LDBLE *) PhreeqcPtr->PHRQ_malloc((count + 1) * sizeof(LDBLE));
+				if (coef_varrec->UU.U0.arr == NULL)
+					PhreeqcPtr->malloc_error();
+
+				// first position not used
+				elts_varrec->UU.U1.sarr[0] = NULL;
+				coef_varrec->UU.U0.arr[0] = 0;
+
+				// set dims for Basic array
+				for (i = 0; i < maxdims; i++)
+				{
+					elts_varrec->dims[i] = 0;
+					coef_varrec->dims[i] = 0;
+				}
+				// set dims for first dimension and number of dims
+				elts_varrec->dims[0] = (long) (count + 1);
+				coef_varrec->dims[0] = (long) (count + 1);
+				elts_varrec->numdims = 1;
+				coef_varrec->numdims = 1;
+
+				// fill in arrays
+				i = 1;
+				for (cxxNameDouble::iterator it = stoichiometry.begin(); it != stoichiometry.end(); it++)
+				{
+					elts_varrec->UU.U1.sarr[i] = PhreeqcPtr->string_duplicate((it->first).c_str());
+					coef_varrec->UU.U0.arr[i] = it->second;
+					i++;
+				}
+
+			}
+			break;
+		}
 	case tokphase_formula:
 	case tokphase_formula_:
 		{
@@ -3470,6 +3583,9 @@ factor(struct LOC_exec * LINK)
 	case tokrho:
 		n.UU.val = (parse_all) ? 1 : PhreeqcPtr->calc_dens();
 		break;
+	case tokrho_0:
+		n.UU.val = (parse_all) ? 1 : PhreeqcPtr->rho_0;
+		break;
 /* VP: Density End */
 	case tokcell_volume:
 		{
@@ -3552,9 +3668,6 @@ factor(struct LOC_exec * LINK)
  			n.UU.val = (parse_all) ? 1 : PhreeqcPtr->viscos_0;
 		}
  		break;
-	case tokrho_0:
-		n.UU.val = (parse_all) ? 1 : PhreeqcPtr->rho_0;
-		break;
 	case tokcurrent_a:
 		//n.UU.val = (parse_all) ? 1 : PhreeqcPtr->current_x;
 		n.UU.val = (parse_all) ? 1 : PhreeqcPtr->current_A;
@@ -3576,6 +3689,13 @@ factor(struct LOC_exec * LINK)
 			//}
 		}
 		break;
+  	case tokphase_vm:
+		{
+			const char * str = stringfactor(STR1, LINK);
+ 			n.UU.val = (parse_all) ? 1 : PhreeqcPtr->phase_vm(str);
+		}
+ 		break;
+
 	case toksin:
 		n.UU.val = sin(realfactor(LINK));
 		break;
@@ -7178,6 +7298,9 @@ const std::map<const std::string, PBasic::BASIC_TOKEN>::value_type temp_tokens[]
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("viscos",             PBasic::tokviscos),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("viscos_0",           PBasic::tokviscos_0),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("rho_0",              PBasic::tokrho_0),
+	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("kinetics_formula",   PBasic::tokkinetics_formula),
+	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("kinetics_formula$",  PBasic::tokkinetics_formula),
+	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("phase_vm",           PBasic::tokphase_vm),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("current_a",          PBasic::tokcurrent_a),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("pot_v",              PBasic::tokpot_v),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("setdiff_c",          PBasic::toksetdiff_c)

PBasic.h

@@ -279,6 +279,7 @@ public:
 		tokgamma,
 		toklg,
 		tokrho,
+		tokrho_0,
 		tokcell_volume,
 		tokcell_pore_volume,
 		tokcell_porosity,
@@ -294,6 +295,8 @@ public:
 		tokeol_,
 		tokceil,
 		tokfloor,
+		tokkinetics_formula,
+		tokkinetics_formula_,
 		tokphase_formula,
 		tokphase_formula_,
 		tokspecies_formula,
@@ -328,9 +331,9 @@ public:
 		tokedl_species,
 		tokviscos,
 		tokviscos_0,
-		tokrho_0,
 		tokcurrent_a,
-		tokpot_v
+		tokpot_v,
+		tokphase_vm
 	};
 
 #if !defined(PHREEQCI_GUI)

Phreeqc.cpp

@@ -276,6 +276,13 @@ void Phreeqc::init(void)
 	punch.charge_balance	= FALSE;
 	punch.percent_error		= FALSE;
 #endif
+
+	MIN_LM = -30.0;			    /* minimum log molality allowed before molality set to zero */
+	LOG_ZERO_MOLALITY = -30;	/* molalities <= LOG_ZERO_MOLALITY are considered equal to zero */
+	MIN_RELATED_LOG_ACTIVITY = -30;
+	MIN_TOTAL = 1e-25;
+	MIN_TOTAL_SS = MIN_TOTAL/100;
+	MIN_RELATED_SURFACE = MIN_TOTAL*100;
 	// auto Rxn_temperature_map;
 	// auto Rxn_pressure_map;
 
@@ -1173,6 +1180,12 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc)
 	current_mu                      = pSrc->current_mu;
 	mu_terms_in_logk                = pSrc->mu_terms_in_logk;
 
+	MIN_LM = pSrc->MIN_LM;			    /* minimum log molality allowed before molality set to zero */
+	LOG_ZERO_MOLALITY = pSrc->LOG_ZERO_MOLALITY;	/* molalities <= LOG_ZERO_MOLALITY are considered equal to zero */
+	MIN_RELATED_LOG_ACTIVITY = pSrc->MIN_RELATED_LOG_ACTIVITY;
+	MIN_TOTAL = pSrc->MIN_TOTAL;
+	MIN_TOTAL_SS = pSrc->MIN_TOTAL_SS;
+	MIN_RELATED_SURFACE = pSrc->MIN_RELATED_SURFACE;
 	/* ----------------------------------------------------------------------
 	*   STRUCTURES
 	* ---------------------------------------------------------------------- */

Phreeqc.h

@@ -100,6 +100,7 @@ public:
 	LDBLE activity_coefficient(const char *species_name);
 	LDBLE log_activity_coefficient(const char *species_name);
 	LDBLE aqueous_vm(const char *species_name);
+	LDBLE phase_vm(const char *phase_name);
 	LDBLE diff_c(const char *species_name);
 	LDBLE setdiff_c(const char *species_name, double d);
 	LDBLE sa_declercq(double type, double sa, double d, double m, double m0, double gfw);
@@ -148,6 +149,7 @@ public:
 	static int system_species_compare(const void *ptr1, const void *ptr2);
 	LDBLE system_total(const char *total_name, LDBLE * count, char ***names,
 		char ***types, LDBLE ** moles);
+	std::string kinetics_formula(std::string kinetics_name, cxxNameDouble &stoichiometry);
 	std::string phase_formula(std::string phase_name, cxxNameDouble &stoichiometry);
 	std::string species_formula(std::string phase_name, cxxNameDouble &stoichiometry);
 	LDBLE list_ss(std::string ss_name, cxxNameDouble &composition);
@@ -158,6 +160,7 @@ public:
 	int system_total_surf(void);
 	int system_total_gas(void);
 	int system_total_equi(void);
+	int system_total_kin(void);
 	int system_total_ss(void);
 	int system_total_elt(const char *total_name);
 	int system_total_elt_secondary(const char *total_name);
@@ -1605,6 +1608,13 @@ protected:
 	int stop_program;
 	int incremental_reactions;
 
+	double MIN_LM;
+	double LOG_ZERO_MOLALITY;
+	double MIN_TOTAL;
+	double MIN_TOTAL_SS;
+	double MIN_RELATED_SURFACE;
+	double MIN_RELATED_LOG_ACTIVITY;
+
 	int count_strings;
 	int max_strings;
 
@@ -2008,7 +2018,11 @@ public:
 #endif /*defined (PHREEQ98) || defined (_MSC_VER)*/
 
 #if defined(HAVE_ISFINITE)
-#  define PHR_ISFINITE(x) isfinite(x)
+#  if __GNUC__ && (__cplusplus >= 201103L)
+#    define PHR_ISFINITE(x) std::isfinite(x)
+#  else
+#    define PHR_ISFINITE(x) isfinite(x)
+#  endif
 #elif defined(HAVE_FINITE)
 #  define PHR_ISFINITE(x) finite(x)
 #elif defined(HAVE_ISNAN)

Solution.cxx

@@ -1120,6 +1120,7 @@ cxxSolution::Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge, const cxxNameDouble
 	this->total_h = h_tot;
 	this->total_o = o_tot;
 	this->cb = charge;
+	this->mass_water = o_tot / 55.5;
 
 	// Don`t bother to update activities?
 	this->Update(const_nd);

basicsubs.cpp

@@ -98,7 +98,26 @@ aqueous_vm(const char *species_name)
 	}
 	return (g);
 }
+/* ---------------------------------------------------------------------- */
+LDBLE Phreeqc::
+phase_vm(const char *phase_name)
+/* ---------------------------------------------------------------------- */
+{
+	struct phase *phase_ptr;
+	int l;
+	LDBLE g;
 
+	phase_ptr = phase_bsearch(phase_name, &l, FALSE);
+	if (phase_ptr == NULL)
+	{
+		g = 0.0;
+	}
+	else
+	{
+		g = phase_ptr->logk[vm0];
+	}
+	return (g);
+}
 /* ---------------------------------------------------------------------- */
 LDBLE Phreeqc::
 sa_declercq(double sa_type, double Sa, double d, double m, double m0, double gfw)
@@ -2954,6 +2973,10 @@ system_total(const char *total_name, LDBLE * count, char ***names,
 	{
 		system_total_equi();
 	}
+	else if (strcmp_nocase(total_name, "kin") == 0)
+	{
+		system_total_kin();
+	}
 	else
 	{
 		if (strstr(total_name, "(") == NULL)
@@ -3014,6 +3037,64 @@ system_total(const char *total_name, LDBLE * count, char ***names,
 	return (sys_tot);
 }
 
+/* ---------------------------------------------------------------------- */
+std::string Phreeqc::
+kinetics_formula(std::string kin_name, cxxNameDouble &stoichiometry)
+/* ---------------------------------------------------------------------- */
+{
+/*
+ *   Returns formula of kinetic reactant
+ *   Also returns arrays of elements and stoichiometry in stoichiometry
+ */
+	stoichiometry.clear();
+	std::string formula;
+
+	if (use.Get_kinetics_ptr() == NULL)
+		return (formula);
+	std::vector <cxxKineticsComp> comps = use.Get_kinetics_ptr()->Get_kinetics_comps();
+	count_elts = 0;
+	paren_count = 0;
+	for  (size_t i=0 ; i < comps.size(); i++)
+	{
+			cxxKineticsComp *comp_ptr = &comps[i]; 
+			if (kin_name == comp_ptr->Get_rate_name().c_str())
+			{
+				cxxNameDouble nd = comp_ptr->Get_namecoef();
+				cxxNameDouble::iterator it = nd.begin();
+				for ( ; it != nd.end(); it++)
+				{
+					// Try Phases
+					int l; 
+					struct phase *phase_ptr = phase_bsearch(it->first.c_str(), &l, FALSE);
+					if (phase_ptr != NULL)
+					{
+						add_elt_list(phase_ptr->next_elt, it->second);
+					}
+					else
+					{
+						// add formula
+						std::string name = it->first;
+						LDBLE coef = it->second;
+						char * temp_name = string_duplicate(name.c_str());
+						char *ptr = temp_name;
+						get_elts_in_species(&ptr, coef);
+						free_check_null(temp_name);
+					}
+				}
+				formula.append(kin_name);
+				//elt_list[count_elts].elt = NULL;
+				if (count_elts > 0)
+				{
+					qsort(elt_list, (size_t) count_elts,
+						(size_t) sizeof(struct elt_list), elt_list_compare);
+					elt_list_combine();
+				}
+				stoichiometry = elt_list_NameDouble();
+				break;
+			}
+	}
+	return (formula);
+}
 /* ---------------------------------------------------------------------- */
 std::string Phreeqc::
 phase_formula(std::string phase_name, cxxNameDouble &stoichiometry)
@@ -3363,7 +3444,8 @@ system_total_equi(void)
 			int l;
 			struct phase *phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE);
 			sys[count_sys].name = string_duplicate(phase_ptr->name);
-			sys[count_sys].moles = comp_ptr->Get_moles();
+			//sys[count_sys].moles = comp_ptr->Get_moles();
+			sys[count_sys].moles = equi_phase(sys[count_sys].name);
 			sys_tot += sys[count_sys].moles;
 			sys[count_sys].type = string_duplicate("equi");
 			count_sys++;
@@ -3374,6 +3456,30 @@ system_total_equi(void)
 }
 /* ---------------------------------------------------------------------- */
 int Phreeqc::
+system_total_kin(void)
+/* ---------------------------------------------------------------------- */
+{
+/*
+ *  Equilibrium phases
+ */
+	if (use.Get_kinetics_ptr() == NULL)
+		return (OK);
+	std::vector <cxxKineticsComp> comps = use.Get_kinetics_ptr()->Get_kinetics_comps();
+	for  (size_t i=0 ; i < comps.size(); i++)
+	{
+			cxxKineticsComp *comp_ptr = &comps[i]; 
+			sys[count_sys].name = string_duplicate(comp_ptr->Get_rate_name().c_str());
+			sys[count_sys].moles = comp_ptr->Get_m();
+			sys_tot += sys[count_sys].moles;
+			sys[count_sys].type = string_duplicate("kin");
+			count_sys++;
+			space((void **) ((void *) &sys), count_sys, &max_sys,
+				  sizeof(struct system_species));
+	}
+	return (OK);
+}
+/* ---------------------------------------------------------------------- */
+int Phreeqc::
 system_total_ss(void)
 /* ---------------------------------------------------------------------- */
 {

dw.cpp

@@ -1,415 +0,0 @@
-#include "Phreeqc.h"
-
-#ifdef SKIP
-/* ---------------------------------------------------------------------- */
-int Phreeqc::
-DW(LDBLE T)
-/* ---------------------------------------------------------------------- */
-/*
-C
-C      SUBROUTINE TO CALCULATE THE DENSITY OF WATER AS A FUNCTION OF
-C      TEMPERATURE.  T IS IN KELVIN, P IS IN PASCALS, DW0 IS IN G/CM^3
-C
-C     FROM L. HAAR, J. S. GALLAGHER, AND G. S. KELL, (1984)
-C
-*/
-{
-	LDBLE FP = 9.869232667e0, P, DGSS, D;
-
-	BB(T);
-	P = 1.0e0 / FP;
-	if (T > 373.149e0)
-		P = PS(T);
-	DGSS = P / T / .4e0;
-	if (T < TZ)
-	{
-		DGSS = 1.0e0 / (VLEST(T));
-	}
-	DFIND(&D, P, DGSS, T);
-	DW0 = D;
-	VP = P * FP;
-	return OK;
-}
-
-/* ---------------------------------------------------------------------- */
- int Phreeqc::
-BB(LDBLE T)
-/* ---------------------------------------------------------------------- */
-/*
-
-C
-C     THIS SUBROUTINE CALCULATES THE B'S NEEDED FOR FUNCTION DW.
-C     THE B'S CALCULATED HERE ARE IN CM3/G.
-C
-C     FROM L. HAAR, J. S. GALLAGHER, AND G. S. KELL, (1984)
-C
-*/
-{
-	LDBLE V[11];
-	int I;
-	/* COMMON /BCONST/ */
-	LDBLE P[11] =
-		{ 0, 0.7478629e0, -.3540782e0, 0.e0, 0e0, .007159876e0, 0.e0,
-		-.003528426e0, 0., 0., 0.
-	};
-	LDBLE Q[11] =
-		{ 0, 1.1278334e0, 0.e0, -.5944001e0, -5.010996e0, 0.e0, .63684256e0,
-		0., 0., 0., 0.
-	};
-
-	V[1] = 1.0;
-	for (I = 2; I <= 10; I++)
-	{
-		V[I] = V[I - 1] * TZ / T;
-	}
-	B1 = P[1] + P[2] * log(1.e0 / V[2]);
-	B2 = Q[1];
-	B1T = P[2] * V[2] / TZ;
-	B2T = 0.e0;
-	B1TT = 0.e0;
-	B2TT = 0.e0;
-	for (I = 3; I <= 10; I++)
-	{
-		B1 = B1 + P[I] * V[I - 1];
-		B2 = B2 + Q[I] * V[I - 1];
-		B1T = B1T - (I - 2) * P[I] * V[I - 1] / T;
-		B2T = B2T - (I - 2) * Q[I] * V[I - 1] / T;
-		B1TT = B1TT + P[I] * (I - 2) * (I - 2) * V[I - 1] / T / T;
-		B2TT = B2TT + Q[I] * (I - 2) * (I - 2) * V[I - 1] / T / T;
-	}
-	B1TT = B1TT - B1T / T;
-	B2TT = B2TT - B2T / T;
-	return OK;
-}
-
-/* ---------------------------------------------------------------------- */
- LDBLE Phreeqc::
-PS(LDBLE T)
-/* ---------------------------------------------------------------------- */
-/*
-C
-C     THIS FUNCTION CALCULATES AN APPROXIMATION TO THE VAPOR PRESSURE, P
-C     AS A FUNCTION OF THE INPUT TEMPERATURE. THE VAPOR PRESSURE
-C     CALCULATED AGREES WITH THE VAPOR PRESSURE PREDICTED BY THE SURFACE
-C     TO WITHIN .02% TO WITHIN A DEGREE OR SO OF THE CRITICAL TEMPERATUR
-C     AND CAN SERVE AS AN INITIAL GUESS FOR FURTHER REFINEMENT BY
-C     IMPOSING THE CONDITION THAT GL=GV.
-C
-C     FROM L. HAAR, J. S. GALLAGHER, AND G. S. KELL, (1984)
-C
-*/
-{
-	LDBLE A[9] = { 0, -7.8889166e0, 2.5514255e0, -6.716169e0,
-		33.239495e0, -105.38479e0, 174.35319e0, -148.39348e0,
-		48.631602e0
-	};
-	LDBLE PL, V, W, B, L_Z, Q;
-	int I;
-	if (T <= 314.e0)
-	{
-		PL = 6.3573118e0 - 8858.843e0 / T + 607.56335e0 * pow(T, (LDBLE) -.6e0);
-		return (.1e0 * exp(PL));
-	}
-	V = T / 647.25e0;
-	W = fabs(1.e0 - V);
-	B = 0.e0;
-	for (I = 1; I <= 8; I++)
-	{
-		L_Z = I;
-		B = B + A[I] * pow(W, ((L_Z + 1.e0) / 2.e0));
-	}
-	Q = B / V;
-	return (22.093e0 * exp(Q));
-}
-
-/* ---------------------------------------------------------------------- */
- LDBLE Phreeqc::
-VLEST(LDBLE T)
-/* ---------------------------------------------------------------------- */
-/*
-C
-C     FROM L. HAAR, J. S. GALLAGHER, AND G. S. KELL, (1984)
-C
-*/
-{
-	LDBLE A = -1.59259e1, B = 6.57886e-2, C = -1.12666e-4, D = 7.33191e-8,
-		E = 1.60229e3, F = 2.88572e0, G = 650.0e0;
-
-	return (A + B * T + C * T * T + D * T * T * T + E / T + F / (G - T));
-}
-
-/* ---------------------------------------------------------------------- */
- int Phreeqc::
-DFIND(LDBLE * DOUT, LDBLE P, LDBLE D, LDBLE T)
-/* ---------------------------------------------------------------------- */
-/*
-C
-C     ROUTINE TO FIND DENSITY CORRESPONDING TO INPUT PRESSURE P(MPA), AN
-C     TEMPERATURE T(K), USING INITIAL GUESS DENSITY D(G/CM3). THE OUTPUT
-C     DENSITY IS IN G/CM3.
-C
-C     FROM L. HAAR, J. S. GALLAGHER, AND G. S. KELL, (1984)
-C
-*/
-{
-	int L;
-	LDBLE DD, RT, PP_dfind, DPD, DPDX, DP, X;
-	/*      LDBLE DD, RT, PP, DPD, DPDX, DP, X; */
-
-	DD = D;
-	RT = GASCON * T;
-	if (DD <= 0.e0)
-		DD = 1.e-8;
-	if (DD > 1.9e0)
-		DD = 1.9e0;
-	L = 0;
-	for (L = 1; L <= 30; L++)
-	{
-		if (DD <= 0.e0)
-			DD = 1.e-8;
-		if (DD > 1.9e0)
-			DD = 1.9e0;
-		QQ(T, DD);
-		PP_dfind = RT * DD * BASE(DD) + Q0;
-		DPD = RT * (Z + Y * DZ) + Q5;
-		/*
-		   C
-		   C  THE FOLLOWING 3 LINES CHECK FOR NEGATIVE DP/DRHO, AND IF SO ASSUME
-		   C    GUESS TO BE IN 2-PHASE REGION, AND CORRECT GUESS ACCORDINGLY.
-		   C
-		 */
-		if (DPD <= 0.e0)
-		{
-			if (D >= .2967e0)
-				DD = DD * 1.02e0;
-			if (D < .2967e0)
-				DD = DD * .98e0;
-			if (L <= 10)
-				continue;
-		}
-		else
-		{
-/* 13 */
-			DPDX = DPD * 1.1e0;
-			if (DPDX < 0.1e0)
-				DPDX = 0.1e0;
-			DP = fabs(1.e0 - PP_dfind / P);
-			if (DP < 1.e-8)
-				break;
-			if (D > .3e0 && DP < 1.e-7)
-				break;
-			if (D > .7e0 && DP < 1.e-6)
-				break;
-			X = (P - PP_dfind) / DPDX;
-			if (fabs(X) > .1e0)
-				X = X * .1e0 / fabs(X);
-			DD = DD + X;
-			if (DD < 0.e0)
-				DD = 1.e-8;
-		}
-	}
-	if (L > 30)
-		error_msg("In subroutine DFIND", STOP);
-	//cite Jonathan Toner, remove error_msg
-	//DD = 1;
-/*   20 CONTINUE */
-	*DOUT = DD;
-	return OK;
-}
-
-/* ---------------------------------------------------------------------- */
- int Phreeqc::
-QQ(LDBLE T, LDBLE D)
-/* ---------------------------------------------------------------------- */
-/*
-C
-C     THIS ROUTINE CALCULATES, FOR A GIVEN T(K) AND D(G/CM3), THE RESIDUL
-C     CONTRIBUTIONS TO: PRESSURE (Q), DP/DRHO (Q5)
-C     THIS SUBROUTINE IS USED IN DENSITY OF WATER CALCULATION.
-C
-C     FROM L. HAAR, J. S. GALLAGHER, AND G. S. KELL, (1984)
-C
-*/
-{
-	/* COMMON /NCONST/ */
-	LDBLE G[41] =
-		{ 0, -.53062968529023e3, .22744901424408e4, .78779333020687e3,
-		-.69830527374994e2, .17863832875422e5, -.39514731563338e5,
-		.33803884280753e5,
-		-.13855050202703e5, -.25637436613260e6, .48212575981415e6,
-		-.34183016969660e6,
-		.12223156417448e6, .11797433655832e7, -.21734810110373e7,
-		.10829952168620e7,
-		-.25441998064049e6, -.31377774947767e7, .52911910757704e7,
-		-.13802577177877e7,
-		-.25109914369001e6, .46561826115608e7, -.72752773275387e7,
-		.41774246148294e6,
-		.14016358244614e7, -.31555231392127e7, .47929666384584e7,
-		.40912664781209e6,
-		-.13626369388386e7, .69625220862664e6, -.10834900096447e7,
-		-.22722827401688e6,
-		.38365486000660e6, .68833257944332e4, .21757245522644e5,
-		-.26627944829770e4,
-		-.70730418082074e5, -.225e0, -1.68e0, .055e0, -93.0e0
-	};
-	/*int II[41]={0,4*0,4*1,4*2,4*3,4*4,4*5,4*6,4*8,2*2,0,4,3*2,4}; */
-	int II[41] =
-		{ 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5,
-		5,
-		5, 6, 6, 6, 6, 8, 8, 8, 8, 2, 2, 0, 4, 2, 2, 2, 4
-	};
-	/*int JJ[41]={0,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,1,3*4,0,2,0,0}; */
-	int JJ[41] =
-		{ 0, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3,
-		5,
-		7, 2, 3, 5, 7, 2, 3, 5, 7, 1, 4, 4, 4, 0, 2, 0, 0
-	};
-	int NC = 36;
-	/* COMMON /ADDCON/ */
-	LDBLE ATZ[5] = { 0, 64.e1, 64.e1, 641.6e0, 27.e1 }, ADZ[5] =
-	{
-	0, .319e0, .319e0, .319e0, 1.55e0}, AAT[5] =
-	{
-	0, 2.e4, 2.e4, 4.e4, 25.e0}, AAD[5] =
-	{
-	0, 34.e0, 4.e1, 3.e1, 1.05e3};
-	LDBLE *QZT;
-	LDBLE QR[12], QT[11] /*, QZT[10] */ ;
-	/*EQUIVALENCE (QT(2),QZT(1)) */
-
-	LDBLE E, Q10, Q20, V, QP, DDZ, DEL, EX1, DEX, ATT, TX,
-		TAU, EX2, TEX, QM, FCT, Q5T;
-	int I, K, L, J, KM;
-	QZT = &(QT[1]);
-	QR[1] = 0.e0;
-	Q5 = 0.e0;
-	Q0 = 0.e0;
-	E = exp(-AA * D);
-	Q10 = D * D * E;
-	Q20 = 1.e0 - E;
-	QR[2] = Q10;
-	V = TZ / T;
-	QT[1] = T / TZ;
-	/*DO 4 I=2,10 */
-	for (I = 2; I <= 10; I++)
-	{
-		QR[I + 1] = QR[I] * Q20;
-		/*   4  QT[I]=QT[I-1]*V */
-		QT[I] = QT[I - 1] * V;
-	}
-	/* DO 10 I=1,NC */
-	for (I = 1; I <= NC; I++)
-	{
-		K = II[I] + 1;
-		L = JJ[I];
-		QP = G[I] * AA * QR[K + 1] * QZT[L];
-		Q0 = Q0 + QP;
-		/*10 Q5 = Q5 + AA*(2.e0/D-AA*(1.e0-E*(K-1)/Q20))*QP */
-		Q5 = Q5 + AA * (2.e0 / D - AA * (1.e0 - E * (K - 1) / Q20)) * QP;
-	}
-	QP = 0.e0;
-	/* DO 20 J=37,40 */
-	for (J = 37; J <= 40; J++)
-	{
-		if (fabs(G[J]) < 1.0e-20)
-			continue;
-		K = II[J];
-		KM = JJ[J];
-		DDZ = ADZ[J - 36];
-		DEL = D / DDZ - 1.e0;
-		if (fabs(DEL) < 1.e-10)
-			DEL = 1.e-10;
-		EX1 = -AAD[J - 36] * pow(DEL, K);
-		if (EX1 <= -88.028e0)
-		{
-			DEX = 0.e0;
-		}
-		else
-		{
-			DEX = exp(EX1) * pow(DEL, KM);
-		}
-		ATT = AAT[J - 36];
-		TX = ATZ[J - 36];
-		TAU = T / TX - 1.e0;
-		EX2 = -ATT * TAU * TAU;
-		if (EX2 <= -88.028e0)
-		{
-			TEX = 0.e0;
-		}
-		else
-		{
-			TEX = exp(EX2);
-		}
-		Q10 = DEX * TEX;
-		QM = KM / DEL - K * AAD[J - 36] * pow(DEL, (K - 1));
-		FCT = QM * D * D * Q10 / DDZ;
-		Q5T =
-			FCT * (2.e0 / D + QM / DDZ) - pow((D / DDZ),
-											  2) * Q10 * (KM / DEL / DEL +
-														  K * (K -
-															   1) * AAD[J -
-																		36] *
-														  pow(DEL, (K - 2)));
-		Q5 = Q5 + Q5T * G[J];
-		QP = QP + G[J] * FCT;
-		/*  20  CONTINUE */
-	}
-	Q0 = Q0 + QP;
-	return OK;
-}
-
-/* ---------------------------------------------------------------------- */
- LDBLE Phreeqc::
-BASE(LDBLE D)
-/* ---------------------------------------------------------------------- */
-/*
-C
-C     THIS FUNCTION CALCULATES THE Z (=PBASE/(DRT)) NEEDED FOR FUNCTION
-C
-C     FROM L. HAAR, J. S. GALLAGHER, AND G. S. KELL, (1984)
-C
-*/
-{
-	LDBLE X, Z0, DZ0;
-/*
-C
-C     G1,G2 AND GF ARE THE ALPHA, BETA AND GAMMA FOR DENSITY OF WATER
-C     CALCULATIONS.  B1 AND B2 ARE THE 'EXCLUDED VOLUME' AND '2ND VIRIAL
-C     SUPPLIED BY THE SUBROUTINE BB(T), WHICH ALSO SUPPLIES THE 1ST AND
-C     2ND DERIVATIVES WITH RESPECT TO T (B1T,B2T,B1TT,B2TT).
-C
-*/
-	Y = .25e0 * B1 * D;
-	X = 1.e0 - Y;
-	Z0 = (1.e0 + G1 * Y + G2 * Y * Y) / pow(X, 3);
-	Z = Z0 + 4.e0 * Y * (B2 / B1 - GF);
-	DZ0 =
-		(G1 + 2.e0 * G2 * Y) / pow(X,
-								   3) + 3.e0 * (1.e0 + G1 * Y +
-												G2 * Y * Y) / pow(X, 4);
-	DZ = DZ0 + 4.e0 * (B2 / B1 - GF);
-	return (Z);
-}
-
-/* ---------------------------------------------------------------------- */
-LDBLE Phreeqc::
-DC(LDBLE T)
-/* ---------------------------------------------------------------------- */
-/*
-C
-C     THIS FUNCTION CALCULATES THE RELATIVE DIELECTRIC CONSTANT AS A
-C     FUNCTION OF TEMPERATURE, ASSUMING ONE ATMOSPHERE PRESSURE
-C     ACCORDING TO D. J. BRADLEY AND K. S. PITZER, (1979)
-C
-*/
-{
-	LDBLE D1000, C, B;
-	LDBLE U[10] = { 0, 3.4279e2, -5.0866e-3, 9.4690e-7, -2.0525e0, 3.1159e3,
-		-1.8289e2, -8.0325e3, 4.2142e6, 2.1417e0
-	};
-	D1000 = U[1] * exp(U[2] * T + U[3] * T * T);
-	C = U[4] + U[5] / (U[6] + T);
-	B = U[7] + U[8] / T + U[9] * T;
-	return (D1000 + C * log((B + VP * 1.01325e0) / (B + 1000.0e0)));
-}
-#endif

global_structures.h

@@ -132,19 +132,19 @@
 #define MAX_LM 3.0				/* maximum log molality allowed in intermediate iterations */
 #define MAX_M 1000.0
 #ifdef USE_DECIMAL128
-#define MIN_LM -80.0			/* minimum log molality allowed before molality set to zero */
+//#define MIN_LM -80.0			/* minimum log molality allowed before molality set to zero */
-#define LOG_ZERO_MOLALITY -80	/* molalities <= LOG_ZERO_MOLALITY are considered equal to zero */
+//#define LOG_ZERO_MOLALITY -80	/* molalities <= LOG_ZERO_MOLALITY are considered equal to zero */
-#define MIN_TOTAL 1e-60
+//#define MIN_TOTAL 1e-60
-#define MIN_TOTAL_SS MIN_TOTAL/100
+//#define MIN_TOTAL_SS MIN_TOTAL/100
-#define MIN_RELATED_SURFACE MIN_TOTAL*100
+//#define MIN_RELATED_SURFACE MIN_TOTAL*100
-#define MIN_RELATED_LOG_ACTIVITY -60
+//#define MIN_RELATED_LOG_ACTIVITY -60
 #else
-#define MIN_LM -30.0			/* minimum log molality allowed before molality set to zero */
+//#define MIN_LM -30.0			/* minimum log molality allowed before molality set to zero */
-#define LOG_ZERO_MOLALITY -30	/* molalities <= LOG_ZERO_MOLALITY are considered equal to zero */
+//#define LOG_ZERO_MOLALITY -30	/* molalities <= LOG_ZERO_MOLALITY are considered equal to zero */
-#define MIN_TOTAL 1e-25
+//#define MIN_TOTAL 1e-25
-#define MIN_TOTAL_SS MIN_TOTAL/100
+//#define MIN_TOTAL_SS MIN_TOTAL/100
-#define MIN_RELATED_SURFACE MIN_TOTAL*100
+//#define MIN_RELATED_SURFACE MIN_TOTAL*100
-#define MIN_RELATED_LOG_ACTIVITY -30
+//#define MIN_RELATED_LOG_ACTIVITY -30
 #endif
 #define REF_PRES_PASCAL 1.01325E5   /* Reference pressure: 1 atm */
 /*

kinetics.cpp

@@ -230,8 +230,12 @@ RESTART:   // if limiting rates, jump to here
 						if (0.9 * surface_comp_ptr->Get_phase_proportion() *
 							(kinetics_comp_ptr->Get_m()) < MIN_RELATED_SURFACE)
 						{
-							master_ptr = master_bsearch(ptr);
+							//master_ptr = master_bsearch(ptr);
-							master_ptr->total = 0.0;
+							master_ptr = master_bsearch(surface_comp_ptr->Get_master_element().c_str());
+							if (master_ptr != NULL) 
+							{
+								master_ptr->total = 0.0;
+							}
 						}
 						else
 						{
@@ -1249,9 +1253,15 @@ set_and_run_wrapper(int i, int use_mix, int use_kinetics, int nsaver,
 	int old_diag, old_itmax;
 	LDBLE old_tol, old_min_value, old_step, old_pe, old_pp_column_scale;
 	LDBLE small_pe_step, small_step;
+#if (__GNUC__ && (__cplusplus >= 201103L)) || (_MSC_VER >= 1600)
+	std::unique_ptr<cxxPPassemblage> pp_assemblage_save=NULL;
+	std::unique_ptr<cxxSSassemblage> ss_assemblage_save=NULL;
+	std::unique_ptr<cxxKinetics> kinetics_save=NULL;
+#else
 	std::auto_ptr<cxxPPassemblage> pp_assemblage_save(NULL);
 	std::auto_ptr<cxxSSassemblage> ss_assemblage_save(NULL);
 	std::auto_ptr<cxxKinetics> kinetics_save(NULL);
+#endif
 
 	
 	small_pe_step = 5.;

mainsubs.cpp

@@ -2391,7 +2391,7 @@ run_simulations(void)
 {
 	char token[MAX_LENGTH];
 //#ifdef SKIP_KEEP
-#if defined(WIN32)
+#if defined(_MSC_VER) && (_MSC_VER < 1900)  // removed in vs2015
 	unsigned int old_exponent_format;
 	old_exponent_format = _set_output_format(_TWO_DIGIT_EXPONENT);
 #endif

model.cpp

@@ -1679,19 +1679,19 @@ ineq(int in_kode)
 	}
 #define SHRINK_ARRAY
 #ifdef SHRINK_ARRAY
-	if (sit_model && full_pitzer == FALSE)
+	if ((sit_model || pitzer_model) && full_pitzer == FALSE)
 	{
 		n = count_unknowns - (int) s_list.size();
 		for (int i = 0; i < l_count_rows; i++)
 		{
-			//for (int j = 0; j < n; j++)
+			for (int j = 0; j < n; j++)
-			//{
-			//	ineq_array[i*(n+2) + j] = ineq_array[i*(count_unknowns+2) +j];
-			//}
-			if (i > 0)
 			{
-				memcpy((void *) &ineq_array[i*(n+2)], (void *) &ineq_array[i*(count_unknowns+2)], (size_t) (n) * sizeof(LDBLE));
+				ineq_array[i*(n+2) + j] = ineq_array[i*(count_unknowns+2) +j];
 			}
+			//if (i > 0)
+			//{
+			//	memcpy((void *) &ineq_array[i*(n+2)], (void *) &ineq_array[i*(count_unknowns+2)], (size_t) (n) * sizeof(LDBLE));
+			//}
 			ineq_array[i*(n+2) + n] = ineq_array[i*(count_unknowns+2) + count_unknowns];
 		}
 	}

pitzer.cpp

@@ -1812,12 +1812,14 @@ pitzer_revise_guesses(void)
 	int l_iter, max_iter, repeat, fail;
 	LDBLE weight, f;
 
-	max_iter = 10;
+	max_iter = 100;
 	/* gammas(mu_x); */
 	l_iter = 0;
 	repeat = TRUE;
-	fail = FALSE;;
+	fail = FALSE;
-	while (repeat == TRUE)
+	double d = 2;
+	double logd = log10(d);
+	while (repeat == TRUE && fail == FALSE)
 	{
 		l_iter++;
 		if (debug_set == TRUE)
@@ -1892,18 +1894,14 @@ pitzer_revise_guesses(void)
 				else if (f == 0)
 				{
 					repeat = TRUE;
-					x[i]->master[0]->s->la += 5;
+					x[i]->master[0]->s->la += logd;
 /*!!!!*/ if (x[i]->master[0]->s->la < -999.)
 						x[i]->master[0]->s->la = MIN_RELATED_LOG_ACTIVITY;
 				}
-				else if (fail == TRUE && f < 1.5 * fabs(x[i]->moles))
+				else if (f > d * fabs(x[i]->moles)
+					|| f < 1.0/d * fabs(x[i]->moles))
 				{
-					continue;
+					weight = (f < 1.0/d * fabs(x[i]->moles)) ? 0.3 : 1.0;
-				}
-				else if (f > 1.5 * fabs(x[i]->moles)
-						 || f < 1e-5 * fabs(x[i]->moles))
-				{
-					weight = (f < 1e-5 * fabs(x[i]->moles)) ? 0.3 : 1.0;
 					if (x[i]->moles <= 0)
 					{
 						x[i]->master[0]->s->la = MIN_RELATED_LOG_ACTIVITY;
@@ -1932,10 +1930,10 @@ pitzer_revise_guesses(void)
 					continue;
 				}
 				if (f > 1.5 * fabs(x[i]->moles)
-					|| f < 1e-5 * fabs(x[i]->moles))
+					|| f < 1.0/d * fabs(x[i]->moles))
 				{
 					repeat = TRUE;
-					weight = (f < 1e-5 * fabs(x[i]->moles)) ? 0.3 : 1.0;
+					weight = (f < 1.0/d * fabs(x[i]->moles)) ? 0.3 : 1.0;
 					x[i]->master[0]->s->la += weight *
 						log10(fabs(x[i]->moles / x[i]->sum));
 					if (debug_set == TRUE)
@@ -2266,7 +2264,10 @@ model_pz(void)
 			{
 				full_pitzer = FALSE;
 			}
-			molalities(TRUE);
+			if (molalities(FALSE) == ERROR)
+			{
+				pitzer_revise_guesses();
+			}
 			if (use.Get_surface_ptr() != NULL &&
 				use.Get_surface_ptr()->Get_dl_type() != cxxSurface::NO_DL &&
 				use.Get_surface_ptr()->Get_related_phases())
@@ -2274,6 +2275,17 @@ model_pz(void)
 			mb_sums();
 			mb_gases();
 			mb_ss();
+/*
+ *   Switch bases if necessary
+ */
+			if (switch_bases() == TRUE)
+			{
+				
+				count_basis_change++;
+				count_unknowns -= (int) s_list.size();
+				reprep();
+				full_pitzer = false;
+			}
 			/* debug
 			   species_list_sort();
 			   sum_species();

prep.cpp

@@ -4099,6 +4099,7 @@ calc_PR(std::vector<struct phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
 			P = 1.;
 	} else
 	{
+		if (P < 1e-10) P = 1e-10;
 		r3[1] = b_sum - R_TK / P;
 		r3_12 = r3[1] * r3[1];
 		r3[2] = -3.0 * b2 + (a_aa_sum - R_TK * 2.0 * b_sum) / P;

print.cpp

@@ -3185,7 +3185,7 @@ punch_identifiers(void)
 					(double) (100 * cb_x / total_ions_x));
 		}
 	}
-
+	punch_flush();
 	return (OK);
 }
 

read.cpp

@@ -4807,12 +4807,17 @@ read_selected_output(void)
 	else if (n_user == 1 && so == SelectedOutput_map.end())
 	{
 		// n_user = 1, new definition, do nothing use; constructor default
+		temp_selected_output.Set_new_def(true);
 	}
 	else 
 	{
 		// n_user != 1 then reset false
 
 		temp_selected_output.Reset(false);
+		if (so == SelectedOutput_map.end())
+		{
+			temp_selected_output.Set_new_def(true);
+		}
 	}
 
 	CParser parser(this->phrq_io);
@@ -8346,9 +8351,11 @@ read_debug(void)
 		"try",						       /* 17 */
 		"numerical_fixed_volume",          /* 18 */
 		"force_numerical_fixed_volume",    /* 19 */
-		"equi_delay"                       /* 20 */
+		"equi_delay",                      /* 20 */
+		"minimum_total",                   /* 21 */  
+		"min_total"                        /* 22 */   
 	};
-	int count_opt_list = 21;
+	int count_opt_list = 23;
 /*
  *   Read parameters:
  *	ineq_tol;
@@ -8444,6 +8451,12 @@ read_debug(void)
 		case 20:				/* equi_delay */
 			sscanf(next_char, "%d", &equi_delay);
 			break;
+		case 21:				/* minimum_total */
+		case 22:				/* min_total */
+			sscanf(next_char, SCANFORMAT, &MIN_TOTAL);
+			MIN_TOTAL_SS = MIN_TOTAL/100;
+			MIN_RELATED_SURFACE = MIN_TOTAL*100;
+			break;
 		}
 		if (return_value == EOF || return_value == KEYWORD)
 			break;

tally.cpp

@@ -397,6 +397,7 @@ zero_tally_table(void)
 	int i, j, k;
 	for (i = 0; i < count_tally_table_columns; i++)
 	{
+		tally_table[i].moles = 0.0;
 		for (j = 0; j < count_tally_table_rows; j++)
 		{
 			for (k = 0; k < 3; k++)