Merge commit '6e50709a45ab4533b39686ae23f147b1d08d455d'

src/phast/PhreeqcRM/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Solution.cxx

@@ -48,6 +48,7 @@ cxxSolution::cxxSolution(PHRQ_io * io)
 	this->cb = 0.0;
 	this->density = 1.0;
 	this->viscosity = 1.0;
+	this->viscos_0 = 1.0;
 	this->mass_water = 1.0;
 	this->soln_vol = 1.0;
 	this->total_alkalinity = 0.0;
@@ -82,6 +83,7 @@ cxxSolution::operator =(const cxxSolution &rhs)
 		this->total_o                    = rhs.total_o;
 		this->density                    = rhs.density;
 		this->viscosity                  = rhs.viscosity;
+		this->viscos_0                   = rhs.viscos_0;
 		this->cb                         = rhs.cb;
 		this->mass_water                 = rhs.mass_water;
 		this->soln_vol                   = rhs.soln_vol;
@@ -274,6 +276,8 @@ cxxSolution::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) con
 	// new identifier
 	s_oss << indent1;
 	s_oss << "-viscosity                 " << this->viscosity << "\n";
+	s_oss << indent1;
+	s_oss << "-viscos_0                  " << this->viscos_0 << "\n";
 
 	// soln_total conc structures
 	s_oss << indent1;
@@ -1086,6 +1090,16 @@ cxxSolution::read_raw(CParser & parser, bool check)
 			}
 			opt_save = CParser::OPT_DEFAULT;
 			break;
+		case 29:				// viscos_0
+			if (!(parser.get_iss() >> this->viscos_0))
+			{
+				this->viscos_0 = 1.0;
+				parser.incr_input_error();
+				parser.error_msg("Expected numeric value for viscos_0.",
+					PHRQ_io::OT_CONTINUE);
+			}
+			opt_save = CParser::OPT_DEFAULT;
+			break;
 		}
 		if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD)
 			break;
@@ -1415,6 +1429,7 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive)
 	this->cb += addee.cb * extensive;
 	this->density = f1 * this->density + f2 * addee.density;
 	this->viscosity = f1 * this->viscosity + f2 * addee.viscosity;
+	this->viscos_0 = f1 * this->viscos_0 + f2 * addee.viscos_0;
 	this->patm = f1 * this->patm + f2 * addee.patm;
 	// this->potV = f1 * this->potV + f2 * addee.potV; // appt
 	this->mass_water += addee.mass_water * extensive;
@@ -1773,6 +1788,7 @@ const std::vector< std::string >::value_type temp_vopts[] = {
 	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("viscosity")                             // 28
+	std::vector< std::string >::value_type("viscosity"),                            // 28
+	std::vector< std::string >::value_type("viscos_0")                              // 29
 };									   
 const std::vector< std::string > cxxSolution::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]);	

src/phast/PhreeqcRM/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/Solution.h

@@ -51,6 +51,8 @@ class cxxSolution:public cxxNumKeyword
 	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_viscos_0() const                          { return this->viscos_0; }
+	void Set_viscos_0(LDBLE l_viscos_0)                 { this->viscos_0 = l_viscos_0; }
 	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;}
@@ -134,6 +136,7 @@ class cxxSolution:public cxxNumKeyword
 	LDBLE mass_water;
 	LDBLE density;
 	LDBLE viscosity;
+	LDBLE viscos_0;
 	LDBLE soln_vol;
 	LDBLE total_alkalinity;
 	cxxNameDouble totals;

src/phast/PhreeqcRM/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/global_structures.h

@@ -109,7 +109,7 @@
 #define TRANSPORT		 8
 #define PHAST		     9
 
-/* constants in mass balance */
+/* constraints in mass balance */
 #define EITHER 0
 #define DISSOLVE 1
 #define PRECIPITATE -1
@@ -1486,7 +1486,9 @@ class spec
 		c = 0;
 		// charge number
 		z = 0;
-		// temperature corrected free water diffusion coefficient, m2/s
+		// free water diffusion coefficient, m2/s
+		Dw = 0;
+		// temperature and viscosity corrected free water diffusion coefficient, m2/s
 		Dwt = 0;
 		// temperature factor for Dw
 		dw_t = 0;
@@ -1503,6 +1505,7 @@ class spec
 	LDBLE lg;
 	LDBLE c;
 	LDBLE z;
+	LDBLE Dw;
 	LDBLE Dwt;
 	LDBLE dw_t;
 	LDBLE dw_a_v_dif;
@@ -1521,17 +1524,17 @@ class sol_D
 		count_exch_spec = 0;
 		// total moles of X-, max X- in transport step in sol_D[1], tk
 		exch_total = 0, x_max = 0, tk_x = 0;
-		// (tk_x * viscos_0_25) / (298 * viscos_0) 
-		viscos_f0 = 0;
-		// (viscos_0) / (298 * viscos) 
-		viscos_f = 0;
+		// viscos_0 at I = 0
+		viscos_0 = 0;
+		// viscosity of solution
+		viscos = 0;
 		spec = NULL;
 		spec_size = 0;
 	}
 	int count_spec;
 	int count_exch_spec;
 	LDBLE exch_total, x_max, tk_x;
-	LDBLE viscos_f0, viscos_f;
+	LDBLE viscos_0, viscos;
 	class spec* spec;
 	int spec_size;
 };

src/phast/PhreeqcRM/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/kinetics.cpp

@@ -1581,7 +1581,8 @@ set_and_run(int i, int use_mix, int use_kinetics, int nsaver,
 	sum_species();
 	viscos = viscosity(NULL);
 	use.Get_solution_ptr()->Set_viscosity(viscos);
-	if (use.Get_surface_ptr() != NULL && dl_type_x != cxxSurface::NO_DL)
+	use.Get_solution_ptr()->Set_viscos_0(viscos_0);
+	if (use.Get_surface_ptr() != NULL && dl_type_x != cxxSurface::NO_DL && use.Get_surface_ptr()->Get_calc_viscosity())
 		use.Get_surface_ptr()->Set_DDL_viscosity(viscosity(use.Get_surface_ptr()));
 	return (converge);
 }

src/phast/PhreeqcRM/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/mainsubs.cpp

@@ -438,6 +438,7 @@ initial_solutions(int print)
 			sum_species();
 			viscos = viscosity(NULL);
 			use.Get_solution_ptr()->Set_viscosity(viscos);
+			use.Get_solution_ptr()->Set_viscos_0(viscos_0);
 			if (use.Get_surface_ptr() != NULL && dl_type_x != cxxSurface::NO_DL)
 				use.Get_surface_ptr()->Set_DDL_viscosity(viscosity(use.Get_surface_ptr()));
 			add_isotopes(solution_ref);
@@ -1267,6 +1268,7 @@ xsolution_save(int n_user)
 	// the subroutine is called at the start of a new simulation, and the following 2 go wrong since s_x is not updated 
 	temp_solution.Set_density(density_x);
 	temp_solution.Set_viscosity(viscos);
+	temp_solution.Set_viscos_0(viscos_0);
 	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 after sum_species */

src/phast/PhreeqcRM/src/IPhreeqcPhast/IPhreeqc/phreeqcpp/step.cpp

@@ -328,6 +328,7 @@ xsolution_zero(void)
 	mu_x = 0.0;
 	ah2o_x = 0.0;
 	viscos = 0.0;
+	viscos_0 = 0.0;
 	density_x = 0.0;
 	total_h_x = 0.0;
 	total_o_x = 0.0;
@@ -381,6 +382,7 @@ add_solution(cxxSolution *solution_ptr, LDBLE extensive, LDBLE intensive)
 	mu_x += solution_ptr->Get_mu() * intensive;
 	ah2o_x += solution_ptr->Get_ah2o() * intensive;
 	viscos += solution_ptr->Get_viscosity() * intensive;
+	viscos_0 += solution_ptr->Get_viscos_0() * intensive;
 	density_x += solution_ptr->Get_density() * intensive;
 
 	total_h_x += solution_ptr->Get_total_h() * extensive;