Merge commit '8494cef9c5aca1120ce55dbf0be46f87beb49b96'

src/phreeqcpp/ChartObject.cpp

@@ -11,7 +11,8 @@
 #include "ChartObject.h"
 #include "Parser.h"
 #include <fstream>
-#include <math.h>
+//#include <math.h>
+#include <cmath>
 #include <iomanip>
 #include "phqalloc.h"
 

src/phreeqcpp/PBasic.cpp

@@ -7215,7 +7215,7 @@ _NilCheck(void)
 	return _Escape(-3);
 }
 
-#ifdef SKIP
+#ifdef NPP
 /* The following is suitable for the HP Pascal operating system.
    It might want to be revised when emulating another system. */
 

src/phreeqcpp/PBasic.h

@@ -7,7 +7,8 @@
 #include <stdio.h>
 #include <limits.h>
 #include <ctype.h>
-#include <math.h>
+//#include <math.h>
+#include <cmath>
 #include <setjmp.h>
 #include "phrqtype.h"
 #include "PHRQ_base.h"

src/phreeqcpp/Phreeqc.h

@@ -283,12 +283,12 @@ public:
 	int sum_diffuse_layer(cxxSurfaceCharge* surface_charge_ptr1);
 	int calc_all_donnan(void);
 	int calc_init_donnan(void);
+	LDBLE calc_psi_avg(cxxSurfaceCharge * charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std::vector<LDBLE> &zcorr);
 	LDBLE g_function(LDBLE x_value);
 	LDBLE midpnt(LDBLE x1, LDBLE x2, int n);
 	void polint(LDBLE* xa, LDBLE* ya, int n, LDBLE xv, LDBLE* yv,
 		LDBLE* dy);
 	LDBLE qromb_midpnt(cxxSurfaceCharge* charge_ptr, LDBLE x1, LDBLE x2);
-	LDBLE calc_psi_avg(cxxSurfaceCharge* charge_ptr, LDBLE surf_chrg_eq);
 
 	// inverse.cpp -------------------------------
 	int inverse_models(void);
@@ -1851,7 +1851,7 @@ isfinite handling
 #elif defined(HAVE_FINITE)
 #  define PHR_ISFINITE(x) finite(x)
 #elif defined(HAVE_ISNAN)
-#  define PHR_ISFINITE(x) ( ((x) == 0.0) || ((!isnan(x)) && ((x) != (2.0 * (x)))) )
+#  define PHR_ISFINITE(x) ( ((x) == 0.0) || ((!std::isnan(x)) && ((x) != (2.0 * (x)))) )
 #else
 #  define PHR_ISFINITE(x) ( ((x) == 0.0) || (((x) == (x)) && ((x) != (2.0 * (x)))) )
 #endif

src/phreeqcpp/Solution.cxx

@@ -47,6 +47,7 @@ cxxSolution::cxxSolution(PHRQ_io * io)
 	this->total_o = 55.55;
 	this->cb = 0.0;
 	this->density = 1.0;
+	this->viscosity = 1.0;
 	this->mass_water = 1.0;
 	this->soln_vol = 1.0;
 	this->total_alkalinity = 0.0;
@@ -80,6 +81,7 @@ cxxSolution::operator =(const cxxSolution &rhs)
 		this->total_h                    = rhs.total_h;
 		this->total_o                    = rhs.total_o;
 		this->density                    = rhs.density;
+		this->viscosity = rhs.viscosity;
 		this->cb                         = rhs.cb;
 		this->mass_water                 = rhs.mass_water;
 		this->soln_vol                   = rhs.soln_vol;
@@ -269,6 +271,10 @@ cxxSolution::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) con
 	s_oss << indent1;
 	s_oss << "-density                   " << this->density << "\n";
 
+	// new identifier
+	s_oss << indent1;
+	s_oss << "-viscosity                 " << this->viscosity << "\n";
+
 	// soln_total conc structures
 	s_oss << indent1;
 	s_oss << "-totals" << "\n";
@@ -1070,6 +1076,16 @@ cxxSolution::read_raw(CParser & parser, bool check)
 			opt_save = 27;
 		}
 		break;
+		case 28:				// viscosity
+			if (!(parser.get_iss() >> this->viscosity))
+			{
+				this->viscosity = 1.0;
+				parser.incr_input_error();
+				parser.error_msg("Expected numeric value for viscosity.",
+					PHRQ_io::OT_CONTINUE);
+			}
+			opt_save = CParser::OPT_DEFAULT;
+			break;
 		}
 		if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD)
 			break;
@@ -1365,6 +1381,7 @@ cxxSolution::zero()
 	this->total_o = 0.0;
 	this->cb = 0.0;
 	this->density = 1.0;
+	this->viscosity = 1.0;
 	this->mass_water = 0.0;
 	this->soln_vol = 0.0;
 	this->total_alkalinity = 0.0;
@@ -1397,6 +1414,7 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive)
 	this->total_o += addee.total_o * extensive;
 	this->cb += addee.cb * extensive;
 	this->density = f1 * this->density + f2 * addee.density;
+	this->viscosity = f1 * this->viscosity + f2 * addee.viscosity;
 	this->patm = f1 * this->patm + f2 * addee.patm;
 	// this->potV = f1 * this->potV + f2 * addee.potV; // appt
 	this->mass_water += addee.mass_water * extensive;
@@ -1574,6 +1592,7 @@ cxxSolution::Serialize(Dictionary & dictionary, std::vector < int >&ints,
 	doubles.push_back(this->cb);
 	doubles.push_back(this->mass_water);
 	doubles.push_back(this->density);
+	doubles.push_back(this->viscosity);
 	doubles.push_back(this->soln_vol);
 	doubles.push_back(this->total_alkalinity);
 /*
@@ -1660,6 +1679,7 @@ cxxSolution::Deserialize(Dictionary & dictionary, std::vector < int >&ints, std:
 	this->cb = doubles[dd++];
 	this->mass_water = doubles[dd++];
 	this->density = doubles[dd++];
+	this->viscosity = doubles[dd++];
 	this->soln_vol = doubles[dd++];
 	this->total_alkalinity = doubles[dd++];
 /*
@@ -1752,6 +1772,7 @@ const std::vector< std::string >::value_type temp_vopts[] = {
 	std::vector< std::string >::value_type("species_map"), 	                        // 24
 	std::vector< std::string >::value_type("log_gamma_map"), 	                    // 25
 	std::vector< std::string >::value_type("potential"), 	                        // 26
-	std::vector< std::string >::value_type("log_molalities_map")                    // 27
+	std::vector< std::string >::value_type("log_molalities_map"),                   // 27
+	std::vector< std::string >::value_type("viscosity")                             // 28
 };									   
 const std::vector< std::string > cxxSolution::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]);	

src/phreeqcpp/Solution.h

@@ -49,6 +49,8 @@ class cxxSolution:public cxxNumKeyword
 	void Set_cb(LDBLE l_cb)                             {this->cb = l_cb;}
 	LDBLE Get_density() const                           {return this->density;}
 	void Set_density(LDBLE l_density)                   {this->density = l_density;}
+	LDBLE Get_viscosity() const { return this->viscosity; }
+	void Set_viscosity(LDBLE l_viscos) { this->viscosity = l_viscos; }
 	LDBLE Get_mass_water() const                        {return this->mass_water;}
 	void Set_mass_water(LDBLE l_mass_water)             {this->mass_water = l_mass_water;}
 	LDBLE Get_total_alkalinity() const                  {return this->total_alkalinity;}
@@ -131,6 +133,7 @@ class cxxSolution:public cxxNumKeyword
 	LDBLE cb;
 	LDBLE mass_water;
 	LDBLE density;
+	LDBLE viscosity;
 	LDBLE soln_vol;
 	LDBLE total_alkalinity;
 	cxxNameDouble totals;

src/phreeqcpp/SolutionIsotope.cxx

@@ -70,7 +70,8 @@ cxxSolutionIsotope::dump_xml(std::ostream & s_oss, unsigned int indent) const
 #ifdef NPP
 	if (!isnan(this->ratio_uncertainty))
 #else
-	if (this->ratio_uncertainty != NAN)
+	//if (this->ratio_uncertainty != NAN)
+	if (!std::isnan(this->ratio_uncertainty))
 #endif
 	{
 		s_oss << indent1;

src/phreeqcpp/Surface.cxx

@@ -36,6 +36,7 @@ cxxSurface::cxxSurface(PHRQ_io *io)
 	dl_type = NO_DL;
 	sites_units = SITES_ABSOLUTE;
 	only_counter_ions = false;
+	correct_GC = false;
 	thickness = 1e-8;
 	debye_lengths = 0.0;
 	DDL_viscosity = 1.0;
@@ -55,6 +56,7 @@ cxxNumKeyword(io)
 	dl_type = NO_DL;
 	sites_units = SITES_ABSOLUTE;
 	only_counter_ions = false;
+	correct_GC = false;
 	thickness = 1e-8;
 	debye_lengths = 0.0;
 	DDL_viscosity = 1.0;
@@ -128,6 +130,8 @@ cxxSurface::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) cons
 	s_oss << indent1;
 	s_oss << "-only_counter_ions         " << this->only_counter_ions << "\n";
 	s_oss << indent1;
+	s_oss << "-correct_GC                " << this->correct_GC << "\n";
+	s_oss << indent1;
 	s_oss << "-thickness                 " << this->thickness << "\n";
 	s_oss << indent1;
 	s_oss << "-debye_lengths             " << this->debye_lengths << "\n";
@@ -189,6 +193,7 @@ cxxSurface::read_raw(CParser & parser, bool check)
 	this->Set_tidied(true);
 
 	bool only_counter_ions_defined(false);
+	//bool correct_GC_defined(false);
 	bool thickness_defined(false);
 	bool type_defined(false);
 	bool dl_type_defined(false);
@@ -468,6 +473,17 @@ cxxSurface::read_raw(CParser & parser, bool check)
 					PHRQ_io::OT_CONTINUE);
 			}
 			break;
+		case 19:				// correct_GC
+			if (!(parser.get_iss() >> this->correct_GC))
+			{
+				this->correct_GC = false;
+				parser.incr_input_error();
+				parser.
+					error_msg("Expected boolean value for correct_GC.",
+						PHRQ_io::OT_CONTINUE);
+			}
+			//correct_GC_defined = true;
+			break;
 		}
 		if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD)
 			break;
@@ -482,6 +498,13 @@ cxxSurface::read_raw(CParser & parser, bool check)
 				error_msg("Only_counter_ions not defined for SURFACE_RAW input.",
 				PHRQ_io::OT_CONTINUE);
 		}
+		//if (correct_GC_defined == false)
+		//{
+		//	parser.incr_input_error();
+		//	parser.
+		//		error_msg("correct_GC not defined for SURFACE_RAW input.",
+		//		PHRQ_io::OT_CONTINUE);
+		//}
 		if (thickness_defined == false)
 		{
 			parser.incr_input_error();
@@ -559,6 +582,7 @@ cxxSurface::add(const cxxSurface & addee_in, LDBLE extensive)
 	if (this->surface_comps.size() == 0)
 	{
 		this->only_counter_ions = addee.only_counter_ions;
+		this->correct_GC = addee.correct_GC;
 		this->dl_type = addee.dl_type;
 		this->type = addee.type;
 		this->sites_units = addee.sites_units;
@@ -730,6 +754,7 @@ cxxSurface::Serialize(Dictionary & dictionary, std::vector < int >&ints,
 	doubles.push_back(this->debye_lengths);
 	doubles.push_back(this->DDL_viscosity);
 	doubles.push_back(this->DDL_limit);
+	ints.push_back(this->correct_GC ? 1 : 0);
 	ints.push_back(this->transport ? 1 : 0);
 	this->totals.Serialize(dictionary, ints, doubles);
 	ints.push_back(this->solution_equilibria ? 1 : 0);
@@ -776,6 +801,7 @@ cxxSurface::Deserialize(Dictionary & dictionary, std::vector < int >&ints,
 	this->debye_lengths = doubles[dd++];
 	this->DDL_viscosity = doubles[dd++];
 	this->DDL_limit = doubles[dd++];
+	this->correct_GC = (ints[ii++] != 0);
 	this->transport = (ints[ii++] != 0);
 	this->totals.Deserialize(dictionary, ints, doubles, ii, dd);
 	this->solution_equilibria = (ints[ii++] != 0);
@@ -803,6 +829,7 @@ const std::vector< std::string >::value_type temp_vopts[] = {
 	std::vector< std::string >::value_type("solution_equilibria"),	// 15
 	std::vector< std::string >::value_type("n_solution"),	        // 16
 	std::vector< std::string >::value_type("totals"), 	            // 17
-	std::vector< std::string >::value_type("tidied") 	            // 18
+	std::vector< std::string >::value_type("tidied"),	            // 18
+	std::vector< std::string >::value_type("correct_gc")	        // 19
 };
 const std::vector< std::string > cxxSurface::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]);	

src/phreeqcpp/Surface.h

@@ -67,6 +67,9 @@ public:
 	void Set_DDL_viscosity(LDBLE t) {DDL_viscosity = t;}
 	LDBLE Get_DDL_limit(void) const {return DDL_limit;}
 	void Set_DDL_limit(LDBLE t) {DDL_limit = t;}
+	bool Get_correct_GC(void) const { return correct_GC; }
+	void Set_correct_GC(bool tf) { correct_GC = tf; }
+	std::vector<LDBLE> Donnan_factors;
 	bool Get_transport(void) const {return transport;}
 	void Set_transport(bool tf) {transport = tf;}
 	cxxNameDouble & Get_totals() {return this->totals;}
@@ -91,6 +94,7 @@ protected:
 	LDBLE debye_lengths;
 	LDBLE DDL_viscosity;
 	LDBLE DDL_limit;
+	bool correct_GC;
 	bool transport;
 	cxxNameDouble totals;
 	bool solution_equilibria;

src/phreeqcpp/basicsubs.cpp

@@ -187,21 +187,60 @@ diff_c(const char *species_name)
 /* ---------------------------------------------------------------------- */
 {
 	class species *s_ptr;
-	LDBLE g;
+	LDBLE ka, l_z, Dw, ff, sqrt_mu;
+	sqrt_mu = sqrt(mu_x);
 
 	s_ptr = s_search(species_name);
-	if (s_ptr != NULL /*&& s_ptr->in != FALSE && s_ptr->type < EMINUS*/)
+	//LDBLE g;
+	//if (s_ptr != NULL /*&& s_ptr->in != FALSE && s_ptr->type < EMINUS*/)
+	//{
+	//	g = s_ptr->dw;
+	//	if (s_ptr->dw_t)
+	//			g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
+	//	g *= viscos_0_25 / viscos * tk_x / 298.15;
+	//}
+	//else
+	//{
+	//	g = 0;
+	//}
+	//return (g);
+	if (s_ptr == NULL)
+		return(0);
+	if ((Dw = s_ptr->dw) == 0)
 	{
-		g = s_ptr->dw;
+		if (correct_Dw)
-		if (s_ptr->dw_t)
+			Dw = default_Dw;
-				g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
+		else
-		g *= viscos_0_25 / viscos * tk_x / 298.15;
+			return(0);
+	}
+	if ((l_z = fabs(s_ptr->z)) == 0)
+	{
+		//l_z = 1; // only a 1st approximation for correct_Dw in electrical field
 	}
 	else
 	{
-		g = 0;
+		if (s_ptr->dw_a2)
+			ka = DH_B * s_ptr->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75));
+		else
+			ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75));
+		if (s_ptr->dw_a)
+		{
+			ff = exp(-s_ptr->dw_a * DH_A * l_z * sqrt_mu / (1 + ka));
+			//if (print_viscosity && s_ptr->dw_a_visc)
+			//	ff *= pow((viscos_0 / viscos), s_ptr->dw_a_visc);
+		}
+		else
+		{
+			ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka));
+		}
+		Dw *= ff;
 	}
-	return (g);
+
+	if (tk_x != 298.15 && s_ptr->dw_t)
+		Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
+
+	s_ptr->dw_corr = Dw;
+	return (Dw * viscos_0_25 / viscos_0);
 }
 /* ---------------------------------------------------------------------- */
 LDBLE Phreeqc::
@@ -209,22 +248,58 @@ setdiff_c(const char *species_name, double d)
 /* ---------------------------------------------------------------------- */
 {
 	class species *s_ptr;
-	LDBLE g;
+	LDBLE ka, l_z, Dw, ff, sqrt_mu;
+	sqrt_mu = sqrt(mu_x);
 
 	s_ptr = s_search(species_name);
-	if (s_ptr != NULL)
+
+	//LDBLE g;
+	//s_ptr = s_search(species_name);
+	//if (s_ptr != NULL)
+	//{
+
+	//	s_ptr->dw = d;
+	//	g = s_ptr->dw;
+	//	if (s_ptr->dw_t)
+	//			g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
+	//	g *= viscos_0_25 / viscos * tk_x / 298.15;;
+	//}
+	//else
+	//{
+	//	g = 0;
+	//}
+	//return (g);
+	if (s_ptr == NULL)
+		return(0);
+	Dw = s_ptr->dw = d;
+	if ((l_z = fabs(s_ptr->z)) == 0)
 	{
-		s_ptr->dw = d;
+		//l_z = 1; // only a 1st approximation for correct_Dw in electrical field
-		g = s_ptr->dw;
-		if (s_ptr->dw_t)
-				g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
-		g *= viscos_0_25 / viscos * tk_x / 298.15;;
 	}
 	else
 	{
-		g = 0;
+		if (s_ptr->dw_a2)
+			ka = DH_B * s_ptr->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75));
+		else
+			ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75));
+		if (s_ptr->dw_a)
+		{
+			ff = exp(-s_ptr->dw_a * DH_A * l_z * sqrt_mu / (1 + ka));
+			//if (print_viscosity && s_ptr->dw_a_visc)
+			//	ff *= pow((viscos_0 / viscos), s_ptr->dw_a_visc);
+		}
+		else
+		{
+			ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka));
+		}
+		Dw *= ff;
 	}
-	return (g);
+
+	if (tk_x != 298.15 && s_ptr->dw_t)
+			Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
+
+	s_ptr->dw_corr = Dw;
+	return (Dw * viscos_0_25 / viscos_0);
 }
 /* ---------------------------------------------------------------------- */
 LDBLE Phreeqc::
@@ -314,19 +389,18 @@ calc_SC(void)
 
 	SC = 0;
 	//LDBLE ta1, ta2, ta3, ta4;
-	for (i = 0; i < (int)this->s_x.size(); i++)
+	//for (i = 0; i < (int)this->s_x.size(); i++)
-	{
+	//{
-		// ** for optimizing, get numbers from -analyt for H+ = H+...
+	//	// ** for optimizing, get numbers from -analyt for H+ = H+...
-		//if (!strcmp(s_x[i]->name, "H+"))
+	//	if (!strcmp(s_x[i]->name, "H+"))
-		//{
+	//	{
-		//	ta1 = s_x[i]->logk[2];
+	//		ta1 = s_x[i]->logk[2] * 1e15;
-		//	ta2 = s_x[i]->logk[3];
+	//		ta2 = s_x[i]->logk[3] * 1e15;
-		//	ta3 = s_x[i]->logk[4];
+	//		ta3 = s_x[i]->logk[4] * 1e15;
-		//	ta4 = s_x[i]->logk[5];
+	//		ta4 = s_x[i]->logk[5] * 1e15;
-		//	break;
+	//		break;
-		//}
+	//	}
-		//
+	//}
-	}
 	for (i = 0; i < (int)this->s_x.size(); i++)
 	{
 		if (s_x[i]->type != AQ && s_x[i]->type != HPLUS)
@@ -338,43 +412,51 @@ calc_SC(void)
 			else
 				continue;
 		}
+		if (s_x[i]->lm < min_dif_LM)
+			continue;
 		if ((l_z = fabs(s_x[i]->z)) == 0)
-			l_z = 1; // only a 1st approximation for correct_Dw in electrical field
+		{
-		if (s_x[i]->dw_t)
+			//l_z = 1; // only a 1st approximation for correct_Dw in electrical field
-			Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); // the viscosity multiplier is done in SC
+		}
-		if (s_x[i]->dw_a2)
-			ka = DH_B * s_x[i]->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75));
 		else
 		{
-			ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x , 0.75));
+			if (s_x[i]->dw_a2)
-			//ka = DH_B * ta1 * sqrt_mu / (1 + pow(mu_x, ta2));
+				ka = DH_B * s_x[i]->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75));
-			//ka = DH_B * ta1 * sqrt_mu / (1 + mu_x / ta2);
+			else
-		}
-		if (s_x[i]->dw_a)
-		{
-			ff = exp(-s_x[i]->dw_a * DH_A * l_z * sqrt_mu / (1 + ka));
-			if (print_viscosity)
 			{
-				ff *= pow((viscos_0 / viscos), s_x[i]->dw_a_visc);
+				ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75));
+				//ka = DH_B * ta1 * sqrt_mu / (1 + pow(mu_x, ta2));
+				//ka = DH_B * ta1 * sqrt_mu / (1 + mu_x / ta2);
 			}
+			if (s_x[i]->dw_a)
+			{
+				ff = exp(-s_x[i]->dw_a * DH_A * l_z * sqrt_mu / (1 + ka));
+				//if (print_viscosity && s_x[i]->dw_a_visc)
+				//	ff *= pow((viscos_0 / viscos), s_x[i]->dw_a_visc);
+			}
+			else
+			{
+				ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka));
+				//ff = exp(-ta3 * DH_A * l_z * sqrt_mu / (1 + ka));
+			}
+			Dw *= ff;
 		}
-		else
+		if (tk_x != 298.15)
 		{
-			ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka));
+			if (s_x[i]->dw_t)
-			//ff = exp(-ta3 * DH_A * l_z * sqrt_mu / (1 + ka));
+				Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15);
+			//else
+			//{
+			//	Dw *= exp(ta1 / tk_x - ta1 / 298.15);
+			//}
 		}
-		Dw *= ff;
 		s_x[i]->dw_corr = Dw;
 
-		if (s_x[i]->z == 0)
-			continue;
 		s_x[i]->dw_t_SC = s_x[i]->moles / mass_water_aq_x * l_z * l_z * Dw;
 		SC += s_x[i]->dw_t_SC;
 	}
 	SC *= 1e7 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298150.0);
-	/* correct for temperature dependency...
+	/* correct for viscosity dependency...
-	SC_T = SC_298 * (Dw_T / T) * (298 / Dw_298) and
-	Dw_T = Dw_298 * (T / 298) * (viscos_298 / viscos_T) give:
 	SC_T = SC_298 * (viscos_298 / viscos_T)
 	*/
 	SC *= viscos_0_25 / viscos_0;
@@ -1520,7 +1602,8 @@ get_calculate_value(const char *name)
 #ifdef NPP
 	if (isnan(rate_moles))
 #else
-	if (rate_moles == NAN)
+	//if (rate_moles == NAN)
+	if(std::isnan(rate_moles))
 #endif
 	{
 		error_string = sformatf( "Calculated value not SAVEed for %s.",
@@ -1721,22 +1804,22 @@ pressure(void)
 
 /* ---------------------------------------------------------------------- */
 LDBLE Phreeqc::
-pr_phi(const char* phase_name)
+pr_phi(const char *phase_name)
 /* ---------------------------------------------------------------------- */
 {
 	cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
 	if (gas_phase_ptr != NULL)
 	{
-		int l;
+	int l;
 		class phase* phase_ptr = phase_bsearch(phase_name, &l, FALSE);
-		if (phase_ptr == NULL)
+	if (phase_ptr == NULL)
-		{
+	{
-			error_string = sformatf("Gas %s, not found.", phase_name);
+		error_string = sformatf( "Gas %s, not found.", phase_name);
-			warning_msg(error_string);
+		warning_msg(error_string);
-			return (1e-99);
+		return (1e-99);
-		}
+	}
 		for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
-		{
+	{
 			const cxxGasComp* gas_comp_ptr = &(gas_phase_ptr->Get_gas_comps()[i]);
 			int j;
 			class phase* phase_ptr_gas = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE);
@@ -1744,8 +1827,8 @@ pr_phi(const char* phase_name)
 			{
 				if (gas_phase_ptr->Get_pr_in())
 				{
-					return phase_ptr->pr_phi;
+		return phase_ptr->pr_phi;
-				}
+	}
 				else
 				{
 					return gas_comp_ptr->Get_phi();
@@ -1753,7 +1836,7 @@ pr_phi(const char* phase_name)
 			}
 		}
 	}
-	return(1.0);
+	return (1.0);
 }
 /* ---------------------------------------------------------------------- */
 LDBLE Phreeqc::

src/phreeqcpp/cl1.cpp

@@ -1,6 +1,7 @@
 #include <stdio.h>
 #include <string.h>
-#include <math.h>
+//#include <math.h>
+#include <cmath>
 #include <stdlib.h>
 #include "Phreeqc.h"
 #include "phqalloc.h"

src/phreeqcpp/class_main.cpp

@@ -323,7 +323,7 @@ write_banner(void)
 
 	/* date */
 #ifdef NPP
-	len = snprintf(buffer, sizeof(buffer), "%s", "July 5, 2021");
+	len = snprintf(buffer, sizeof(buffer), "%s", "March 20, 2023");
 #else
 	len = snprintf(buffer, sizeof(buffer), "%s", "@VER_DATE@");
 #endif
@@ -507,7 +507,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
 		output_msg(sformatf("   Input file: %s\n", in_file.c_str()));
 		output_msg(sformatf("  Output file: %s\n", out_file.c_str()));
 #ifdef NPP
-		output_msg(sformatf("Using PHREEQC: version 3.7.1, compiled July 5, 2021\n"));
+		output_msg(sformatf("Using PHREEQC: version 3.7.3, compiled March 20, 2023\n"));
 #endif
 		output_msg(sformatf("Database file: %s\n\n", token.c_str()));
 #ifdef NPP

src/phreeqcpp/common/Utils.cxx

@@ -10,7 +10,8 @@
 #include "Utils.h"
 #include "Parser.h"
 #include "float.h"
-#include "math.h"
+//#include <math.h>
+#include <cmath>
 
 #if defined(PHREEQCI_GUI)
 #ifdef _DEBUG

src/phreeqcpp/gases.cpp

@@ -624,6 +624,7 @@ int Phreeqc::
 calc_fixed_volume_gas_pressures(void)
 /* ---------------------------------------------------------------------- */
 {
+	//int n_g = 0;
 	LDBLE lp;
 	class rxn_token *rxn_ptr;
 	class phase *phase_ptr;
@@ -644,6 +645,7 @@ calc_fixed_volume_gas_pressures(void)
 		{
 			if (!PR && phase_ptr->t_c > 0 && phase_ptr->p_c > 0)
 				PR = true;
+			//n_g++;
 		}
 		gas_phase_ptr->Set_total_moles(gas_phase_ptr->Get_total_moles() + gas_unknowns[i]->moles);
 	}

src/phreeqcpp/global_structures.h

@@ -5,8 +5,8 @@
 /* ----------------------------------------------------------------------
  *   #define DEFINITIONS
  * ---------------------------------------------------------------------- */
-#ifndef NAN
+#if !defined(NAN)
-#   define NAN -99999999
+#define NAN nan("1")
 #endif
 #define MISSING -9999.999            
 #include "NA.h"   /* NA = not available */

src/phreeqcpp/integrate.cpp

@@ -25,9 +25,9 @@ calc_all_g(void)
 
 	if (use.Get_surface_ptr() == NULL)
 		return (OK);
-/*
+	/*
- *   calculate g for each surface
+	 *   calculate g for each surface
- */
+	 */
 	epsilon = convergence_tolerance;
 	if (convergence_tolerance >= 1e-8)
 	{
@@ -45,7 +45,7 @@ calc_all_g(void)
 		if (x[j]->type != SURFACE_CB)
 			continue;
 		if (debug_diffuse_layer == TRUE)
-			output_msg(sformatf( "Calc_all_g, X[%d]\n", j));
+			output_msg(sformatf("Calc_all_g, X[%d]\n", j));
 		cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
 		std::map<LDBLE, cxxSurfDL> temp_g_map;
 		cxxSurfDL temp_g;
@@ -56,11 +56,11 @@ calc_all_g(void)
 		/* 1000 J/kJ and 1000 L/m**3 */
 		//alpha_global =	sqrt(EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 *
 		//		 tk_x * 0.5);
-		alpha_global =	sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 *
+		alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 *
-				 tk_x * 0.5);
+			tk_x * 0.5);
-/*
+		/*
- *   calculate g for given surface for each species
+		 *   calculate g for given surface for each species
- */
+		 */
 		for (int i = 0; i < (int)this->s_x.size(); i++)
 		{
 			if (s_x[i]->type > HPLUS)
@@ -73,7 +73,7 @@ calc_all_g(void)
 
 				if ((use.Get_surface_ptr()->Get_only_counter_ions() == false) ||
 					(((x[j]->master[0]->s->la > 0) && (z_global < 0))
-					 || ((x[j]->master[0]->s->la < 0) && (z_global > 0))))
+						|| ((x[j]->master[0]->s->la < 0) && (z_global > 0))))
 				{
 					if (xd_global > 0.1)
 					{
@@ -180,11 +180,11 @@ calc_all_g(void)
 				if (debug_diffuse_layer == TRUE)
 				{
 					output_msg(sformatf(
-							   "\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
+						"\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
-							   (double) z_global,
+						(double)z_global,
-							   (double) charge_ptr->Get_g_map()[z_global].Get_g(),
+						(double)charge_ptr->Get_g_map()[z_global].Get_g(),
-							   (double) new_g,
+						(double)new_g,
-							   (double) (new_g - charge_ptr->Get_g_map()[z_global].Get_g())));
+						(double)(new_g - charge_ptr->Get_g_map()[z_global].Get_g())));
 				}
 			}
 			charge_ptr->Get_g_map()[z_global].Set_g(new_g);
@@ -201,7 +201,7 @@ calc_all_g(void)
 						g_function(xd_global) / F_C_MOL;
 					dg *=
 						-2. / (exp(x[j]->master[0]->s->la * LOG_10) *
-							   exp(x[j]->master[0]->s->la * LOG_10));
+							exp(x[j]->master[0]->s->la * LOG_10));
 					if ((xd_global - 1) < 0.0)
 					{
 						dg *= -1.0;
@@ -225,16 +225,16 @@ calc_all_g(void)
 		if (debug_diffuse_layer == TRUE)
 		{
 			output_msg(sformatf("\nSurface component %d: charge,\tg,\tdg/dlny,\txd\n",
-					  (int) charge_ptr->Get_g_map().size()));
+				(int)charge_ptr->Get_g_map().size()));
 			std::map<LDBLE, cxxSurfDL>::iterator it;
 			for (it = charge_ptr->Get_g_map().begin(); it != charge_ptr->Get_g_map().end(); it++)
 			{
 				output_msg(sformatf(
-						   "\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
+					"\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
-						   (double) it->first,
+					(double)it->first,
-						   (double) it->second.Get_g(),
+					(double)it->second.Get_g(),
-						   (double) it->second.Get_dg(),
+					(double)it->second.Get_dg(),
-						   (double) xd_global));
+					(double)xd_global));
 			}
 		}
 	}
@@ -256,7 +256,7 @@ g_function(LDBLE x_value)
 
 	cxxSurfaceCharge *charge_ptr = &(use.Get_surface_ptr()->Get_surface_charges()[0]);
 	std::map<LDBLE, cxxSurfDL>::iterator it = charge_ptr->Get_g_map().begin();
-	for ( ; it != charge_ptr->Get_g_map().end(); it++)
+	for (; it != charge_ptr->Get_g_map().end(); it++)
 	{
 		it->second.Set_psi_to_z(exp(ln_x_value * it->first) - 1.0);
 	}
@@ -272,31 +272,31 @@ g_function(LDBLE x_value)
 		sum = 0.0;
 		sum1 = 0.0;
 		output_msg(sformatf(
-				   "Species\tmoles\tX**z-1\tsum\tsum charge\n"));
+			"Species\tmoles\tX**z-1\tsum\tsum charge\n"));
 		for (i = 0; i < (int)this->s_x.size(); i++)
 		{
 			if (s_x[i]->type < H2O && s_x[i]->z != 0.0)
 			{
 				sum += s_x[i]->moles * (pow(x_value, s_x[i]->z) - 1.0);
 				sum1 += s_x[i]->moles * s_x[i]->z;
-				output_msg(sformatf( "%s\t%e\t%e\t%e\t%e\n",
+				output_msg(sformatf("%s\t%e\t%e\t%e\t%e\n",
-						   s_x[i]->name, (double) s_x[i]->moles,
+					s_x[i]->name, (double)s_x[i]->moles,
-						   (double) (pow((LDBLE) x_value, (LDBLE) s_x[i]->z) -
+					(double)(pow((LDBLE)x_value, (LDBLE)s_x[i]->z) -
-									 1.0), (double) sum, (double) sum1));
+						1.0), (double)sum, (double)sum1));
 			}
 		}
-		error_string = sformatf( "Negative sum in g_function, %e\t%e.",
+		error_string = sformatf("Negative sum in g_function, %e\t%e.",
-				(double) sum, (double) x_value);
+			(double)sum, (double)x_value);
 		error_msg(error_string, CONTINUE);
 		error_string = sformatf(
-				"Solutions must be charge balanced, charge imbalance is %e\n",
+			"Solutions must be charge balanced, charge imbalance is %e\n",
-				(double) sum1);
+			(double)sum1);
 		error_msg(error_string, STOP);
 	}
 
 	return_value =
 		(exp(ln_x_value * z_global) -
-		 1) / sqrt((x_value * x_value * mass_water_aq_x * sum));
+			1) / sqrt((x_value * x_value * mass_water_aq_x * sum));
 	return (return_value);
 }
 /* ---------------------------------------------------------------------- */
@@ -309,9 +309,9 @@ polint(LDBLE * xa, LDBLE * ya, int n, LDBLE xv, LDBLE * yv, LDBLE * dy)
 
 	ns = 1;
 	dif = fabs(xv - xa[1]);
-/*
+	/*
- *   Malloc work space
+	 *   Malloc work space
- */
+	 */
 	std::vector<double> c, d;
 	c.resize((size_t)n + 1);
 	d.resize((size_t)n + 1);
@@ -354,7 +354,7 @@ polint(LDBLE * xa, LDBLE * ya, int n, LDBLE xv, LDBLE * yv, LDBLE * dy)
 		}
 		*yv += *dy;
 
-/*		*yv += (*dy = (2 * ns < (n-m) ? c[ns+1] : d[ns--])); */
+		/*		*yv += (*dy = (2 * ns < (n-m) ? c[ns+1] : d[ns--])); */
 	}
 	return;
 }
@@ -376,7 +376,7 @@ midpnt(LDBLE x1, LDBLE x2, int n)
 	{
 		for (it = 1, j = 1; j < n - 1; j++)
 			it *= 3;
-		tnm = (LDBLE) it;
+		tnm = (LDBLE)it;
 		del = (x2 - x1) / (3 * tnm);
 		ddel = del + del;
 		xv = x1 + 0.5 * del;
@@ -420,7 +420,7 @@ qromb_midpnt(cxxSurfaceCharge *charge_ptr, LDBLE x1, LDBLE x2)
 			if (debug_diffuse_layer == TRUE)
 			{
 				output_msg(sformatf(
-						   "Iterations in qromb_midpnt: %d\n", j));
+					"Iterations in qromb_midpnt: %d\n", j));
 			}
 			return (sv[j]);
 		}
@@ -436,7 +436,7 @@ qromb_midpnt(cxxSurfaceCharge *charge_ptr, LDBLE x1, LDBLE x2)
 				if (debug_diffuse_layer == TRUE)
 				{
 					output_msg(sformatf(
-							   "Iterations in qromb_midpnt: %d\n", j));
+						"Iterations in qromb_midpnt: %d\n", j));
 				}
 				return (ss);
 			}
@@ -444,7 +444,7 @@ qromb_midpnt(cxxSurfaceCharge *charge_ptr, LDBLE x1, LDBLE x2)
 
 	}
 	error_string = sformatf(
-			"\nToo many iterations integrating diffuse layer.\n");
+		"\nToo many iterations integrating diffuse layer.\n");
 	error_msg(error_string, STOP);
 	return (-999.9);
 }
@@ -455,9 +455,9 @@ calc_init_g(void)
 {
 	if (use.Get_surface_ptr() == NULL)
 		return (OK);
-/*
+	/*
- *   calculate g for each surface
+	 *   calculate g for each surface
- */
+	 */
 	for (int j = 0; j < count_unknowns; j++)
 	{
 		if (x[j]->type != SURFACE_CB)
@@ -468,7 +468,7 @@ calc_init_g(void)
 		/*  second 1000 is liters/m**3 */
 		//alpha_global =	sqrt(EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
 		//	1000.0 * tk_x * 0.5);
-		alpha_global =	sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
+		alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
 			1000.0 * tk_x * 0.5);
 
 		if (charge_ptr->Get_g_map().size() == 0)
@@ -476,9 +476,9 @@ calc_init_g(void)
 			cxxSurfDL temp_g;
 			charge_ptr->Get_g_map()[0.0] = temp_g;
 		}
-/*
+		/*
- *   calculate g for given surface for each species
+		 *   calculate g for given surface for each species
- */
+		 */
 		for (int i = 0; i < (int)this->s_x.size(); i++)
 		{
 			if (s_x[i]->type > HPLUS)
@@ -511,7 +511,7 @@ calc_init_g(void)
 
 			{
 				int is = s_x[i]->number;
-				assert (is < (int) s_diff_layer.size());
+				assert(is < (int)s_diff_layer.size());
 				// species found in diff_layer
 				s_diff_layer[is][charge_ptr->Get_name()].Set_g_moles(0);
 				s_diff_layer[is][charge_ptr->Get_name()].Set_dg_g_moles(0);
@@ -520,15 +520,15 @@ calc_init_g(void)
 		if (debug_diffuse_layer == TRUE)
 		{
 			output_msg(sformatf(
-					   "\nSurface component %d: charge,\tg,\tdg\n",
+				"\nSurface component %d: charge,\tg,\tdg\n",
-					   (int) charge_ptr->Get_g_map().size()));
+				(int)charge_ptr->Get_g_map().size()));
 			std::map<LDBLE, cxxSurfDL>::iterator it;
 			for (it = charge_ptr->Get_g_map().begin(); it != charge_ptr->Get_g_map().end(); it++)
 			{
-				output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n",
+				output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\n",
-						   (double) it->first,
+					(double)it->first,
-						   (double) it->second.Get_g(),
+					(double)it->second.Get_g(),
-						   (double) it->second.Get_dg()));
+					(double)it->second.Get_dg()));
 			}
 		}
 	}
@@ -539,21 +539,21 @@ int Phreeqc::
 initial_surface_water(void)
 /* ---------------------------------------------------------------------- */
 {
-/*
+	/*
- *   In initial surface calculation, need to calculate
+	 *   In initial surface calculation, need to calculate
- *   mass of water in diffuse layer.
+	 *   mass of water in diffuse layer.
- *   diffuse layer water + aqueous solution water = bulk water.
+	 *   diffuse layer water + aqueous solution water = bulk water.
- *   Ionic strength is fixed, so diffuse-layer water will not change
+	 *   Ionic strength is fixed, so diffuse-layer water will not change
- */
+	 */
 	LDBLE debye_length, b, r, rd, ddl_limit, rd_limit, fraction, sum_surfs, l_s;
 	LDBLE damp_aq;
-/*
+	/*
- *   Debye  length = 1/k = sqrt[eta*eta_zero*R*T/(2*F**2*mu_x*1000)], Dzombak and Morel, p 36
+	 *   Debye  length = 1/k = sqrt[eta*eta_zero*R*T/(2*F**2*mu_x*1000)], Dzombak and Morel, p 36
- *
+	 *
- *   1000 converts kJ to J; 1000 converts Liters to meter**3; debye_length is in meters.
+	 *   1000 converts kJ to J; 1000 converts Liters to meter**3; debye_length is in meters.
- */
+	 */
-	//debye_length = (EPSILON * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x)
+	 //debye_length = (EPSILON * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x)
-	//	/ (2. * F_C_MOL * F_C_MOL * mu_x * 1000.);
+	 //	/ (2. * F_C_MOL * F_C_MOL * mu_x * 1000.);
 	debye_length = (eps_r * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x)
 		/ (2. * F_C_MOL * F_C_MOL * mu_x * 1000.);
 	debye_length = sqrt(debye_length);
@@ -561,10 +561,10 @@ initial_surface_water(void)
 	/*  ddl is at most the fraction ddl_limit of bulk water */
 	ddl_limit = use.Get_surface_ptr()->Get_DDL_limit();
 
-/*
+	/*
- *   Loop through all surface components, calculate each H2O surface (diffuse layer),
+	 *   Loop through all surface components, calculate each H2O surface (diffuse layer),
- *   H2O aq, and H2O bulk (diffuse layers plus aqueous).
+	 *   H2O aq, and H2O bulk (diffuse layers plus aqueous).
- */
+	 */
 
 	if (use.Get_surface_ptr()->Get_debye_lengths() > 0)
 	{
@@ -606,14 +606,14 @@ initial_surface_water(void)
 				mass_water_surfaces_x =
 					use.Get_solution_ptr()->Get_mass_water() * ddl_limit / (1 - ddl_limit);
 				r = 0.002 * (mass_water_surfaces_x +
-							 use.Get_solution_ptr()->Get_mass_water()) / sum_surfs;
+					use.Get_solution_ptr()->Get_mass_water()) / sum_surfs;
 				rd_limit = (1 - sqrt(1 - ddl_limit)) * r;
 				rd = rd_limit;
 				use.Get_surface_ptr()->Set_thickness(rd);
 			}
 			else
 				mass_water_surfaces_x =
-					(r * r / pow(r - rd, 2) - 1) * use.Get_solution_ptr()->Get_mass_water();
+				(r * r / pow(r - rd, 2) - 1) * use.Get_solution_ptr()->Get_mass_water();
 			for (int i = 0; i < count_unknowns; i++)
 			{
 				if (x[i]->type != SURFACE_CB)
@@ -692,14 +692,14 @@ sum_diffuse_layer(cxxSurfaceCharge *charge_ptr)
 
 	if (use.Get_surface_ptr() == NULL)
 		return (OK);
-/*
+	/*
- *   Find position of component in list of components
+	 *   Find position of component in list of components
- */
+	 */
 
-/*
+	 /*
- *   Loop through all surface components, calculate each H2O surface (diffuse layer),
+	  *   Loop through all surface components, calculate each H2O surface (diffuse layer),
- *   H2O aq, and H2O bulk (diffuse layers plus aqueous).
+	  *   H2O aq, and H2O bulk (diffuse layers plus aqueous).
- */
+	  */
 	count_elts = 0;
 	paren_count = 0;
 	mass_water_surface = charge_ptr->Get_mass_water();
@@ -717,9 +717,9 @@ sum_diffuse_layer(cxxSurfaceCharge *charge_ptr)
 		}
 		moles_excess = mass_water_aq_x * molality * g;
 		moles_surface = mass_water_surface * molality + moles_excess;
-/*
+		/*
- *   Accumulate elements in diffuse layer
+		 *   Accumulate elements in diffuse layer
- */
+		 */
 		add_elt_list(s_x[j]->next_elt, moles_surface);
 	}
 	add_elt_list(s_h2o->next_elt, mass_water_surface / gfw_water);
@@ -733,20 +733,31 @@ calc_all_donnan(void)
 {
 	bool converge;
 	int cd_m;
-	LDBLE new_g, f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq, ratio_aq_tot;
+	LDBLE new_g, f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq, ratio_aq_tot, co_ion;
 	LDBLE new_g2, f_psi2, surf_chrg_eq2, psi_avg2, dif, var1;
 
 	if (use.Get_surface_ptr() == NULL)
 		return (OK);
-	//f_sinh = sqrt(8000.0 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
-	//		 tk_x * mu_x);
 	f_sinh = sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
-			 tk_x * mu_x);
+		tk_x * mu_x);
-/*
+	bool only_count = use.Get_surface_ptr()->Get_only_counter_ions();
- *   calculate g for each surface...
+	bool correct_GC = use.Get_surface_ptr()->Get_correct_GC();
- */
+	/*   calculate g for each surface...
-	if (!calculating_deriv || use.Get_surface_ptr()->Get_debye_lengths()) // DL_pitz
+	 */
+	if (!calculating_deriv || use.Get_surface_ptr()->Get_debye_lengths() ||
+		correct_GC) // DL_pitz && correct_GC
 		initial_surface_water();
+	LDBLE nDbl = 1;
+	if (correct_GC)
+	{
+		if ((nDbl = use.Get_surface_ptr()->Get_debye_lengths()) == 0)
+		{
+			LDBLE debye_length = f_sinh / (F_C_MOL * mu_x * 4e3);
+			nDbl = use.Get_surface_ptr()->Get_thickness() / debye_length;
+			//use.Get_surface_ptr()->Set_debye_lengths(nDbl);
+		}
+	}
+
 	converge = TRUE;
 	for (int j = 0; j < count_unknowns; j++)
 	{
@@ -755,10 +766,10 @@ calc_all_donnan(void)
 		cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
 
 		if (debug_diffuse_layer == TRUE)
-			output_msg(sformatf( "Calc_all_g, X[%d]\n", j));
+			output_msg(sformatf("Calc_all_g, X[%d]\n", j));
-/*
+		/*
- *  sum eq of each charge number in solution...
+		 *  sum eq of each charge number in solution...
- */
+		 */
 		std::map<LDBLE, LDBLE>::iterator it;
 		for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
 		{
@@ -778,15 +789,18 @@ calc_all_donnan(void)
 			f_psi = x[(size_t)j + 2]->master[0]->s->la * LOG_10;	/* -FPsi/RT */
 			f_psi = f_psi / 2;
 			cd_m = 1;
-		} else
+		}
+		else
 		{
 			f_psi = x[j]->master[0]->s->la * LOG_10;
 			cd_m = -1;
 		}
 		surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;
-		if (fabs(surf_chrg_eq) > 5e3)
+		LDBLE lim_seq = 5e3;
+		if (correct_GC) lim_seq = 5e1;
+		if (fabs(surf_chrg_eq) > lim_seq)
 		{
-			surf_chrg_eq = (surf_chrg_eq < 0 ? -5e3 : 5e3);
+			surf_chrg_eq = (surf_chrg_eq < 0 ? -lim_seq : lim_seq);
 			var1 = surf_chrg_eq / (A_surf * f_sinh / F_C_MOL);
 			var1 = (var1 + sqrt(var1 * var1 + 1));
 			f_psi = (var1 > 1e-8 ? log(var1) : -18.4);
@@ -799,8 +813,11 @@ calc_all_donnan(void)
 		surf_chrg_eq2 = A_surf * f_sinh * sinh(f_psi2) / F_C_MOL;
 
 		/* find psi_avg that matches surface charge... */
-		psi_avg = calc_psi_avg(charge_ptr, surf_chrg_eq);
+		std::vector<LDBLE> zcorr(charge_group_map.size());
-		psi_avg2 = calc_psi_avg(charge_ptr, surf_chrg_eq2);
+		std::vector<LDBLE> zcorr2(charge_group_map.size());
+		//LDBLE fD = 0;
+		psi_avg = calc_psi_avg(charge_ptr, surf_chrg_eq, nDbl, zcorr);
+		psi_avg2 = calc_psi_avg(charge_ptr, surf_chrg_eq2, nDbl, zcorr2);
 
 		/*output_msg(sformatf( "psi's  %e %e %e\n", f_psi, psi_avg, surf_chrg_eq)); */
 
@@ -808,9 +825,15 @@ calc_all_donnan(void)
 		ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
 		ratio_aq_tot = charge_ptr->Get_mass_water() / mass_water_bulk_x;
 
+		int z_iter = 0;
 		for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
 		{
-			LDBLE z = it->first;
+			LDBLE z = it->first, z1 = z;
+			co_ion = surf_chrg_eq * z;
+			if (correct_GC)
+				z1 = zcorr[z_iter];
+			//z1 *= cgc[0] * pow(z_factor, abs(z));
+
 			if (!ratio_aq)
 			{
 				charge_ptr->Get_g_map()[z].Set_g(0);
@@ -819,17 +842,13 @@ calc_all_donnan(void)
 				converge = true;
 				continue;
 			}
-			new_g = ratio_aq * (exp(cd_m * z * psi_avg) - 1);
+			new_g = ratio_aq * (exp(cd_m * z1 * psi_avg) - 1);
-			if (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0)
+			if (only_count && co_ion > 0)
-				//((surf_chrg_eq < 0 && z < 0)
-				// || (surf_chrg_eq > 0 && z > 0)))
 				new_g = -ratio_aq;
 			if (new_g <= -ratio_aq)
 				new_g = -ratio_aq + G_TOL * 1e-3;
-			new_g2 = ratio_aq * (exp(cd_m * z * psi_avg2) - 1);
+			new_g2 = ratio_aq * (exp(cd_m * z1 * psi_avg2) - 1);
-			if (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0)
+			if (only_count && co_ion > 0)
-				//((surf_chrg_eq < 0 && z < 0)
-				// || (surf_chrg_eq > 0 && z > 0)))
 				new_g2 = -ratio_aq;
 			if (new_g2 <= -ratio_aq)
 				new_g2 = -ratio_aq + G_TOL * 1e-3;
@@ -859,9 +878,9 @@ calc_all_donnan(void)
 			/* save Boltzmann factor * water fraction for MCD calc's in transport */
 			if (converge)
 			{
-				if (use.Get_surface_ptr()->Get_only_counter_ions())
+				if (only_count)
 				{
-					if (surf_chrg_eq * z > 0) // co-ions are not in the DL
+					if (co_ion > 0) // co-ions are not in the DL
 						charge_ptr->Get_z_gMCD_map()[z] = 0;
 					else // assume that counter-ions have the free water conc for diffusion
 						charge_ptr->Get_z_gMCD_map()[z] = ratio_aq_tot;
@@ -869,23 +888,27 @@ calc_all_donnan(void)
 				else
 					charge_ptr->Get_z_gMCD_map()[z] = (new_g / ratio_aq + 1) * ratio_aq_tot;
 			}
+			z_iter++;
 		}
 		if (debug_diffuse_layer == TRUE)
 		{
-			std::string name =  x[j]->master[0]->elt->name;
+			std::string name = x[j]->master[0]->elt->name;
 			Utilities::replace("_psi", "", name);
 			output_msg(sformatf(
-					   "\nDonnan all on %s (%d): charge, \tg, \tdg, Psi_surface = %8f V. \n",
+				"\nDonnan all on %s (%d): charge, \tg, \tdg, \tzcorr, \tPsi_surface = %8f V, \tDebye lengths = %8f. \n",
-					   name.c_str(), (int) charge_ptr->Get_g_map().size(),
+				name.c_str(), (int)charge_ptr->Get_g_map().size(),
-					   x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL *
+				x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ,
-					   tk_x / F_KJ_V_EQ));
+				nDbl));
+			int i = 0;
 			for (std::map<LDBLE, cxxSurfDL>::iterator i_it = charge_ptr->Get_g_map().begin();
 				i_it != charge_ptr->Get_g_map().end(); i_it++)
 			{
-				output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n",
+				output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
-						   (double) i_it->first,
+					(double)i_it->first,
-						   (double) i_it->second.Get_g(),
+					(double)i_it->second.Get_g(),
-						   (double) i_it->second.Get_dg()));
+					(double)i_it->second.Get_dg(),
+					(double)zcorr[i]));
+				i++;
 			}
 		}
 	}
@@ -904,7 +927,7 @@ calc_init_donnan(void)
 		//sqrt(8000.0 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
 		//	 tk_x * mu_x);
 		sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
-			 tk_x * mu_x);
+			tk_x * mu_x);
 	if (convergence_tolerance >= 1e-8)
 	{
 		G_TOL = 1e-9;
@@ -913,9 +936,9 @@ calc_init_donnan(void)
 	{
 		G_TOL = 1e-13;
 	}
-/*
+	/*
- *  sum eq of each charge number in solution...
+	 *  sum eq of each charge number in solution...
- */
+	 */
 	charge_group_map.clear();
 	charge_group_map[0.0] = 0.0;
 
@@ -932,9 +955,10 @@ calc_init_donnan(void)
 			charge_group_map[s_x[i]->z] = s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl;
 		}
 	}
-/*
+	std::vector<LDBLE> zcorr(charge_group_map.size());
- *   calculate g for each surface...
+	/*
- */
+	 *   calculate g for each surface...
+	 */
 	for (int j = 0; j < count_unknowns; j++)
 	{
 		if (x[j]->type != SURFACE_CB)
@@ -948,14 +972,13 @@ calc_init_donnan(void)
 		{
 			f_psi = x[(size_t)j + 2]->master[0]->s->la * LOG_10;	/* -FPsi/RT */
 			f_psi = f_psi / 2;
-		} else
+		}
+		else
 			f_psi = x[j]->master[0]->s->la * LOG_10;
 		surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;
 
 		/* find psi_avg that matches surface charge... */
-/*    psi_avg = calc_psi_avg (0);
+		psi_avg = calc_psi_avg(charge_ptr, 0 * surf_chrg_eq, 0, zcorr);
-    appt 7/9/8... may have to change above one */
-		psi_avg = calc_psi_avg(charge_ptr, 0 * surf_chrg_eq);
 
 		/* fill in g's */
 		ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
@@ -986,7 +1009,7 @@ calc_init_donnan(void)
 			for (int i = 0; i < (int)this->s_x.size(); i++)
 			{
 				int is = s_x[i]->number;
-				assert (is < (int) s_diff_layer.size());
+				assert(is < (int)s_diff_layer.size());
 
 				s_diff_layer[is][charge_ptr->Get_name()].Set_g_moles(0.0);
 				s_diff_layer[is][charge_ptr->Get_name()].Set_dg_g_moles(0.0);
@@ -997,17 +1020,17 @@ calc_init_donnan(void)
 			std::string name = x[j]->master[0]->elt->name;
 			Utilities::replace("_psi", "", name);
 			output_msg(sformatf(
-					   "\nDonnan init on %s : charge, \tg, \tdg, Psi_surface = %8f V. \n",
+				"\nDonnan init on %s : charge, \tg, \tdg, \tzcorr, Psi_surface = %8f V. \n",
-					   name.c_str(),
+				name.c_str(),
-					   x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL *
+				x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL *
-					   tk_x / F_KJ_V_EQ));
+				tk_x / F_KJ_V_EQ));
 			for (std::map<LDBLE, cxxSurfDL>::iterator i_it = charge_ptr->Get_g_map().begin();
 				i_it != charge_ptr->Get_g_map().end(); i_it++)
 			{
-				output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n",
+				output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
-						   (double) i_it->first,
+					(double)i_it->first,
-						   (double) i_it->second.Get_g(),
+					(double)i_it->second.Get_g(),
-						   (double) i_it->second.Get_dg()));
+					(double)i_it->second.Get_dg()));
 			}
 		}
 	}
@@ -1015,13 +1038,13 @@ calc_init_donnan(void)
 }
 /* ---------------------------------------------------------------------- */
 LDBLE Phreeqc::
-calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq)
+calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std::vector<LDBLE> &zcorr)
 /* ---------------------------------------------------------------------- */
 {
-/*
+	/*
- * calculate the average (F * Psi / RT) that lets the DL charge counter the surface charge
+	 * calculate the average (F * Psi / RT) that lets the DL charge counter the surface charge
- */
+	 */
-	LDBLE fd, fd1, p, temp, ratio_aq;
+	LDBLE fd, fd1, p, /*psi_DL, */p_psi = R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ, temp, ratio_aq, z, z1, z1_c, eq, co_ion, sum_counter, sum_co;
 
 	ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
 	p = 0;
@@ -1031,30 +1054,79 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq)
 		p = -0.5 * log(-surf_chrg_eq * ratio_aq / mu_x + 1);
 	else if (surf_chrg_eq > 0)
 		p = 0.5 * log(surf_chrg_eq * ratio_aq / mu_x + 1);
-/*
+	/*
- * Optimize p in SS{s_x[i]->moles * z_i * g(p)} = -surf_chrg_eq
+	 * Optimize p in SS{s_x[i]->moles * z_i * g(p)} = -surf_chrg_eq
- *  g(p) = exp(-p * z_i) * ratio_aq
+	 *  g(p) = exp(-p * z_i) * ratio_aq
- * Elsewhere in PHREEQC, g is the excess, after subtraction of conc's for p = 0:
+	 * Elsewhere in PHREEQC, g is the excess, after subtraction of conc's for p = 0:
- *		      g(p) = (exp(-p *z_i) - 1) * ratio_aq
+	 *		      g(p) = (exp(-p *z_i) - 1) * ratio_aq
- */
+	 * with correct_GC true:
-	int l_iter = 0;
+	 * correct ions to better match the integrated concentrations:
+		 z == 1? z *= 0.285 cgc[6]
+		 z == 2? z *= 0.372 cgc[7]
+		 z == -1? z *= cgc[0] * (mu_x**(         cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * I ** cgc[4])
+		 z == -2? z *= cgc[0] * (mu_x**(cgc[5] * cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * I ** cgc[4])
+	 */
+
+	cxxSurface *surf_p = use.Get_surface_ptr();
+	bool correct_GC = surf_p->Get_correct_GC(),  local_correct_GC = correct_GC;
+	bool only_count = surf_p->Get_only_counter_ions();
+	LDBLE Gamma = fabs(surf_chrg_eq) / (charge_ptr->Get_specific_area() * charge_ptr->Get_grams()) / 1e-6,
+		cgc[10] = { 0.36, 0.1721, 0.798, 0.287, 0.1457, 1.2, 0.285, 0.372 };
+
+	if (!surf_p->Donnan_factors.empty())
+		std::copy(std::begin(surf_p->Donnan_factors), std::end(surf_p->Donnan_factors), cgc);
+
+	cgc[1] *= pow(nDbl, cgc[2]) * pow(Gamma, cgc[3]) * pow(mu_x, cgc[4]);
+
+	int l_iter = 0, z_iter;
+	sum_co = sum_counter = 0;
 	do
 	{
 		fd = surf_chrg_eq;
 		fd1 = 0.0;
+		z_iter = 0;
+		if (l_iter == 1 && local_correct_GC && fabs(sum_counter) < fabs(sum_co))
+		{
+			local_correct_GC = false;
+			l_iter = 0;
+		}
 		std::map<LDBLE, LDBLE>::iterator it;
 		for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
 		{
-			LDBLE z = it->first;
+			z = it->first;
-			if (!z || (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0))
+			z1 = z;
+			if (l_iter == 0) zcorr[z_iter] = z;
+			co_ion = surf_chrg_eq * z;
+			if (!z || (only_count && co_ion > 0))
+			{
+				z_iter++;
 				continue;
-			LDBLE eq = it->second;
+			}
-			/*  multiply with ratio_aq for multiplier options cp and cm
+			if (nDbl && local_correct_GC)
-				in calc_all_donnan (not used now)...  */
+			{
-			temp = exp(-z * p) * ratio_aq;
+				/*psi_DL = fabs(p * p_psi);*/
+				if (co_ion < 0)
+				{//counter-ion
+					if (fabs(z) > 1) temp = cgc[7];
+					else             temp = cgc[6];
+					sum_counter += z * temp;
+				}
+				else
+				{// co-ion
+					if (fabs(z) > 1) temp = cgc[0] * pow(mu_x, cgc[1] * cgc[5]);
+					else             temp = cgc[0] * pow(mu_x, cgc[1]);
+					sum_co += z * temp;
+				}
+				zcorr[z_iter] = z * temp;
+			}
+			z1 = zcorr[z_iter];
+			eq = it->second;
+			temp = exp(-z1 * p) * ratio_aq;
 
 			fd += eq * temp;
-			fd1 -= z * eq * temp;
+			fd1 -= z1 * eq * temp;
+			if (z == 1) z1_c = z1;
+			z_iter++;
 		}
 		fd /= -fd1;
 		p += (fd > 1) ? 1 : ((fd < -1) ? -1 : fd);
@@ -1076,12 +1148,11 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq)
 				(double)surf_chrg_eq, (double)charge_ptr->Get_mass_water());
 			error_msg(error_string, STOP);
 		}
-	}
+	} while (fabs(fd) > 1e-12 && p != 0.0);
-	while (fabs(fd) > 1e-12 && p != 0.0);
 	if (debug_diffuse_layer == TRUE)
 		output_msg(sformatf(
-				   "iter in calc_psi_avg = %d. g(+1) = %8f. surface charge = %12.4e.\n",
+			"iter in calc_psi_avg = %d. g(+1) = %8f, surface charge = %12.4e, psi_DL = %12.3e V.\n",
-				   l_iter, (double) (exp(-p) - 1), (double) surf_chrg_eq));
+			l_iter, (double)(exp(-p) - 1), (double)surf_chrg_eq, (double)(p * z1_c * p_psi)));
 
 	return (p);
 }

src/phreeqcpp/inverse.cpp

@@ -3586,8 +3586,10 @@ check_isotopes(class inverse *inv_ptr)
 			if (j < inv_ptr->i_u[i].uncertainties.size()
 				&& !isnan(inv_ptr->i_u[i].uncertainties[j]))
 #else
+			//if (j < inv_ptr->i_u[i].uncertainties.size()
+			//	&& inv_ptr->i_u[i].uncertainties[j] != NAN)
 			if (j < inv_ptr->i_u[i].uncertainties.size()
-				&& inv_ptr->i_u[i].uncertainties[j] != NAN)
+				&& !std::isnan(inv_ptr->i_u[i].uncertainties[j]))
 #endif
 			{
 				kit->second.Set_x_ratio_uncertainty(inv_ptr->i_u[i].uncertainties[j]);
@@ -3598,8 +3600,10 @@ check_isotopes(class inverse *inv_ptr)
 			else if (inv_ptr->i_u[i].uncertainties.size() > 0
 				&& !isnan(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1]))
 #else
+			//else if (inv_ptr->i_u[i].uncertainties.size() > 0
+			//	&& inv_ptr->i_u[i].uncertainties[(size_t)inv_ptr->i_u[i].uncertainties.size() - 1] != NAN)
 			else if (inv_ptr->i_u[i].uncertainties.size() > 0
-				&& inv_ptr->i_u[i].uncertainties[(size_t)inv_ptr->i_u[i].uncertainties.size() - 1] != NAN)
+				&& !std::isnan(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1]))
 #endif
 			{
 				kit->second.Set_x_ratio_uncertainty(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1]);
@@ -3609,7 +3613,8 @@ check_isotopes(class inverse *inv_ptr)
 #ifdef NPP
 			else if (!isnan(kit->second.Get_ratio_uncertainty()))
 #else
-			else if (kit->second.Get_ratio_uncertainty() != NAN)
+			//else if (kit->second.Get_ratio_uncertainty() != NAN)
+			else if (!std::isnan(kit->second.Get_ratio_uncertainty()))
 #endif
 			{
 				kit->second.Set_x_ratio_uncertainty(
@@ -3641,7 +3646,8 @@ check_isotopes(class inverse *inv_ptr)
 #ifdef NPP
 			if (isnan(kit->second.Get_x_ratio_uncertainty()))
 #else
-			if (kit->second.Get_x_ratio_uncertainty() == NAN)
+			//if (kit->second.Get_x_ratio_uncertainty() == NAN)
+			if (std::isnan(kit->second.Get_x_ratio_uncertainty()))
 #endif
 			{
 				error_string = sformatf(

src/phreeqcpp/isotopes.cpp

@@ -903,7 +903,8 @@ punch_calculate_values(void)
 #ifdef NPP
 			if (isnan(rate_moles))
 #else
-			if (rate_moles == NAN)
+			//if (rate_moles == NAN)
+			if (std::isnan(rate_moles))
 #endif
 			{
 				error_string = sformatf( "Calculated value not SAVEed for %s.",
@@ -1136,7 +1137,8 @@ calculate_values(void)
 #ifdef NPP
 				if (isnan(rate_moles))
 #else
-				if (rate_moles == NAN)
+				//if (rate_moles == NAN)
+				if (std::isnan(rate_moles))
 #endif
 				{
 					error_string = sformatf( "Calculated value not SAVEed for %s.",
@@ -1203,7 +1205,8 @@ calculate_values(void)
 #ifdef NPP
 				if (isnan(rate_moles))
 #else
-				if (rate_moles == NAN)
+				//if (rate_moles == NAN)
+				if (std::isnan(rate_moles))
 #endif
 				{
 					error_string = sformatf( "Calculated value not SAVEed for %s.",

src/phreeqcpp/kinetics.cpp

@@ -118,7 +118,8 @@ calc_kinetic_reaction(cxxKinetics *kinetics_ptr, LDBLE time_step)
 #ifdef NPP
 			if (isnan(rate_moles))
 #else
-			if (rate_moles == NAN)
+			//if (rate_moles == NAN)
+			if (std::isnan(rate_moles))
 #endif
 			{
 				error_string = sformatf( "Moles of reaction not SAVEed for %s.",

src/phreeqcpp/mainsubs.cpp

@@ -540,6 +540,8 @@ initial_exchangers(int print)
 			viscosity();
 			species_list_sort();
 			print_exchange();
+			if (pr.user_print)
+				print_user_print();
 			xexchange_save(n_user);
 			punch_all();
 			/* free_model_allocs(); */
@@ -671,7 +673,9 @@ initial_gas_phases(int print)
 			}
 
 			print_gas_phase();
- 			if (PR /*&& use.Get_gas_phase_ptr()->total_p > 1.0*/)
+			if (pr.user_print)
+				print_user_print();
+			if (PR /*&& use.Get_gas_phase_ptr()->total_p > 1.0*/)
  				warning_msg("While initializing gas phase composition by equilibrating:\n"
 				"         Found definitions of gas` critical temperature and pressure.\n"
 				"         Going to use Peng-Robinson in subsequent calculations.\n");
@@ -745,7 +749,8 @@ initial_surfaces(int print)
 			set_and_run_wrapper(-1, FALSE, FALSE, -1, 0.0);
 			species_list_sort();
 			print_surface();
-			/*print_all(); */
+			if (pr.user_print)
+				print_user_print();
 			punch_all();
 			xsurface_save(n_user);
 			/* free_model_allocs(); */
@@ -1257,6 +1262,7 @@ xsolution_save(int n_user)
 	temp_solution.Set_ah2o(ah2o_x);
 	//temp_solution.Set_density(density_x);
 	temp_solution.Set_density(calc_dens());
+	temp_solution.Set_viscosity(this->viscos);
 	temp_solution.Set_total_h(total_h_x);
 	temp_solution.Set_total_o(total_o_x);
 	temp_solution.Set_cb(cb_x);	/* cb_x does not include surface charge sfter sum_species */

src/phreeqcpp/model.cpp

@@ -2543,6 +2543,7 @@ int Phreeqc::
 calc_gas_pressures(void)
 /* ---------------------------------------------------------------------- */
 {
+	//int n_g = 0;
 	LDBLE lp, V_m = 0;
 	class rxn_token *rxn_ptr;
 	std::vector<class phase *> phase_ptrs;
@@ -2575,6 +2576,7 @@ calc_gas_pressures(void)
 			phase_ptrs.push_back(phase_ptr);
 			if (!PR && phase_ptr->t_c > 0 && phase_ptr->p_c > 0)
 				PR = true;
+			//n_g++;
 		}
 		if (iterations > 2 && gas_phase_ptr->Get_type() == cxxGasPhase::GP_VOLUME)
 		{
@@ -3776,6 +3778,7 @@ residuals(void)
 	int converge;
 
 	LDBLE l_toler;
+	//LDBLE sum_residual;
 	LDBLE sinh_constant;
 	LDBLE sum, sum1;
 	class master *master_ptr, *master_ptr1, *master_ptr2;
@@ -3783,6 +3786,7 @@ residuals(void)
 	int print_fail;
 	std::vector<LDBLE> cd_psi;
 	print_fail = FALSE;
+	//sum_residual = 0.0;
 	sigmaddl = 0;
 	sum = 0;
 /*
@@ -4526,6 +4530,7 @@ residuals(void)
  *   Store residuals in array
  */
 		my_array[((size_t)i + 1) * (count_unknowns + 1) - 1] = residual[i];
+		//sum_residual += fabs(residual[i]);
 	}
 /*
  *   Return

src/phreeqcpp/pitzer.cpp

@@ -52,6 +52,7 @@ pitzer_tidy(void)
 	const char *string1, *string2;
 	int i, j, order;
 	int i0, i1, i2;
+	//int count_pos, count_neg, count_neut, count[3], jj;
 	int count_neut, count[3], jj;
 	LDBLE z0, z1;
 	class pitz_param *pzp_ptr;
@@ -339,21 +340,22 @@ pitzer_tidy(void)
 			i0 = pitz_params[i]->ispec[0];
 			i1 = pitz_params[i]->ispec[1];
 			i2 = pitz_params[i]->ispec[2];
+			//count_pos = count_neg = count_neut = 0;
 			count_neut = 0;
 			for (j = 0; j <= 2; j++)
 			{
-				//if (spec[pitz_params[i]->ispec[j]]->z > 0)
+				if (spec[pitz_params[i]->ispec[j]]->z > 0)
-				//{
+				{
-				//	count_pos++;
+					//count_pos++;
-				//}
+				}
 				if (spec[pitz_params[i]->ispec[j]]->z == 0)
 				{
 					count_neut++;
 				}
-				//if (spec[pitz_params[i]->ispec[j]]->z < 0)
+				if (spec[pitz_params[i]->ispec[j]]->z < 0)
-				//{
+				{
-				//	count_neg++;
+					//count_neg++;
-				//}
+				}
 			}
 			/* All neutral */
 			if (count_neut == 3)

src/phreeqcpp/prep.cpp

@@ -1980,6 +1980,7 @@ get_list_master_ptrs(const char* cptr, class master *master_ptr)
  *   Output: space is allocated and a list of master species pointers is
  *           returned.
  */
+	//int j, l, count_list;
 	int j, l;
 	char token[MAX_LENGTH];
 	std::vector<class master*> master_ptr_list;
@@ -1987,6 +1988,7 @@ get_list_master_ptrs(const char* cptr, class master *master_ptr)
 /*
  *   Make list of master species pointers
  */
+	//count_list = 0;
 	//master_ptr_list = unknown_alloc_master();
 	master_ptr0 = master_ptr;
 	if (master_ptr0 == master_ptr->s->primary)
@@ -2146,6 +2148,7 @@ mb_for_species_aq(int n)
  *                             by coef, usually moles.
  *        mb_unknowns.coef - coefficient of s[n] in equation or relation
  */
+	//int i, j;
 	int i;
 	class master *master_ptr;
 	class unknown *unknown_ptr;
@@ -2222,6 +2225,7 @@ mb_for_species_aq(int n)
  */
 	if (use.Get_surface_ptr() != NULL && s[n]->type < H2O && dl_type_x != cxxSurface::NO_DL)
 	{
+		//j = 0;
 		for (i = 0; i < count_unknowns; i++)
 		{
 			if (x[i]->type == SURFACE_CB)
@@ -2233,6 +2237,7 @@ mb_for_species_aq(int n)
 
 				store_mb_unknowns(unknown_ptr, s_diff_layer[n][charge_ptr->Get_name()].Get_g_moles_address(),
 								  s[n]->z, s_diff_layer[n][charge_ptr->Get_name()].Get_dg_g_moles_address());
+				//j++;
 			}
 		}
 	}

src/phreeqcpp/print.cpp

@@ -339,8 +339,12 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr)
 		{
 			LDBLE exp_g = charge_ptr->Get_g_map()[1].Get_g() * mass_water_aq_x / mass_water_surface + 1;
 			LDBLE psi_DL = -log(exp_g) * R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ;
-			output_msg(sformatf(
+			if (use.Get_surface_ptr()->Get_correct_GC())
-				"\n\tTotal moles in diffuse layer (excluding water), Donnan calculation."));
+				output_msg(sformatf(
+					"\n\tTotal moles in diffuse layer (excluding water), Donnan corrected to match Poisson-Boltzmann."));
+			else
+				output_msg(sformatf(
+					"\n\tTotal moles in diffuse layer (excluding water), Donnan calculation."));
 			output_msg(sformatf(
 				"\n\tDonnan Layer potential, psi_DL = %10.3e V.\n\tBoltzmann factor, exp(-psi_DL * F / RT) = %9.3e (= c_DL / c_free if z is +1).\n\n",
 				psi_DL, exp_g));
@@ -2234,7 +2238,7 @@ print_totals(void)
 			   (double) calc_solution_volume()));
 	}
 /* VP: Density End */
-#ifdef NPP
+//#ifdef NPP
 	if (print_viscosity)
 	{
 		output_msg(sformatf("%45s%9.5f", "Viscosity (mPa s)  = ",
@@ -2251,7 +2255,7 @@ print_totals(void)
 		}
 		else output_msg(sformatf("\n"));
 	}
-#endif
+//#endif
 	output_msg(sformatf("%45s%7.3f\n", "Activity of water  = ",
 			   exp(s_h2o->la * LOG_10)));
 	output_msg(sformatf("%45s%11.3e\n", "Ionic strength (mol/kgw)  = ",

src/phreeqcpp/read.cpp

@@ -6283,7 +6283,7 @@ read_surface(void)
 	 *   ERROR   if error occurred reading data
 	 *
 	 */
-	int n_user;
+	int n_user, i1;
 	LDBLE conc;
 	const char* cptr, *cptr1;
 	std::string token, token1, name;
@@ -6306,9 +6306,10 @@ read_surface(void)
 		"ccm",                  /* 13 */
         "equilibrium",          /* 14 */
 		"site_units",           /* 15 */
-		"ddl"                   /* 16 */
+		"ddl",                  /* 16 */
+		"donnan_factors"        /* 17 */
 	};
-	int count_opt_list = 17;
+	int count_opt_list = 18;
 	/*
 	 * kin_surf is for Surfaces, related to kinetically reacting minerals
 	 *    they are defined if "sites" is followed by mineral name:
@@ -6413,7 +6414,7 @@ read_surface(void)
 							if (thickness != 0)
 							{
 								error_msg
-									("You must enter EITHER thickness OR Debye lengths (1/k),\n	   and relative DDL viscosity, DDL limit.\nCorrect is (for example): -donnan 1e-8 viscosity 0.5\n or (default values):     -donnan debye_lengths 1 viscosity 1 limit 0.8",
+									("You must enter EITHER thickness OR Debye lengths (1/k),\n	   and relative DDL viscosity, DDL limit.\nCorrect is (for example): -donnan 1e-8 viscosity 0.5 limit 0.9 correct_GC true\n or (default values):     -donnan debye_lengths 1 viscosity 1 limit 0.8 correct_GC false",
 									CONTINUE);
 								error_msg(line_save, CONTINUE);
 								input_error++;
@@ -6436,6 +6437,23 @@ read_surface(void)
 								break;
 							}
 						}
+						else if (token[0] == 'C' || token[0] == 'c')
+						{
+							copy_token(token1, &next_char);
+							if (token1[0] == 'T' || token1[0] == 't' || token1[0] == 'F' || token1[0] == 'f')
+							{
+								temp_surface.Set_correct_GC(get_true_false(token1.c_str(), TRUE) == TRUE);
+								continue;
+							} else
+							{
+								error_msg
+									("Expected True or False for correct_GC (which brings co-ion concentrations closer to their integrated double layer value).",
+									CONTINUE);
+								error_msg(line_save, CONTINUE);
+								input_error++;
+								break;
+							}
+						}
 						else if (token[0] == 'V' || token[0] == 'v')
 						{
 							int j = copy_token(token1, &next_char);
@@ -6570,6 +6588,31 @@ read_surface(void)
 		case 16:				/* ddl */
 			temp_surface.Set_type(cxxSurface::DDL);
 			break;
+		case 17:			/* Donnan_factors */
+			temp_surface.Donnan_factors.clear();
+			i1 = 0;
+			for (;;)
+			{
+				int i = copy_token(token, &next_char);
+				if (i == DIGIT && i1 < 8)
+				{
+					(void)sscanf(token.c_str(), SCANFORMAT, &dummy);
+					temp_surface.Donnan_factors.push_back(dummy);
+					i1++;
+					continue;
+				}
+				else if (i != EMPTY || i1 > 8)
+				{
+					error_msg
+					("Expected at most 8 numbers for the Donnan_factors for co- and counter-ions,\n 	z *= cgc[0] * (mu_x**(cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * mu_x**(cgc[4])",
+						CONTINUE);
+					error_msg(line_save, CONTINUE);
+					input_error++;
+					break;
+				}
+				break;
+			}
+			break;
 		case OPTION_DEFAULT:
 			/*
 			 *   Read surface component

src/phreeqcpp/spread.cpp

@@ -1141,6 +1141,7 @@ copy_token_tab(std::string& token, const char **cptr)
  *      EOL,
  *      UNKNOWN.
  */
+	//int i, return_value;
 	int return_value;
 	char c;
 /*
@@ -1180,6 +1181,7 @@ copy_token_tab(std::string& token, const char **cptr)
 /*
  *   Begin copying to token
  */
+	//i = 0;
 	for (;;)
 	{
 		c = **cptr;
@@ -1196,6 +1198,7 @@ copy_token_tab(std::string& token, const char **cptr)
 		{
 			token.push_back(c);
 			(*cptr)++;
+			//i++;
 		}
 	}
 	return (return_value);

src/phreeqcpp/sundialsmath.cpp

@@ -59,7 +59,8 @@
 
 
 #include <stdio.h>
-#include <math.h>
+ //#include <math.h>
+#include <cmath>
 #include "sundialsmath.h"
 #include "sundialstypes.h"
 

src/phreeqcpp/tally.cpp

@@ -798,6 +798,7 @@ build_tally_table(void)
  */
 	int j, k, l, p, save_print_use;
 	size_t n;
+	//int count_tt_pure_phase, count_tt_ss_phase, count_tt_kinetics;
 	class phase *phase_ptr;
 	char token[MAX_LENGTH];
 	const char* cptr;
@@ -870,6 +871,7 @@ build_tally_table(void)
 /*
  *   Count pure phases
  */
+	//count_tt_pure_phase = 0;
 	if (Rxn_pp_assemblage_map.size() > 0)
 	{
 		/* 
@@ -902,6 +904,7 @@ build_tally_table(void)
 				/*
 				 * Add to table
 				 */
+				//count_tt_pure_phase++;
 				n = count_tally_table_columns;
 				extend_tally_table();
 				tally_table[n].name = phase_ptr->name;
@@ -928,6 +931,7 @@ build_tally_table(void)
 /*
  *   Add solid-solution pure phases
  */
+	//count_tt_ss_phase = 0;
 	if (Rxn_ss_assemblage_map.size() > 0)
 	{
 		/* 
@@ -960,6 +964,7 @@ build_tally_table(void)
 					/*
 					 * Add to table
 					 */
+					//count_tt_ss_phase++;
 					n = count_tally_table_columns;
 					extend_tally_table();
 					tally_table[n].name = phase_ptr->name;
@@ -977,6 +982,7 @@ build_tally_table(void)
 /*
  *   Add kinetic reactants
  */
+	//count_tt_kinetics = 0;
 	if (Rxn_kinetics_map.size() > 0)
 	{
 		std::map<int, cxxKinetics>::iterator it;
@@ -1000,6 +1006,7 @@ build_tally_table(void)
 				/*
 				 * Add to table
 				 */
+				//count_tt_kinetics++;
 				n = count_tally_table_columns;
 				extend_tally_table();
 				tally_table[n].name = string_hsave(kinetics_comp_ptr->Get_rate_name().c_str());

src/phreeqcpp/tidy.cpp

@@ -982,7 +982,8 @@ tidy_gas_phase(void)
 #ifdef NPP
 					if (!isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()))
 #else
-					if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN)
+					//if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN)
+					if (!std::isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()))
 #endif
 					{
 						P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read();
@@ -1015,7 +1016,8 @@ tidy_gas_phase(void)
 #ifdef NPP
 						if (!isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()))
 #else
-						if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN)
+						//if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN)
+						if (!std::isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()))
 #endif
 						{
 							P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read();
@@ -1690,7 +1692,8 @@ tidy_ss_assemblage(void)
 #ifdef NPP
 				if (isnan(comp_ptr->Get_moles()))
 #else
-				if (comp_ptr->Get_moles() == NAN)
+				//if (comp_ptr->Get_moles() == NAN)
+				if (std::isnan(comp_ptr->Get_moles()))
 #endif
 				{
 					input_error++;
@@ -3021,7 +3024,8 @@ tidy_isotopes(void)
 #ifdef NPP
 				if (!isnan(master[k]->isotope_ratio_uncertainty))
 #else
-				if (master[k]->isotope_ratio_uncertainty != NAN)
+				//if (master[k]->isotope_ratio_uncertainty != NAN)
+				if (!std::isnan(master[k]->isotope_ratio_uncertainty))
 #endif
 				{
 					temp_iso.Set_ratio_uncertainty_defined(true);

src/phreeqcpp/transport.cpp

@@ -1888,16 +1888,17 @@ fill_spec(int l_cell_no, int ref_cell)
 	if (por_il < interlayer_Dpor_lim)
 		por_il = viscos_il_f = 0.0;
 	/*
-	* correct diffusion coefficient for temperature and viscosity, D_T = D_298 * Tk * viscos_298 / (298 * viscos)
+	* correct diffusion coefficient for temperature and viscosity, D_T = D_298 * viscos_298 / viscos
 	*   modify viscosity effect: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15), SC data from Robinson and Stokes, 1959
 	*/
 	viscos = viscos_0;
 	/*
 	* put temperature factor in por_factor which corrects for porous medium...
 	*/
-	viscos_f *= tk_x * viscos_0_25 / (298.15 * viscos);
+	dum = viscos_0_25 / viscos;
-	viscos_il_f *= tk_x * viscos_0_25 / (298.15 * viscos);
+	viscos_f *= dum;
-	sol_D[l_cell_no].viscos_f = tk_x * viscos_0_25 / (298.15 * viscos);
+	viscos_il_f *= dum;
+	sol_D[l_cell_no].viscos_f = dum;
 
 	count_spec = count_exch_spec = 0;
 	/*
@@ -1908,6 +1909,8 @@ fill_spec(int l_cell_no, int ref_cell)
 		qsort(&species_list[0], species_list.size(),
 			sizeof(class species_list), sort_species_name);
 	}
+	if (correct_Dw)
+		calc_SC();
 
 	for (i = 0; i < (int)species_list.size(); i++)
 	{
@@ -2097,7 +2100,7 @@ fill_spec(int l_cell_no, int ref_cell)
 			}
 			if (correct_Dw)
 			{
-				calc_SC(); // note that neutral species are corrected as if z = 1, but is viscosity-dependent
+				//calc_SC(); // removed that neutral species are corrected as if z = 1, but is viscosity-dependent
 				sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw_corr * viscos_f;
 			}
 			if (l_cell_no <= count_cells + 1 && sol_D[l_cell_no].spec[count_spec].Dwt * pow(por, multi_Dn) > diffc_max)
@@ -2962,20 +2965,20 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 			if (!strcmp(ct[icell].m_s[cp].name, "H"))
 			{
 				dummy = ct[icell].m_s[cp].tot1;
-				if (dV_dcell || (icell > 0 && icell <= ilast))
+				if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
 					sptr1->Set_total_h(sptr1->Get_total_h() - dummy);
-				if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
+				if (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
 					sptr2->Set_total_h(sptr2->Get_total_h() + dummy);
 				if (sptr_stag)
 				{
 					dummy = ct[icell].m_s[cp].tot_stag;
-					if (dV_dcell || (icell > 0 && icell <= ilast))
+					if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
 						sptr1->Set_total_h(sptr1->Get_total_h() - dummy);
 					if (icell == c)
 					{
 						// mix the boundary solution with the stagnant column end-cell
 						dummy += ct[icell + 1].m_s[cp].tot_stag;
-						if (dV_dcell)
+						if (dV_dcell || fix_current)
 							sptr2->Set_total_h(sptr2->Get_total_h() - ct[icell + 1].m_s[cp].tot_stag);
 					}
 					sptr_stag->Set_total_h(sptr_stag->Get_total_h() + dummy);
@@ -2985,19 +2988,19 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 			if (!strcmp(ct[icell].m_s[cp].name, "O"))
 			{
 				dummy = ct[icell].m_s[cp].tot1;
-				if (dV_dcell || (icell > 0 && icell <= ilast))
+				if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
 					sptr1->Set_total_o(sptr1->Get_total_o() - dummy);
-				if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
+				if (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
 					sptr2->Set_total_o(sptr2->Get_total_o() + dummy);
 				if (sptr_stag)
 				{
 					dummy = ct[icell].m_s[cp].tot_stag;
-					if (dV_dcell || (icell > 0 && icell <= ilast))
+					if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
 						sptr1->Set_total_o(sptr1->Get_total_o() - dummy);
 					if (icell == c)
 					{
 						dummy += ct[icell + 1].m_s[cp].tot_stag;
-						if (dV_dcell)
+						if (dV_dcell || fix_current)
 							sptr2->Set_total_o(sptr2->Get_total_o() - ct[icell + 1].m_s[cp].tot_stag);
 					}
 					sptr_stag->Set_total_o(sptr_stag->Get_total_o() + dummy);
@@ -3018,7 +3021,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 			
 			// check for negative moles, add moles from other redox states and the donnan layer when necessary and available...
 			if (dum1 - ct[icell].m_s[cp].tot1 - ct[icell].m_s[cp].tot_stag < min_mol &&
-			(dV_dcell || (icell > 0 && icell <= ilast)))
+			(dV_dcell || fix_current || (icell > 0 && icell <= ilast)))
 			{
 				dum1 = moles_from_redox_states(sptr1, ct[icell].m_s[cp].name);
 				if (ct[icell].dl_s > 1e-8)
@@ -3039,7 +3042,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 			if (icell == c && sptr_stag && ct[c1].m_s[cp].tot_stag)
 				dum = ct[c1].m_s[cp].tot_stag;
 			if (dum2 + ct[icell].m_s[cp].tot2 - dum < min_mol &&
-			(dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)))
+			(dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)))
 			{
 				dum2 = moles_from_redox_states(sptr2, ct[icell].m_s[cp].name);
 				if (ct[icell + 1].dl_s > 1e-8)
@@ -3057,7 +3060,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 			if (fabs(ct[icell].m_s[cp].tot2) < fabs(ct[icell].m_s[cp].tot1))
 				ct[icell].m_s[cp].tot1 = ct[icell].m_s[cp].tot2;
 			
-			if (dV_dcell || (icell > 0 && icell <= ilast))
+			if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
 			{
 				dum = ct[icell].m_s[cp].tot1;
 				if (stagnant)
@@ -3078,10 +3081,10 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 				//ct[icell].J_ij_sum -= dum * ct[icell].m_s[cp].charge;
 			}
 			
-			if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
+			if (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
 			{
 				dum = ct[icell].m_s[cp].tot1;
-				if (stagnant && icell == c && dV_dcell)
+				if (stagnant && icell == c && (dV_dcell || fix_current))
 					dum -= ct[c1].m_s[cp].tot_stag;
 				dum2 += dum;
 				sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : min_mol);
@@ -3134,7 +3137,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 			}
 			
 			// reduce oscillations in the column-boundary cells, but not for H and O, and current_A is not adjusted...
-			if (dV_dcell && icell == il1 - incr && dV_dcell * ct[0].m_s[cp].charge < 0 && strcmp(ct[0].m_s[cp].name, "H") && strcmp(ct[0].m_s[cp].name, "O") && c > 3 && mixrun > 1)
+			if ((dV_dcell/* || fix_current*/) && icell == il1 - incr && dV_dcell * ct[0].m_s[cp].charge < 0 && strcmp(ct[0].m_s[cp].name, "H") && strcmp(ct[0].m_s[cp].name, "O") && c > 3 && mixrun > 1)
 			{
 				dummy = Utilities::Rxn_find(Rxn_solution_map, 0)->Get_totals()[ct[0].m_s[cp].name] / ct[0].kgw * (1 - ct[0].dl_s);
 				if (dummy > 1e-6)
@@ -4828,10 +4831,10 @@ Step (from cell 1 to count_cells + 1):
 	cxxSurface *surface_ptr1, *surface_ptr2;
 	LDBLE viscos_f;
 	/*
-	* temperature and viscosity correction for MCD coefficient, D_T = D_298 * Tk * viscos_298 / (298 * viscos)
+	* temperature and viscosity correction for MCD coefficient, D_T = D_298 * viscos_298 / viscos
 	*/
 	viscos_f = viscos_0;
-	viscos_f = tk_x * viscos_0_25 / (298.15 * viscos_f);
+	viscos_f = viscos_0_25 / viscos_f;
 
 	//n1 = 0;
 	//n2 = n1 + 1;
@@ -5502,7 +5505,7 @@ diff_stag_surf(int mobile_cell)
 	* temperature and viscosity correction for MCD coefficient, D_T = D_298 * Tk * viscos_298 / (298 * viscos)
 	*/
 	viscos_f = viscos_0;
-	viscos_f = tk_x * viscos_0_25 / (298.15 * viscos_f);
+	viscos_f = viscos_0_25 / viscos_f;
 
 	cxxSurface surface_n1, surface_n2;
 	cxxSurface *surface_n1_ptr = &surface_n1;
@@ -5991,8 +5994,7 @@ viscosity(void)
 			}
 			// find B * m and D * m * mu^d3
 			Bc += (s_x[i]->Jones_Dole[0] +
-				s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) *
+				s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * t1;
-				t1;
 			// define f_I from the exponent of the D * m^d3 term...
 			if (s_x[i]->Jones_Dole[5] >= 1)
 				t2 = mu_x / 3 / s_x[i]->Jones_Dole[5];
@@ -6074,8 +6076,8 @@ viscosity(void)
 		viscos += (viscos_0 * (Bc + Dc));
 	if (viscos < 0)
 	{
-		viscos = 0;
+		viscos = viscos_0;
-		warning_msg("viscosity < 0, reset to 0.");
+		warning_msg("viscosity < 0, reset to viscosity of pure water");
 	}
 	return viscos;