Merge commit '27025d1ccd77fdf67a75e7dc1e49394daa4d6a0d'

Author: github-actions[bot]

IPhreeqcMMS/IPhreeqc/src/phreeqcpp/transport.cpp

@@ -6070,21 +6070,22 @@ viscosity(cxxSurface *surf_ptr)
 	Jones_Dole[6] contains the anion factor, 1 for Cl-, variable for other anions
 	f_z = (z * z + |z|) / 2, the contribution of the ion to mu_x, if z = 0: f_z = mu_x / m_i
 	f_I = variable, depends on d3_i > 1, or d3_i < 1.
+	tc is limited to 200 C.
+
 
 	A from Falkenhagen-Dole for a salt:
 	A = 4.3787e-14 * TK**1.5 / (eps_r**0.5)* (z1 + z2)**-0.5 / (D1 * D2) * psi
 	psi = (D1*z2 + D2*z1)/4 - z1*z2 * (D1-D2)**2 / ((D1*z1 + D2*z2)**0.5 + ((D1 + D2) * (z1 + z2))**0.5)**2
 	D1, z1 for the cation, D2, |z2| for the anion of the salt.
 	We use the harmonic mean of the Dw's, and the arithmetic mean of the z's,
 	both weighted by the equivalent concentration.
-
-	tc is limited to 200 C.
 	*/
 	LDBLE D1, D2, z1, z2, m_plus, m_min, eq_plus, eq_min, eq_dw_plus, eq_dw_min, t1, t2, t3, fan = 1;
-	LDBLE A, psi, Bc = 0, Dc = 0, Dw = 0.0, l_z, f_z, lm, V_an, m_an, V_Cl, l_moles, l_water, l_mu_x, dw_t_visc;
+	LDBLE A, psi, Bc = 0, Dc = 0, Dw = 0.0, l_z, f_z, lm, V_an, m_an, V_Cl, tc, l_moles, l_water, l_mu_x, dw_t_visc;
 
 	m_plus = m_min = eq_plus = eq_min = eq_dw_plus = eq_dw_min = V_an = m_an = V_Cl = 0;
 
+	tc = (tc_x > 200) ? 200 : tc_x;
 	l_water = mass_water_aq_x;
 	l_mu_x = mu_x;
 
@@ -6103,8 +6104,6 @@ viscosity(cxxSurface *surf_ptr)
 	LDBLE ratio_surf_aq = mass_water_surfaces_x / mass_water_aq_x;
 	for (i1 = 0; i1 < i1_last; i1++)
 	{
-		if (tc_x > 200)
-			goto highT;
 		Bc = Dc = Dw = m_plus = m_min = eq_plus = eq_min = eq_dw_plus = eq_dw_min = V_an = m_an = 0;
 		if (surf_ptr)
 		{
@@ -6182,7 +6181,7 @@ viscosity(cxxSurface *surf_ptr)
 				}
 				// find B * m and D * m * mu^d3
 				dw_t_visc = (s_x[i]->Jones_Dole[0] +
-					s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc_x)) * t1;
+					s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * t1;
 				Bc += dw_t_visc;
 				// define f_I from the exponent of the D * m^d3 term...
 				if (s_x[i]->Jones_Dole[5] >= 1)
@@ -6191,7 +6190,7 @@ viscosity(cxxSurface *surf_ptr)
 					t2 = -0.8 / s_x[i]->Jones_Dole[5];
 				else
 					t2 = -1;
-				t3 = (s_x[i]->Jones_Dole[3] * exp(-s_x[i]->Jones_Dole[4] * tc_x)) *
+				t3 = (s_x[i]->Jones_Dole[3] * exp(-s_x[i]->Jones_Dole[4] * tc)) *
 					t1 * (pow(l_mu_x, s_x[i]->Jones_Dole[5])*(1 + t2) + pow(t1 * f_z, s_x[i]->Jones_Dole[5])) / (2 + t2);
 				if (t3 < -1e-5)
 					t3 = 0;
@@ -6260,35 +6259,28 @@ viscosity(cxxSurface *surf_ptr)
 			V_an /= m_an;
 		if (!V_Cl)
 			V_Cl = calc_vm_Cl();
-		if (V_an == V_Cl)
-			fan = 1;
-		else if (V_Cl >= 1)
+		if (V_an && V_Cl)
 			fan = 2 - V_an / V_Cl;
 		else
-			fan = 2 - V_an * mu_x;
-		//if (Dc < 0) 
-		//	Dc = 0; // provisional...
+			fan = 1;
+		if (Dc < 0) 
+			Dc = 0; // provisional...
 		viscos += viscos_0 * fan * (Bc + Dc);
 		if (viscos < 0)
 		{
 			viscos = viscos_0;
 			warning_msg("viscosity < 0, reset to viscosity of pure water");
 		}
-
-	highT:
 		if (!surf_ptr)
 		{
-			t1 = fabs(Bc + Dc);
 			for (i = 0; i < (int)this->s_x.size(); i++)
 			{
 				if (s_x[i]->type != AQ && s_x[i]->type > HPLUS)
 					continue;
-				if (s_x[i]->lm < -9)
+				if ((lm = s_x[i]->lm) < -9)
 					continue;
-				if (tc_x > 200)
-					s_x[i]->dw_t_visc = 0;
-				else if (s_x[i]->Jones_Dole[0] || s_x[i]->Jones_Dole[1] || s_x[i]->Jones_Dole[3])
-					s_x[i]->dw_t_visc /= t1;
+				if (s_x[i]->Jones_Dole[0] || s_x[i]->Jones_Dole[1] || s_x[i]->Jones_Dole[3])
+					s_x[i]->dw_t_visc /= (Bc + Dc);
 			}
 		}
 		else //if (surf_ptr->Get_calc_viscosity())

IPhreeqcMMS/IPhreeqc/src/phreeqcpp/utilities.cpp

@@ -192,12 +192,11 @@ calc_rho_0(LDBLE tc, LDBLE pa)
 	LDBLE p1 = -6.0251348E-06 + tc * ( 3.6696407E-07 + tc * (-9.2056269E-09 + tc * ( 6.7024182E-11 + tc * -1.5947241E-13)));
 	LDBLE p2 = -2.2983596E-09 + tc * (-4.0133819E-10 + tc * ( 1.2619821E-11 + tc * (-9.8952363E-14 + tc *  2.3363281E-16)));
 	LDBLE p3 =  7.0517647E-11 + tc * ( 6.8566831E-12 + tc * (-2.2829750E-13 + tc * ( 1.8113313E-15 + tc * -4.2475324E-18)));
-	// /* The minimal pressure equals the saturation pressure... */
-	   /* 5/12/2025: remove ah2o_x, rho_0 is independent of concentrations*/
-	//if (ah2o_x <= 1.0)
-	//	p_sat = exp(11.6702 - 3816.44 / (T - 46.13)) * ah2o_x;
-	//else
-	p_sat = exp(11.6702 - 3816.44 / (T - 46.13));
+	/* The minimal pressure equals the saturation pressure... */
+	if (ah2o_x <= 1.0)
+		p_sat = exp(11.6702 - 3816.44 / (T - 46.13)) * ah2o_x;
+	else
+		p_sat = exp(11.6702 - 3816.44 / (T - 46.13));
 	//ah2o_x0 = ah2o_x; // for updating rho in model(): compare with new ah2o_x
 	if (pa < p_sat || (use.Get_solution_ptr() && use.Get_solution_ptr()->Get_patm() < p_sat))
 	{