diff --git a/phreeqcpp/ChartHandler.cpp b/phreeqcpp/ChartHandler.cpp index b9236865..df2ee816 100644 --- a/phreeqcpp/ChartHandler.cpp +++ b/phreeqcpp/ChartHandler.cpp @@ -170,18 +170,14 @@ ChartHandler::End_timer() io->error_flush(); } } - size_t i(0), i2(0); for ( ; it != chart_map.end(); it++) { - i = 0; it->second->Rate_free(); if (it->second->Get_form_started()) { #if defined(__cplusplus_cli) while (0 != System::Threading::Interlocked::CompareExchange(it->second->usingResource, 6, 0)) { - //if (i > max_tries) break; - i++; System::Threading::Thread::Sleep(60); } #endif @@ -191,12 +187,8 @@ ChartHandler::End_timer() int n = System::Threading::Interlocked::Exchange(it->second->usingResource, 0); assert(n == 6); #endif - - i2 = 0; while (it->second->Get_done() != true) { - //if (i2 > max_tries) break; - i2++; #if defined(__cplusplus_cli) System::Threading::Thread::Sleep(60); #endif @@ -225,7 +217,6 @@ ChartHandler::dump(std::ostream & oss, unsigned int indent) std::map::iterator it = this->chart_map.begin(); for ( ; it != chart_map.end(); it++) { - size_t i = 0; it->second->dump(oss, indent); } return true; diff --git a/phreeqcpp/ChartObject.cpp b/phreeqcpp/ChartObject.cpp index 018e3b25..ae5004c9 100644 --- a/phreeqcpp/ChartObject.cpp +++ b/phreeqcpp/ChartObject.cpp @@ -11,7 +11,6 @@ #include "ChartObject.h" #include "Parser.h" #include -//#include #include #include #include "phqalloc.h" diff --git a/phreeqcpp/PBasic.cpp b/phreeqcpp/PBasic.cpp index b76897ad..a103b631 100644 --- a/phreeqcpp/PBasic.cpp +++ b/phreeqcpp/PBasic.cpp @@ -1558,6 +1558,9 @@ listtokens(FILE * f, tokenrec * l_buf) case tokt_sc: output_msg("T_SC"); break; + case tokf_visc: + output_msg("F_VISC"); + break; case toktc: output_msg("TC"); break; @@ -3904,6 +3907,13 @@ factor(struct LOC_exec * LINK) } break; + case tokf_visc: + { + const char* str = stringfactor(STR1, LINK); + n.UU.val = (parse_all) ? 1 : PhreeqcPtr->calc_f_visc(str); + } + break; + case toktc: { n.UU.val = PhreeqcPtr->tc_x; @@ -7518,6 +7528,7 @@ const std::map::value_type temp_tokens[] std::map::value_type("surf", PBasic::toksurf), std::map::value_type("sys", PBasic::toksys), std::map::value_type("t_sc", PBasic::tokt_sc), + std::map::value_type("f_visc", PBasic::tokf_visc), std::map::value_type("tc", PBasic::toktc), std::map::value_type("time", PBasic::toktime), std::map::value_type("title", PBasic::toktitle), diff --git a/phreeqcpp/PBasic.h b/phreeqcpp/PBasic.h index 73321cc8..3035ceb8 100644 --- a/phreeqcpp/PBasic.h +++ b/phreeqcpp/PBasic.h @@ -7,7 +7,6 @@ #include #include #include -//#include #include #include #include "phrqtype.h" @@ -338,6 +337,7 @@ public: toktotmol, toktotmoles, toktrim, + tokf_visc, tokviscos, tokviscos_0, tokvm, diff --git a/phreeqcpp/Phreeqc.h b/phreeqcpp/Phreeqc.h index 4173da56..f2ad5d52 100644 --- a/phreeqcpp/Phreeqc.h +++ b/phreeqcpp/Phreeqc.h @@ -996,6 +996,7 @@ public: int reformat_surf(const char* comp_name, LDBLE fraction, const char* new_comp_name, LDBLE new_Dw, int cell); LDBLE viscosity(void); + LDBLE calc_f_visc(const char *name); LDBLE calc_vm_Cl(void); int multi_D(LDBLE DDt, int mobile_cell, int stagnant); LDBLE find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant); diff --git a/phreeqcpp/SolutionIsotope.cxx b/phreeqcpp/SolutionIsotope.cxx index 4059b0b9..222d9110 100644 --- a/phreeqcpp/SolutionIsotope.cxx +++ b/phreeqcpp/SolutionIsotope.cxx @@ -67,12 +67,7 @@ cxxSolutionIsotope::dump_xml(std::ostream & s_oss, unsigned int indent) const s_oss << indent1; s_oss << "iso_ratio=\"" << this->ratio << "\"" << "\n"; -#ifdef NPP - if (!isnan(this->ratio_uncertainty)) -#else - //if (this->ratio_uncertainty != NAN) if (!std::isnan(this->ratio_uncertainty)) -#endif { s_oss << indent1; s_oss << "iso_ratio_uncertainty=\"" << this-> diff --git a/phreeqcpp/basicsubs.cpp b/phreeqcpp/basicsubs.cpp index 3094c18b..4e07a4e5 100644 --- a/phreeqcpp/basicsubs.cpp +++ b/phreeqcpp/basicsubs.cpp @@ -1252,6 +1252,24 @@ calc_t_sc(const char *name) return (-999.99); } +/* ---------------------------------------------------------------------- */ +LDBLE Phreeqc:: +calc_f_visc(const char *name) +/* ---------------------------------------------------------------------- */ +{ + char token[MAX_LENGTH]; + class species *s_ptr; + + if (print_viscosity) + { + strcpy(token, name); + s_ptr = s_search(token); + if (s_ptr != NULL) + return s_ptr->dw_t_visc; + } + return 0; +} + /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: equi_phase(const char *phase_name) @@ -1599,12 +1617,7 @@ get_calculate_value(const char *name) calculate_value_ptr->name); error_msg(error_string, STOP); } -#ifdef NPP - if (isnan(rate_moles)) -#else - //if (rate_moles == NAN) if(std::isnan(rate_moles)) -#endif { error_string = sformatf( "Calculated value not SAVEed for %s.", calculate_value_ptr->name); @@ -1762,18 +1775,17 @@ LDBLE Phreeqc:: pr_pressure(const char* phase_name) /* ---------------------------------------------------------------------- */ { - + int l; + class phase* phase_ptr = phase_bsearch(phase_name, &l, FALSE); + if (phase_ptr == NULL) + { + error_string = sformatf("Gas %s, not found.", phase_name); + warning_msg(error_string); + return (1e-99); + } cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); if (gas_phase_ptr != NULL) { - int l; - class phase* phase_ptr = phase_bsearch(phase_name, &l, FALSE); - if (phase_ptr == NULL) - { - error_string = sformatf("Gas %s, not found.", phase_name); - warning_msg(error_string); - 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]); @@ -1792,6 +1804,10 @@ pr_pressure(const char* phase_name) } } } + else if (phase_ptr->in != FALSE && phase_ptr->pr_in) + { + return phase_ptr->pr_p; + } return(0.0); } /* ---------------------------------------------------------------------- */ @@ -1807,35 +1823,35 @@ LDBLE Phreeqc:: pr_phi(const char *phase_name) /* ---------------------------------------------------------------------- */ { - cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); - if (gas_phase_ptr != NULL) - { int l; - class phase* phase_ptr = phase_bsearch(phase_name, &l, FALSE); + class phase* phase_ptr = phase_bsearch(phase_name, &l, FALSE); if (phase_ptr == NULL) { error_string = sformatf( "Gas %s, not found.", phase_name); warning_msg(error_string); return (1e-99); } - for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++) + cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); + if (gas_phase_ptr != NULL) { + 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); if (phase_ptr == phase_ptr_gas) { if (gas_phase_ptr->Get_pr_in()) - { - return phase_ptr->pr_phi; - } + return phase_ptr->pr_phi; else - { return gas_comp_ptr->Get_phi(); - } } } } + else if (phase_ptr->in != FALSE && phase_ptr->pr_in) + { + return phase_ptr->pr_phi; + } return (1.0); } /* ---------------------------------------------------------------------- */ diff --git a/phreeqcpp/cl1.cpp b/phreeqcpp/cl1.cpp index 4a3a8bbc..22a748e5 100644 --- a/phreeqcpp/cl1.cpp +++ b/phreeqcpp/cl1.cpp @@ -1,6 +1,5 @@ #include #include -//#include #include #include #include "Phreeqc.h" diff --git a/phreeqcpp/common/Utils.cxx b/phreeqcpp/common/Utils.cxx index 8495cd36..72fc8c2d 100644 --- a/phreeqcpp/common/Utils.cxx +++ b/phreeqcpp/common/Utils.cxx @@ -10,7 +10,6 @@ #include "Utils.h" #include "Parser.h" #include "float.h" -//#include #include #if defined(PHREEQCI_GUI) diff --git a/phreeqcpp/gases.cpp b/phreeqcpp/gases.cpp index 72d8fe12..9f71089f 100644 --- a/phreeqcpp/gases.cpp +++ b/phreeqcpp/gases.cpp @@ -607,6 +607,7 @@ calc_PR(void) phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) * log((rz + 2.41421356 * B) / (rz - 0.41421356 * B)); //phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi)); + phi = (phi > 4.44 ? 4.44 : (phi < -4.6 ? -4.6 : phi)); } else phi = -3.0; // fugacity coefficient = 0.05 diff --git a/phreeqcpp/global_structures.h b/phreeqcpp/global_structures.h index 053dbfb4..98f16eef 100644 --- a/phreeqcpp/global_structures.h +++ b/phreeqcpp/global_structures.h @@ -715,6 +715,7 @@ public: dw_a2 = 0; dw_a_visc = 0; // viscosity correction of SC dw_t_SC = 0; // contribution to SC, for calc'ng transport number with BASIC + dw_t_visc = 0; // contribution to viscosity dw_corr = 0; // dw corrected for TK and mu erm_ddl = 0; // enrichment factor in DDL equiv = 0; // equivalents in exchange species @@ -776,6 +777,7 @@ public: LDBLE dw_a2; LDBLE dw_a_visc; LDBLE dw_t_SC; + LDBLE dw_t_visc; LDBLE dw_corr; LDBLE erm_ddl; LDBLE equiv; diff --git a/phreeqcpp/inverse.cpp b/phreeqcpp/inverse.cpp index 7c0b63c9..72dffe40 100644 --- a/phreeqcpp/inverse.cpp +++ b/phreeqcpp/inverse.cpp @@ -3582,40 +3582,21 @@ check_isotopes(class inverse *inv_ptr) i = ii; /* use inverse-defined uncertainties first */ -#ifdef NPP - 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() && !std::isnan(inv_ptr->i_u[i].uncertainties[j])) -#endif { kit->second.Set_x_ratio_uncertainty(inv_ptr->i_u[i].uncertainties[j]); /* use solution-defined uncertainties second */ } -#ifdef NPP - 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 && !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]); /* use solution-defined uncertainties second */ } -#ifdef NPP - else if (!isnan(kit->second.Get_ratio_uncertainty())) -#else - //else if (kit->second.Get_ratio_uncertainty() != NAN) else if (!std::isnan(kit->second.Get_ratio_uncertainty())) -#endif { kit->second.Set_x_ratio_uncertainty( kit->second.Get_ratio_uncertainty()); @@ -3643,12 +3624,7 @@ 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 (std::isnan(kit->second.Get_x_ratio_uncertainty())) -#endif { error_string = sformatf( "In solution %d, isotope ratio uncertainty is needed for element: %g%s.", diff --git a/phreeqcpp/isotopes.cpp b/phreeqcpp/isotopes.cpp index 46e3c318..e09e43cb 100644 --- a/phreeqcpp/isotopes.cpp +++ b/phreeqcpp/isotopes.cpp @@ -900,12 +900,7 @@ punch_calculate_values(void) calculate_value_ptr->name); error_msg(error_string, STOP); } -#ifdef NPP - if (isnan(rate_moles)) -#else - //if (rate_moles == NAN) if (std::isnan(rate_moles)) -#endif { error_string = sformatf( "Calculated value not SAVEed for %s.", calculate_value_ptr->name); @@ -1134,12 +1129,7 @@ calculate_values(void) calculate_value_ptr->name); error_msg(error_string, STOP); } -#ifdef NPP - if (isnan(rate_moles)) -#else - //if (rate_moles == NAN) if (std::isnan(rate_moles)) -#endif { error_string = sformatf( "Calculated value not SAVEed for %s.", calculate_value_ptr->name); @@ -1202,12 +1192,7 @@ calculate_values(void) calculate_value_ptr->name); error_msg(error_string, STOP); } -#ifdef NPP - if (isnan(rate_moles)) -#else - //if (rate_moles == NAN) if (std::isnan(rate_moles)) -#endif { error_string = sformatf( "Calculated value not SAVEed for %s.", calculate_value_ptr->name); diff --git a/phreeqcpp/kinetics.cpp b/phreeqcpp/kinetics.cpp index 6e59bcb1..0472265c 100644 --- a/phreeqcpp/kinetics.cpp +++ b/phreeqcpp/kinetics.cpp @@ -115,12 +115,7 @@ calc_kinetic_reaction(cxxKinetics *kinetics_ptr, LDBLE time_step) kinetics_comp_ptr->Get_rate_name().c_str()); error_msg(error_string, STOP); } -#ifdef NPP - if (isnan(rate_moles)) -#else - //if (rate_moles == NAN) if (std::isnan(rate_moles)) -#endif { error_string = sformatf( "Moles of reaction not SAVEed for %s.", kinetics_comp_ptr->Get_rate_name().c_str()); diff --git a/phreeqcpp/prep.cpp b/phreeqcpp/prep.cpp index f98d584c..cc12fdc9 100644 --- a/phreeqcpp/prep.cpp +++ b/phreeqcpp/prep.cpp @@ -3997,6 +3997,7 @@ calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) * log((rz + 2.41421356 * B) / (rz - 0.41421356 * B)); //phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi)); + phi = (phi > 4.44 ? 4.44 : (phi < -4.6 ? -4.6 : phi)); //if (phi > 4.44) // phi = 4.44; } diff --git a/phreeqcpp/read.cpp b/phreeqcpp/read.cpp index 3052877d..9e875eff 100644 --- a/phreeqcpp/read.cpp +++ b/phreeqcpp/read.cpp @@ -2699,7 +2699,7 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v) /* * Read supcrt parms and Ionic strength terms */ - for (j = 0; j < 11; j++) + for (j = 0; j < 10; j++) { delta_v[j] = 0.0; } @@ -2707,7 +2707,7 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v) /* Vmax, dmax... delta_v[10] = 999.0; delta_v[11] = 1.0; */ - j = sscanf(cptr, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT , + j = sscanf(cptr, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT/* SCANFORMAT SCANFORMAT*/ , /* a1..a4 */ &(delta_v[0]), &(delta_v[1]), &(delta_v[2]), &(delta_v[3]), /* wref */ @@ -2715,9 +2715,9 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v) /* b_Av */ &(delta_v[5]), /* c1..c4 */ - &(delta_v[6]), &(delta_v[7]), &(delta_v[8]), &(delta_v[9]), + &(delta_v[6]), &(delta_v[7]), &(delta_v[8]), &(delta_v[9])/*, //vmP, exP - &(delta_v[10]), &(delta_v[11])); + &(delta_v[10]), &(delta_v[11])*/); if (j < 1) { input_error++; diff --git a/phreeqcpp/sundialsmath.cpp b/phreeqcpp/sundialsmath.cpp index 42487bb7..17143fa5 100644 --- a/phreeqcpp/sundialsmath.cpp +++ b/phreeqcpp/sundialsmath.cpp @@ -59,7 +59,6 @@ #include - //#include #include #include "sundialsmath.h" #include "sundialstypes.h" diff --git a/phreeqcpp/tidy.cpp b/phreeqcpp/tidy.cpp index db5574b1..2b038c22 100644 --- a/phreeqcpp/tidy.cpp +++ b/phreeqcpp/tidy.cpp @@ -979,12 +979,7 @@ tidy_gas_phase(void) error_msg(error_string, CONTINUE); } /* calculate moles */ -#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 (!std::isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read())) -#endif { P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read(); if (!PR) @@ -1013,12 +1008,7 @@ tidy_gas_phase(void) */ if (!gas_phase_ptr->Get_solution_equilibria()) { -#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 (!std::isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read())) -#endif { P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read(); if (!PR) @@ -1689,12 +1679,7 @@ tidy_ss_assemblage(void) phase_ptr->moles_x = 0; phase_ptr->fraction_x = 0; } -#ifdef NPP - if (isnan(comp_ptr->Get_moles())) -#else - //if (comp_ptr->Get_moles() == NAN) if (std::isnan(comp_ptr->Get_moles())) -#endif { input_error++; error_string = sformatf( @@ -3021,12 +3006,7 @@ tidy_isotopes(void) temp_iso.Set_total(0); temp_iso.Set_ratio(master[k]->isotope_ratio); temp_iso.Set_ratio_uncertainty(master[k]->isotope_ratio_uncertainty); -#ifdef NPP - if (!isnan(master[k]->isotope_ratio_uncertainty)) -#else - //if (master[k]->isotope_ratio_uncertainty != NAN) if (!std::isnan(master[k]->isotope_ratio_uncertainty)) -#endif { temp_iso.Set_ratio_uncertainty_defined(true); } diff --git a/phreeqcpp/transport.cpp b/phreeqcpp/transport.cpp index 378d818b..90f51054 100644 --- a/phreeqcpp/transport.cpp +++ b/phreeqcpp/transport.cpp @@ -5963,10 +5963,10 @@ viscosity(void) We use the harmonic mean of the Dw's, and the arithmetic mean of the z's, both weighted by the equivalent concentration. */ - LDBLE D1, D2, z1, z2, m_plus, m_min, eq_plus, eq_min, eq_dw_plus, eq_dw_min, t1, t2, ta; + 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, tc; - m_plus = m_min = eq_plus = eq_min = eq_dw_plus = eq_dw_min = V_an = m_an = V_Cl = ta = 0; + 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; @@ -5978,6 +5978,7 @@ viscosity(void) continue; 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 = 0; t1 = s_x[i]->moles / mass_water_aq_x; l_z = fabs(s_x[i]->z); if (l_z) @@ -5993,8 +5994,9 @@ viscosity(void) s_x[i]->Jones_Dole[8] / exp(-s_x[i]->Jones_Dole[4] * 25.0); } // find B * m and D * m * mu^d3 - Bc += (s_x[i]->Jones_Dole[0] + + s_x[i]->dw_t_visc = (s_x[i]->Jones_Dole[0] + s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * t1; + Bc += s_x[i]->dw_t_visc; // 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]; @@ -6002,8 +6004,10 @@ viscosity(void) t2 = -0.8 / s_x[i]->Jones_Dole[5]; else t2 = -1; - Dc += (s_x[i]->Jones_Dole[3] * exp(-s_x[i]->Jones_Dole[4] * tc)) * + t3 = (s_x[i]->Jones_Dole[3] * exp(-s_x[i]->Jones_Dole[4] * tc)) * t1 * (pow(mu_x, s_x[i]->Jones_Dole[5])*(1 + t2) + pow(t1 * f_z, s_x[i]->Jones_Dole[5])) / (2 + t2); + s_x[i]->dw_t_visc += t3; + Dc += t3; //output_msg(sformatf("\t%s\t%e\t%e\t%e\n", s_x[i]->name, t1, Bc, Dc )); } // parms for A... @@ -6023,13 +6027,11 @@ viscosity(void) { V_Cl = s_x[i]->logk[vm_tc]; V_an += V_Cl * s_x[i]->moles; - ta += s_x[i]->moles; m_an += s_x[i]->moles; } else// if (s_x[i]->Jones_Dole[6]) { V_an += s_x[i]->logk[vm_tc] * s_x[i]->Jones_Dole[6] * s_x[i]->moles; - ta += s_x[i]->moles; m_an += s_x[i]->moles; } if (Dw) @@ -6064,23 +6066,24 @@ viscosity(void) A = 0; viscos = viscos_0 + A * sqrt((eq_plus + eq_min) / 2 / mass_water_aq_x); if (m_an) - { V_an /= m_an; - ta /= m_an; - } if (!V_Cl) V_Cl = calc_vm_Cl(); - if (V_an && V_Cl && ta) - viscos += (viscos_0 * (2 - ta * V_an / V_Cl) * (Bc + Dc)); - else - viscos += (viscos_0 * (Bc + Dc)); + if (V_an) + fan = 2 - V_an / V_Cl; + //else + // fan = 1; + viscos += viscos_0 * fan * (Bc + Dc); if (viscos < 0) { viscos = viscos_0; warning_msg("viscosity < 0, reset to viscosity of pure water"); } + for (i = 0; i < (int)this->s_x.size(); i++) + { + s_x[i]->dw_t_visc /= (Bc + Dc); + } return viscos; - } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc::