Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15406)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15407)
@@ -9336,7 +9336,9 @@
 			if(values[i]>h_max){
 				residual[i]=kappa*(values[i]-h_max);
+				residual[i+numdof2d]=residual[i];
 			}
 			else{
 				residual[i]=0.0;
+				residual[i+numdof2d]=residual[i];
 			}
 		}
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15406)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15407)
@@ -6036,9 +6036,9 @@
 			}
 		}
+		this->inputs->AddInput(new TriaInput(SedimentHeadResidualEnum,residual,P1Enum));
 	}
 
 	/*Add input to the element: */
 	this->inputs->AddInput(new TriaInput(SedimentHeadEnum,values,P1Enum));
-	this->inputs->AddInput(new TriaInput(SedimentHeadResidualEnum,residual,P1Enum));
 
 	/*Free ressources:*/
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 15406)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 15407)
@@ -10,13 +10,16 @@
 void solutionsequence_hydro_nonlinear(FemModel* femmodel){
 	/*solution : */
-	Vector<IssmDouble>* ug_sed=NULL;
+	Vector<IssmDouble>* ug_sed=NULL; 
+	Vector<IssmDouble>* uf_sed=NULL; 
+	Vector<IssmDouble>* uf_sed_sub_iter=NULL; 
+	Vector<IssmDouble>* ug_sed_main_iter=NULL;
+
 	Vector<IssmDouble>* ug_epl=NULL; 
-	Vector<IssmDouble>* uf=NULL;
-	Vector<IssmDouble>* uf_int_iter=NULL; 
-	Vector<IssmDouble>* ug_sed_main_iter=NULL; 
-	Vector<IssmDouble>* ug_epl_main_iter=NULL; 
+	Vector<IssmDouble>* uf_epl=NULL;
+	Vector<IssmDouble>* uf_epl_sub_iter=NULL; 
+	Vector<IssmDouble>* ug_epl_main_iter=NULL;
+ 
 	Vector<IssmDouble>* ys=NULL; 
-	Vector<IssmDouble>* dug=NULL; 
-	Vector<IssmDouble>* old_ug=NULL; 
+	Vector<IssmDouble>* dug=NULL;
 	
 	Matrix<IssmDouble>* Kff=NULL;
@@ -47,7 +50,9 @@
 	femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
 	GetSolutionFromInputsx(&ug_sed, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
+	Reducevectorgtofx(&uf_sed, ug_sed, femmodel->nodes,femmodel->parameters);
 	if(isefficientlayer) {
 		femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
 		GetSolutionFromInputsx(&ug_epl, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
+	Reducevectorgtofx(&uf_epl, ug_epl, femmodel->nodes,femmodel->parameters);
 	}
 
@@ -58,8 +63,12 @@
 		//save pointer to old velocity
 		delete ug_sed_main_iter;
+		delete uf_sed_sub_iter;
 		ug_sed_main_iter=ug_sed;
+		uf_sed_sub_iter=uf_sed;
 		if(isefficientlayer){
 			delete ug_epl_main_iter;
+			delete uf_sed_sub_iter;
 			ug_epl_main_iter=ug_epl;
+			uf_sed_sub_iter=uf_epl;
 		}
 
@@ -74,14 +83,16 @@
 		for(;;){
 			femmodel->HydrologyTransferx();
-			femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &sediment_kmax);
+			femmodel->SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax);
 			CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum);
 			Reduceloadx(pf,Kfs,ys); delete Kfs;
-			delete uf;
-			Solverx(&uf, Kff, pf,uf_int_iter, df, femmodel->parameters);
-			delete uf_int_iter; uf_int_iter=uf->Duplicate();
+			if(sedcount>1)delete uf_sed;
+			Solverx(&uf_sed,Kff,pf,uf_sed_sub_iter,df,femmodel->parameters);
+			delete uf_sed_sub_iter; 
+			uf_sed_sub_iter=uf_sed->Duplicate();
+			uf_sed->Copy(uf_sed_sub_iter);
 			if(sedcount>1)delete ug_sed; /*Not on first time to avoid deleting ug_sed_main_iter*/
 			delete Kff; delete pf; delete df;
 
-			Mergesolutionfromftogx(&ug_sed,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
+			Mergesolutionfromftogx(&ug_sed,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys;
 			InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_sed);
 			ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
@@ -116,14 +127,16 @@
 			for(;;){
 				femmodel->HydrologyTransferx();
-				femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df,NULL);
+				femmodel->SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL);
 				CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
 				Reduceloadx(pf,Kfs,ys); delete Kfs;
-				delete uf;
-				Solverx(&uf, Kff, pf,uf_int_iter, df, femmodel->parameters);
-				delete uf_int_iter; uf_int_iter=uf->Duplicate();
+				if(eplcount>1) delete uf_epl;
+				Solverx(&uf_epl,Kff,pf,uf_epl_sub_iter,df,femmodel->parameters);
+				delete uf_epl_sub_iter;
+				uf_epl_sub_iter=uf_epl->Duplicate();
+				uf_epl->Copy(uf_epl_sub_iter);
 				if(eplcount>1) delete ug_epl; 
 				delete Kff;delete pf;
 				delete df;
-				Mergesolutionfromftogx(&ug_epl,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
+				Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys;
 				InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
 				ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
@@ -200,6 +213,8 @@
 	delete ug_epl;
 	delete ug_sed;
-	delete uf;
-	delete uf_int_iter;
+	delete uf_sed;
+	//	delete uf_epl;
+	delete uf_sed_sub_iter;
+	delete uf_epl_sub_iter;
 	delete ug_sed_main_iter;
 	delete ug_epl_main_iter;