diff --git a/phreeqcpp/PBasic.cpp b/phreeqcpp/PBasic.cpp
index 8fbbd67a..6f863c44 100644
--- a/phreeqcpp/PBasic.cpp
+++ b/phreeqcpp/PBasic.cpp
@@ -1612,6 +1612,9 @@ listtokens(FILE * f, tokenrec * l_buf)
  		case tokdh_a:
  			output_msg("DH_A"); // Debye-Hueckel A
  			break;
+		case tokdebye_length:
+			output_msg("DEBYE_LENGTH"); // Debye-Hueckel length
+			break;
  		case tokdh_b:
  			output_msg("DH_B"); // Debye-Hueckel B
  			break;
@@ -3706,6 +3709,13 @@ factor(struct LOC_exec * LINK)
 			n.UU.val = PhreeqcPtr->DH_A;
 		}
 		break;
+	case tokdebye_length:
+		{
+			double debye_length = (PhreeqcPtr->eps_r * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * PhreeqcPtr->tk_x)
+				/ (2. * F_C_MOL * F_C_MOL * PhreeqcPtr->mu_x * 1000.);
+			n.UU.val = sqrt(debye_length);
+			break;
+		}
 	case tokdh_b:
 		if (PhreeqcPtr->llnl_count_temp > 0)
 		{
@@ -7377,6 +7387,7 @@ const std::map<const std::string, PBasic::BASIC_TOKEN>::value_type temp_tokens[]
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("eps_r",              PBasic::tokeps_r),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("vm",                 PBasic::tokvm),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("dh_a",               PBasic::tokdh_a),
+	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("debye_length",       PBasic::tokdebye_length),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("dh_b",               PBasic::tokdh_b),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("dh_av",              PBasic::tokdh_av),
 	std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("qbrn",               PBasic::tokqbrn),
diff --git a/phreeqcpp/PBasic.h b/phreeqcpp/PBasic.h
index 76728e58..870e5850 100644
--- a/phreeqcpp/PBasic.h
+++ b/phreeqcpp/PBasic.h
@@ -317,6 +317,7 @@ public:
 		tokvm,
 		tokphase_vm,
 		tokdh_a,
+		tokdebye_length,
 		tokdh_b,
 		tokdh_av,
 		tokqbrn,
diff --git a/phreeqcpp/Phreeqc.h b/phreeqcpp/Phreeqc.h
index 254821f7..b1089dfd 100644
--- a/phreeqcpp/Phreeqc.h
+++ b/phreeqcpp/Phreeqc.h
@@ -1070,6 +1070,7 @@ public:
 	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);
+	void calc_b_ij(int icell, int jcell, int k, LDBLE b_i, LDBLE b_j, LDBLE g_i, LDBLE g_j, LDBLE free_i, LDBLE free_j, int stagnant);
 	void diffuse_implicit(LDBLE DDt, int stagnant);
 	int fill_spec(int cell_no, int ref_cell);
 	LDBLE moles_from_redox_states(cxxSolution *sptr, const char *name);
diff --git a/phreeqcpp/tidy.cpp b/phreeqcpp/tidy.cpp
index 2f1f0a94..510ab561 100644
--- a/phreeqcpp/tidy.cpp
+++ b/phreeqcpp/tidy.cpp
@@ -3613,10 +3613,10 @@ tidy_kin_exchange(void)
 			assert(false);
 		}
 		cxxExchange * exchange_ptr = &(it->second);
-		//if (!exchange_ptr->Get_new_def())
-		//	continue;
-		//if (exchange_ptr->Get_n_user() < 0)
-		//	continue;
+		if (!exchange_ptr->Get_new_def())
+			continue;
+		if (exchange_ptr->Get_n_user() < 0)
+			continue;
 		// check elements
 		for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++)
 		{
@@ -3846,10 +3846,10 @@ tidy_min_exchange(void)
 			assert(false);
 		}
 		cxxExchange * exchange_ptr = &(it->second);
-		//if (!exchange_ptr->Get_new_def())
-		//		continue;
-		//if (exchange_ptr->Get_n_user() < 0)
-		//	continue;
+		if (!exchange_ptr->Get_new_def())
+				continue;
+		if (exchange_ptr->Get_n_user() < 0)
+			continue;
 		n = exchange_ptr->Get_n_user();
 		// check elements
 		for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++)
@@ -4166,7 +4166,6 @@ tidy_min_surface(void)
 			assert(false);
 		}
 		cxxSurface *surface_ptr = &(kit->second);
-		if (!surface_ptr->Get_new_def()) continue;
 		if (!surface_ptr->Get_new_def())
 			continue;
 		if (surface_ptr->Get_n_user() < 0)
@@ -4414,7 +4413,20 @@ update_min_surface(void)
 			double comp_moles = 0.0;
 			cxxSurfaceComp* surface_comp_ptr = &(surface_ptr->Get_surface_comps()[j]);
 			if (surface_comp_ptr->Get_phase_name().size() == 0)	continue;
-			cxxSurfaceCharge* surface_charge_ptr = surface_ptr->Find_charge(surface_comp_ptr->Get_charge_name());
+			cxxSurfaceCharge* surface_charge_ptr =  NULL;
+			if (surface_ptr->Get_type() != cxxSurface::SURFACE_TYPE::NO_EDL)
+			{
+				surface_charge_ptr = surface_ptr->Find_charge(surface_comp_ptr->Get_charge_name());
+				if (surface_charge_ptr == NULL)
+				{
+					input_error++;
+					error_string = sformatf("Data structure for surface charge not found "
+						"for %s ",
+						surface_comp_ptr->Get_formula().c_str());
+					error_msg(error_string, CONTINUE);
+					continue;
+				}
+			}
 			int n = surface_ptr->Get_n_user();
 
 			/* First find surface master species */
@@ -4493,7 +4505,11 @@ update_min_surface(void)
 			surface_comp_ptr->Set_phase_name(phase_ptr->name);
 			/* make surface concentration proportional to mineral ... */
 			LDBLE conc = jit->second.Get_moles() * surface_comp_ptr->Get_phase_proportion();
-			double grams = surface_charge_ptr->Get_grams();
+			double grams = 0.0;
+			if (surface_charge_ptr != NULL)
+			{
+				grams = surface_charge_ptr->Get_grams();
+			}
 			if (comp_moles > 0.0)
 			{
 				surface_comp_ptr->multiply(conc / comp_moles);
@@ -4514,7 +4530,7 @@ update_min_surface(void)
 			{
 				surface_charge_ptr->multiply(jit->second.Get_moles() / grams);
 			}
-			else /* need to generate from scratch */
+			else if (surface_charge_ptr != NULL) /* need to generate from scratch */
 			{
 				surface_charge_ptr->Set_grams(jit->second.Get_moles());
 				surface_charge_ptr->Set_charge_balance(0.0);
@@ -4897,10 +4913,23 @@ update_kin_surface(void)
 			/* use database name for rate */
 			comp_ptr->Set_rate_name(kin_comp_ptr->Get_rate_name().c_str());
 			cxxSurfaceCharge* charge_ptr = surface_ptr->Find_charge(comp_ptr->Get_charge_name());
-
+			if (surface_ptr->Get_type() != cxxSurface::SURFACE_TYPE::NO_EDL)
+			{
+				charge_ptr = surface_ptr->Find_charge(comp_ptr->Get_charge_name());
+				if (charge_ptr == NULL)
+				{
+					input_error++;
+					error_string = sformatf("Data structure for surface charge not found "
+						"for %s ",
+						comp_ptr->Get_formula().c_str());
+					error_msg(error_string, CONTINUE);
+					continue;
+				}
+			}
 			/* make surface concentration proportional to mineral ... */
 			LDBLE conc = kin_comp_ptr->Get_m() * comp_ptr->Get_phase_proportion();
-			double grams = charge_ptr->Get_grams();
+			double grams = 0.0;
+			if (charge_ptr != NULL) charge_ptr->Get_grams();
 			if (comp_moles > 0.0)
 			{
 				comp_ptr->multiply(conc / comp_moles);
@@ -4922,7 +4951,7 @@ update_kin_surface(void)
 			{
 				charge_ptr->multiply(kin_comp_ptr->Get_m() / grams);
 			}
-			else /* need to generate from scratch */
+			else if (charge_ptr != NULL) /* need to generate from scratch */
 			{
 				charge_ptr->Set_grams(kin_comp_ptr->Get_m());
 				charge_ptr->Set_charge_balance(0.0);
diff --git a/phreeqcpp/transport.cpp b/phreeqcpp/transport.cpp
index 918a86a3..dffb6046 100644
--- a/phreeqcpp/transport.cpp
+++ b/phreeqcpp/transport.cpp
@@ -585,7 +585,6 @@ transport(void)
 						if (overall_iterations > max_iter)
 							max_iter = overall_iterations;
 						cell_no = i;
-						mixrun = j;
 						if (multi_Dflag)
 							sprintf(token,
 								"Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)",
@@ -746,7 +745,6 @@ transport(void)
 					if (i == first_c && count_cells > 1)
 						kin_time /= 2;
 					cell_no = i;
-					mixrun = 0;
 					if (multi_Dflag)
 						sprintf(token,
 							"Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)",
@@ -794,6 +792,7 @@ transport(void)
 				j = 1;
 			for (; j <= nmix; j++) // loop on j
 			{
+				mixrun = j;
 				if (multi_Dflag && j == nmix && (transport_step % print_modulus == 0))
 				{
 					sprintf(token,
@@ -850,7 +849,6 @@ transport(void)
 					if (overall_iterations > max_iter)
 						max_iter = overall_iterations;
 					cell_no = i;
-					mixrun = j;
 					if (multi_Dflag)
 						sprintf(token,
 							"Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)",
@@ -2827,7 +2825,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 	}
 	current_A = current_x / DDt * F_C_MOL;
 
-	for (i = ifirst; i <= ilast + stagnant + (bcon_last == 2 ? 1 : 0); i++)
+	for (i = ifirst; i <= ilast + stagnant + ((bcon_last == 2 || (dV_dcell && stagnant)) ? 1 : 0); i++)
 	{
 		if (i <= ilast + 1)
 		{
@@ -2879,10 +2877,12 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 		for (icell = if1; icell != il1; icell += incr)
 		{
 			min_mol = min_dif_M * ct[icell].kgw;
+			if (min_mol < 1e-13)
+				min_mol = 1e-13;
 			dum1 = dum2 = 0;
 			sptr1 = Utilities::Rxn_find(Rxn_solution_map, icell);
 			sptr2 = Utilities::Rxn_find(Rxn_solution_map, icell + 1);
-			if (stagnant)
+			if (stagnant && mixf_stag[icell][cp])
 			{
 				i1 = (icell == 0 ? c1 + 1 : icell == c1 ? cc1 : icell + c1);
 				sptr_stag = Utilities::Rxn_find(Rxn_solution_map, i1);
@@ -2890,6 +2890,13 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 			else
 				sptr_stag = NULL;
 
+			//if (!cp)
+			//{
+			//	ct[icell].J_ij_sum = ct[icell + 1].J_ij_sum = 0.0;
+			//	if (sptr_stag)
+			//		ct[i1].J_ij_sum = 0.0;
+			//}
+
 			if (!strcmp(ct[icell].m_s[cp].name, "H"))
 			{
 				dummy = ct[icell].m_s[cp].tot1;
@@ -2933,7 +2940,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 					}
 					sptr_stag->Set_total_o(sptr_stag->Get_total_o() + dummy);
 				}
-				//if (cp == count_m_s - 1) // transport the charge imbalance (but doesn't improve the change of H2 or O2).
+				//if (cp == count_m_s - 1) // transport the charge imbalance
 				//{
 				//	sptr1->Set_cb(sptr1->Get_cb() + ct[icell].J_ij_sum);
 				//	sptr2->Set_cb(sptr2->Get_cb() + ct[icell + 1].J_ij_sum);
@@ -2942,180 +2949,166 @@ diffuse_implicit(LDBLE DDt, int stagnant)
 				//}
 				continue;
 			}
-			// see if icell - incr has negative moles, subtract it, so that it is canceled...
-			if (icell > 0 && icell <= ilast &&
-				(it1 = neg_moles.find(icell - incr)) != neg_moles.end()
-				&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
-			{
-				ct[icell].m_s[cp].tot1 += it2->second;
-				neg_moles.erase(it1);
-				els.erase(it2);
-				neg_moles.insert(std::make_pair(icell - incr, els));
-			}
 			dum1 = sptr1->Get_totals()[ct[icell].m_s[cp].name];
-			if (!dum1)
+			dum2 = sptr2->Get_totals()[ct[icell].m_s[cp].name];
+			if (sptr_stag)
+				dum_stag = sptr_stag->Get_totals()[ct[icell].m_s[cp].name];
+			
+			// 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)))
 			{
 				dum1 = moles_from_redox_states(sptr1, ct[icell].m_s[cp].name);
-				if (dum1)
-					sptr1->Get_totals()[ct[icell].m_s[cp].name] = dum1;
-			}
-			dum2 = sptr2->Get_totals()[ct[icell].m_s[cp].name];
-			if (!dum2)
-			{
-				dum2 = moles_from_redox_states(sptr2, ct[icell].m_s[cp].name);
-				if (dum2)
-					sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2;
-			}
-			if (sptr_stag)
-			{
-				dum_stag = sptr_stag->Get_totals()[ct[icell].m_s[cp].name];
-				if (!dum_stag)
-				{
-					dum_stag = moles_from_redox_states(sptr_stag, ct[icell].m_s[cp].name);
-					if (dum_stag)
-						sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = dum_stag;
-				}
-			}
-
-			// check for negative moles, add moles from the donnan layer when necessary and available...
-			if (ct[icell].m_s[cp].tot1 > dum1 &&
-				(dV_dcell || (icell > 0 && icell <= ilast)))
-			{
-				if (!ct[icell].dl_s)
-					ct[icell].m_s[cp].tot1 = dum1;
-				else
+				if (ct[icell].dl_s > 1e-8)
 				{
 					cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell);
 					if (s_ptr)
 					{
-						dum1 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, ct[icell].m_s[cp].tot1 - dum1 + min_mol);
+						dum1 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, ct[icell].m_s[cp].tot1 + ct[icell].m_s[cp].tot_stag - dum1 + min_mol);
 					}
 				}
 				sptr1->Get_totals()[ct[icell].m_s[cp].name] = dum1;
+				if (dum1 - ct[icell].m_s[cp].tot1 - ct[icell].m_s[cp].tot_stag < min_mol)
+					ct[icell].m_s[cp].tot1 = dum1 - ct[icell].m_s[cp].tot_stag - min_mol;
 			}
 			// check for negative moles, in the other cell...
 			ct[icell].m_s[cp].tot2 = ct[icell].m_s[cp].tot1;
-			if (-ct[icell].m_s[cp].tot2 > dum2 &&
-				(dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)))
+			dum = 0;
+			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)*/))
 			{
-				if (!ct[icell + 1].dl_s)
-					ct[icell].m_s[cp].tot2 = -dum2;
-				else
+				dum2 = moles_from_redox_states(sptr2, ct[icell].m_s[cp].name);
+				if (ct[icell + 1].dl_s > 1e-8)
 				{
 					cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell + 1);
 					if (s_ptr)
 					{
 						dum2 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, -(ct[icell].m_s[cp].tot2 + dum2 - min_mol));
 					}
-					sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2;
-					if (-ct[icell].m_s[cp].tot2 > dum2)
-						ct[icell].m_s[cp].tot2 = -dum2;
 				}
+				sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2;
+				if (dum2 + ct[icell].m_s[cp].tot2 - dum < min_mol)
+					ct[icell].m_s[cp].tot2 = -dum2 + min_mol + dum;
 			}
 			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 (!cp)
-			{
-				ct[icell].J_ij_sum = ct[icell + 1].J_ij_sum = 0.0;
-				if (sptr_stag)
-					ct[i1].J_ij_sum = 0.0;
-			}
-
+			
 			if (dV_dcell || (icell > 0 && icell <= ilast))
 			{
-				dum1 -= ct[icell].m_s[cp].tot1;
-				if (sptr_stag)
-					dum1 -= ct[icell].m_s[cp].tot_stag;
-				sptr1->Get_totals()[ct[icell].m_s[cp].name] = (dum1 > 0 ? dum1 : 0e-16);
-				if (dum1 < -min_mol && incr > 0)
+				dum = ct[icell].m_s[cp].tot1;
+				if (stagnant)
+					dum += ct[icell].m_s[cp].tot_stag;
+				dum1 -= dum;
+				sptr1->Get_totals()[ct[icell].m_s[cp].name] = (dum1 > 0 ? dum1 : min_mol);
+				if (dum1 < 0)
 				{
+					dum += dum1 - min_mol;
+					if ((it1 = neg_moles.find(icell)) != neg_moles.end()
+						&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
+						dum1 += it2->second;
 					els.clear();
 					els.insert(std::make_pair(ct[icell].m_s[cp].name, dum1));
 					neg_moles.erase(icell);
 					neg_moles.insert(std::make_pair(icell, els));
 				}
-				else
-					ct[icell].J_ij_sum += dum1 * ct[icell].m_s[cp].charge;
+				//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 || (icell >= 0 && icell < ilast)/* || (icell == ilast && bcon_last == 2)*/)
 			{
-				dum2 += ct[icell].m_s[cp].tot1;
-				sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : 0e-16);
-				if (dum2 < -min_mol && incr < 0)
+				dum = ct[icell].m_s[cp].tot1;
+				if (stagnant && icell == c && dV_dcell)
+					dum -= ct[c1].m_s[cp].tot_stag;
+				dum2 += dum;
+				sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : min_mol);
+				if (dum2 < 0)
 				{
+					dum -= dum2 - min_mol;
+					if ((it1 = neg_moles.find(icell + 1)) != neg_moles.end()
+						&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
+						dum2 += it2->second;
 					els.clear();
 					els.insert(std::make_pair(ct[icell].m_s[cp].name, dum2));
 					neg_moles.erase(icell + 1);
 					neg_moles.insert(std::make_pair(icell + 1, els));
 				}
-				else
-					ct[icell + 1].J_ij_sum += dum2 * ct[icell].m_s[cp].charge;
+				//ct[icell + 1].J_ij_sum += dum * ct[icell].m_s[cp].charge;
 			}
-
+			
 			if (sptr_stag)
 			{
-				if ((it1 = neg_moles.find(i1)) != neg_moles.end()
-					&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
-				{
-					ct[icell].m_s[cp].tot_stag += it2->second;
-					// the negative moles are canceled...
-					neg_moles.erase(it1);
-					els.erase(it2);
-					neg_moles.insert(std::make_pair(i1, els));
-				}
-				dummy = ct[icell].m_s[cp].tot_stag;
+				dum = ct[icell].m_s[cp].tot_stag;
 				if (icell == c)
+					dum += ct[c1].m_s[cp].tot_stag;
+				if (dum_stag + dum < 0)
 				{
-					if (dV_dcell)
-					{
-						if ((dum = sptr2->Get_totals()[ct[icell + 1].m_s[cp].name] - ct[icell + 1].m_s[cp].tot_stag) > 0)
-						{
-							sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum;
-							dummy += ct[icell + 1].m_s[cp].tot_stag;
-						}
-						else
-						{
-							dummy += sptr2->Get_totals()[ct[icell].m_s[cp].name];
-							sptr2->Get_totals()[ct[icell].m_s[cp].name] = 0.0;
-							if (dum < -min_mol)
-							{
-								els.clear();
-								els.insert(std::make_pair(ct[icell].m_s[cp].name, dum));
-								neg_moles.erase(icell);
-								neg_moles.insert(std::make_pair(icell, els));
-							}
-						}
-					}
-					else
-						dummy += ct[icell + 1].m_s[cp].tot_stag;
-				}
-				if (-dummy > dum_stag)
-				{
+					dum_stag = moles_from_redox_states(sptr_stag, ct[icell].m_s[cp].name);
 					if (ct[i1].dl_s)
 					{
 						cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, i1);
 						if (s_ptr)
 						{
-							dum_stag += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, -(dummy + dum_stag - min_mol));
+							dum_stag += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, -(dum + dum_stag - min_mol));
 						}
 						sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = dum_stag;
 					}
 				}
-				dum_stag += dummy;
-				sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = (dum_stag > 0 ? dum_stag : 0e-16);
-				if (dum_stag < -min_mol)
+				dum_stag += dum;
+				sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = (dum_stag > 0 ? dum_stag : min_mol);
+				if (dum_stag < 0)
 				{
+					dum -= dum_stag - min_mol;
+					if ((it1 = neg_moles.find(i1)) != neg_moles.end()
+						&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
+						dum_stag += it2->second;
 					els.clear();
 					els.insert(std::make_pair(ct[icell].m_s[cp].name, dum_stag));
 					neg_moles.erase(i1);
 					neg_moles.insert(std::make_pair(i1, els));
 				}
-				else
-					ct[i1].J_ij_sum += dum_stag * ct[icell].m_s[cp].charge;
+				//ct[i1].J_ij_sum += dum * ct[icell].m_s[cp].charge;
 			}
-			//if (cp == count_m_s - 1)
+			
+			// reduce oscillations in the column-boundary cells, but not for H and O, and current_A is not adjusted...
+			if (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)
+				{
+					sptr1 = Utilities::Rxn_find(Rxn_solution_map, 1);
+					sptr2 = Utilities::Rxn_find(Rxn_solution_map, 2);
+					dum1 = sptr1->Get_totals()[ct[0].m_s[cp].name] / ct[1].kgw * (1 - ct[1].dl_s) - dummy;
+					dum2 = sptr2->Get_totals()[ct[0].m_s[cp].name] / ct[2].kgw * (1 - ct[2].dl_s) - dummy;
+					if (dum1 / dum2 < 0 || dum1 / dum2 > 1)
+					{
+						dum = cell_data[1].mid_cell_x / cell_data[2].mid_cell_x;
+						//ct[1].J_ij_sum -= sptr1->Get_totals()[ct[0].m_s[cp].name] * ct[1].m_s[cp].charge;
+						dum1 = (dummy + dum * dum2) * ct[1].kgw / (1 - ct[1].dl_s);
+						sptr1->Get_totals()[ct[0].m_s[cp].name] = dum1;
+						//ct[1].J_ij_sum += dum1 * ct[1].m_s[cp].charge;
+					}
+				}
+				dummy = Utilities::Rxn_find(Rxn_solution_map, c1)->Get_totals()[ct[0].m_s[cp].name] / ct[c1].kgw  * (1 - ct[c1].dl_s);
+				if (dummy > 1e-6)
+				{
+					sptr1 = Utilities::Rxn_find(Rxn_solution_map, c);
+					sptr2 = Utilities::Rxn_find(Rxn_solution_map, c_1);
+					dum1 = sptr1->Get_totals()[ct[0].m_s[cp].name] / ct[c].kgw   * (1 - ct[c].dl_s)   - dummy;
+					dum2 = sptr2->Get_totals()[ct[0].m_s[cp].name] / ct[c_1].kgw * (1 - ct[c_1].dl_s) - dummy;
+					if (dum1 / dum2 < 0 || dum1 / dum2 > 1)
+					{
+						dum = (cell_data[c].mid_cell_x - cell_data[c_1].mid_cell_x) /
+							(cell_data[c1].mid_cell_x - cell_data[c_1].mid_cell_x);
+						//ct[c].J_ij_sum -= sptr1->Get_totals()[ct[0].m_s[cp].name] * ct[c].m_s[cp].charge;
+						dum1 = (dummy + (1 - dum) * dum2) * ct[c].kgw / (1 - ct[c].dl_s);
+						sptr1->Get_totals()[ct[0].m_s[cp].name] = dum1;
+						//ct[c].J_ij_sum += dum1 * ct[c].m_s[cp].charge;
+					}
+				}
+			}
+//if (cp == count_m_s - 1)
 			//{
 			//	sptr1->Set_cb(sptr1->Get_cb() + ct[icell].J_ij_sum);
 			//	sptr2->Set_cb(sptr2->Get_cb() + ct[icell + 1].J_ij_sum);
@@ -3664,6 +3657,46 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell, int stagnant)
 	}
 	return (OK);
 }
+/* ---------------------------------------------------------------------- */
+void Phreeqc::
+calc_b_ij(int icell, int jcell, int k, LDBLE b_i, LDBLE b_j, LDBLE g_i, LDBLE g_j, LDBLE free_i, LDBLE free_j, int stagnant)
+/* ---------------------------------------------------------------------- */ 
+{
+	ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) * b_j * (free_j + g_j) / (b_i * (free_i + g_i) + b_j * (free_j + g_j));
+	// At filterends, concentrations of ions change step-wise to the DL. 
+	// We take the harmonic mean for f_free, the average for the DL.
+	if (ct[icell].v_m[k].z)
+	{
+		if (!g_i && g_j)
+		{
+			ct[icell].v_m[k].b_ij = free_j * b_i * b_j / (b_i + b_j) +
+				b_i * (1 - free_j) / 4 + b_j * g_j / 4;
+		}
+		else if (g_i && !g_j)
+			ct[icell].v_m[k].b_ij = free_i * b_i * b_j / (b_i + b_j) +
+			b_j * (1 - free_i) / 4 + b_i * g_i / 4;
+	}
+	// for boundary cells...
+	if (stagnant > 1)
+	{ /* for a diffusion experiment with well-mixed reservoir in cell 3 and the last stagnant cell,
+		   and with the mixf * 2 for the boundary cells in the input... */
+		if (icell == 3 && !g_i && g_j)
+			ct[icell].v_m[k].b_ij = b_j * (free_j + g_j) / 2;
+		else if (jcell == all_cells - 1 && !g_j && g_i)
+			ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) / 2;
+	}
+	else
+	{
+		if (icell == 0 || (icell == count_cells + 1 && jcell == count_cells + count_cells + 1))
+			ct[icell].v_m[k].b_ij = b_j * (free_j + g_j);
+		else if (icell == count_cells && jcell == count_cells + 1)
+			ct[icell].v_m[k].b_ij = b_i * (free_i + g_i);
+	}
+	if (ct[icell].v_m[k].z)
+		ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
+	return;
+}
+
 /* ---------------------------------------------------------------------- */
 LDBLE Phreeqc::
 find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
@@ -3848,7 +3881,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
 	if (dl_aq1 > 0)
 		ct[icell].dl_s = dl_aq1 / t_aq1;
 	if (dl_aq2 > 0)
-		ct[icell].dl_s = ct[jcell].dl_s = dl_aq2 / t_aq2;
+	{
+		ct[jcell].dl_s = dl_aq2 / t_aq2;
+		if (!ct[icell].dl_s)
+			ct[icell].dl_s = 1e-8; // used in implicit appt
+	}
 
 	if (il_calcs)
 	{
@@ -4142,7 +4179,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
 									g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
 								else
 								{
-									dum1 = it_sc->Get_mass_water() / mass_water_bulk_x;
+									dum1 = it_sc->Get_mass_water() / t_aq2;
 									dum2 = it_sc->Get_z_gMCD_map()[1] / dum1;
 									g_j += pow(dum2, ct[icell].v_m[k].z) * dum1;
 								}
@@ -4152,32 +4189,18 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
 					}
 				}
 
-				b_i = A1 * sol_D[icell].spec[i].Dwt * (f_free_i + g_i / ct[icell].visc1);
-				b_j = A2 * (f_free_j + g_j / ct[icell].visc2);
-				if (icell == count_cells && !stagnant)
-					ct[icell].v_m[k].b_ij = b_i;
-				else if (icell == all_cells - 1 && stagnant)
-					ct[icell].v_m[k].b_ij = b_i / 2; /* with the mixf *= 2 for this 'reservoir' cell in the input */
+				b_i = A1 * sol_D[icell].spec[i].Dwt;
+				b_j = A2;
+				if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
+					b_j *= sol_D[icell].spec[i].Dwt;
 				else
 				{
-					if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
-						b_j *= sol_D[icell].spec[i].Dwt;
-					else
-					{
-						dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f;
-						dum2 *= exp(sol_D[icell].spec[i].dw_t / sol_D[jcell].tk_x - sol_D[icell].spec[i].dw_t / sol_D[icell].tk_x);
-						dum2 *= sol_D[jcell].viscos_f;
-						b_j *= dum2;
-					}
-					ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
-					if (icell == 0 && !stagnant)
-							ct[icell].v_m[k].b_ij = b_j;
-					else if (icell == 3 && stagnant && !g_i && g_j)
-						ct[icell].v_m[k].b_ij = b_j / 2; /* with the mixf *= 2 for stagnant cell 3 in the input */
+					dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f;
+					dum2 *= exp(sol_D[icell].spec[i].dw_t / sol_D[jcell].tk_x - sol_D[icell].spec[i].dw_t / sol_D[icell].tk_x);
+					dum2 *= sol_D[jcell].viscos_f;
+					b_j *= dum2;
 				}
-
-				if (ct[icell].v_m[k].z)
-					ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
+				calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
 
 				k++;
 			}
@@ -4249,7 +4272,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
 									g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
 								else
 								{
-									dum1 = it_sc->Get_mass_water() / mass_water_bulk_x;
+									dum1 = it_sc->Get_mass_water() / t_aq1;
 									dum2 = it_sc->Get_z_gMCD_map()[1] / dum1;
 									g_i += pow(dum2, ct[icell].v_m[k].z) * dum1;
 								}
@@ -4266,31 +4289,18 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
 						g_j *= sol_D[jcell].spec[j].erm_ddl;
 					}
 				}
-				b_i = A1 * (f_free_i + g_i / ct[icell].visc1);
-				b_j = A2 * sol_D[jcell].spec[j].Dwt * (f_free_j + g_j / ct[icell].visc2);
-				if (icell == 0 && !stagnant)
-					ct[icell].v_m[k].b_ij = b_j;
-				else if (icell == 3 && stagnant && g_j && !g_i)
-					ct[icell].v_m[k].b_ij = b_j / 2; /* with the mixf *= 2 for 'reservoir' cell 3 in the input */
+				b_i = A1;
+				b_j = A2 * sol_D[jcell].spec[j].Dwt;
+				if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
+					b_i *= sol_D[jcell].spec[j].Dwt;
 				else
 				{
-					if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
-						b_i *= sol_D[jcell].spec[j].Dwt;
-					else
-					{
-						dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f;
-						dum2 *= exp(sol_D[jcell].spec[j].dw_t / sol_D[icell].tk_x - sol_D[jcell].spec[j].dw_t / sol_D[jcell].tk_x);
-						dum2 *= sol_D[icell].viscos_f;
-						b_i *= dum2;
-					}
-					ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
-					if (icell == count_cells && !stagnant)
-						ct[icell].v_m[k].b_ij = b_i;
-					else if (jcell == all_cells - 1 && stagnant && !g_j && g_i)
-						ct[icell].v_m[k].b_ij = b_i / 2; /* with the mixf * 2 for this 'reservoir' cell in the input */
+					dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f;
+					dum2 *= exp(sol_D[jcell].spec[j].dw_t / sol_D[icell].tk_x - sol_D[jcell].spec[j].dw_t / sol_D[jcell].tk_x);
+					dum2 *= sol_D[icell].viscos_f;
+					b_i *= dum2;
 				}
-				if (ct[icell].v_m[k].z)
-					ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
+				calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
 
 				k++;
 			}
@@ -4373,28 +4383,9 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
 						g_j *= sol_D[jcell].spec[j].erm_ddl;
 					}
 				}
-				b_i = A1 * sol_D[icell].spec[i].Dwt * (f_free_i + g_i / ct[icell].visc1);
-				b_j = A2 * sol_D[jcell].spec[j].Dwt * (f_free_j + g_j / ct[icell].visc2);
-				ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
-				// but for boundary cells...
-				if (stagnant > 1)
-				{ /* for a diffusion experiment with well-mixed reservoir in cell 3 and the last stagnant cell,
-					   and with the mixf * 2 for the boundary cells in the input... */
-					if (icell == 3 && !g_i && g_j)
-						ct[icell].v_m[k].b_ij = b_j / 2;
-					else if (jcell == all_cells - 1 && !g_j && g_i)
-						ct[icell].v_m[k].b_ij = b_i / 2;
-				}
-				else
-				{
-					if (icell == 0 || (icell == count_cells + 1 && jcell == count_cells + count_cells + 1))
-						ct[icell].v_m[k].b_ij = b_j;
-					else if (icell == count_cells && jcell == count_cells + 1)
-						ct[icell].v_m[k].b_ij = b_i;
-				}
-
-				if (ct[icell].v_m[k].z)
-					ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
+				b_i = A1 * sol_D[icell].spec[i].Dwt;
+				b_j = A2 * sol_D[jcell].spec[j].Dwt;
+				calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
 
 				//ddlm = sol_D[jcell].spec[j].lm - sol_D[icell].spec[i].lm; // appt: this could give an incorrect large factor for implicit
 				//if (fabs(ddlm) > 1e-10)