Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 15351)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 15352)
@@ -8,17 +8,17 @@
 #include "../modules/modules.h"
 
 void solutionsequence_hydro_nonlinear(FemModel* femmodel){
-
 	/*solution : */
 	Vector<IssmDouble>* ug_sed=NULL;
 	Vector<IssmDouble>* ug_epl=NULL; 
 	Vector<IssmDouble>* uf=NULL;
-	Vector<IssmDouble>* old_uf=NULL; 
-	Vector<IssmDouble>* aged_ug_sed=NULL; 
-	Vector<IssmDouble>* aged_ug_epl=NULL; 
+	Vector<IssmDouble>* uf_int_iter=NULL; 
+	Vector<IssmDouble>* ug_sed_main_iter=NULL; 
+	Vector<IssmDouble>* ug_epl_main_iter=NULL; 
 	Vector<IssmDouble>* ys=NULL; 
 	Vector<IssmDouble>* dug=NULL; 
-
+	Vector<IssmDouble>* old_ug=NULL; 
+	
 	Matrix<IssmDouble>* Kff=NULL;
 	Matrix<IssmDouble>* Kfs=NULL;
 	Vector<IssmDouble>* pf=NULL;
@@ -33,14 +33,14 @@
 	IssmDouble sediment_kmax,time;
 	IssmDouble eps_hyd;
 	IssmDouble ndu_sed,nu_sed;
-	IssmDouble ndu_epl,nu_epl;	
+	IssmDouble ndu_epl,nu_epl;
 
 	/*Recover parameters: */
 	femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);//FIXME
 	femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
 	femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum);
 	femmodel->parameters->FindParam(&time,TimeEnum);
-	hydro_maxiter=100;
+	hydro_maxiter=150;
 	hydrocount=1;
 
 	/*Iteration on the two layers*/
@@ -56,11 +56,11 @@
 		sedcount=1;
 		eplcount=1;
 		//save pointer to old velocity
-		delete aged_ug_sed;
-		aged_ug_sed=ug_sed;
+		delete ug_sed_main_iter;
+		ug_sed_main_iter=ug_sed;
 		if(isefficientlayer){
-			delete aged_ug_epl;
-			aged_ug_epl=ug_epl;
+			delete ug_epl_main_iter;
+			ug_epl_main_iter=ug_epl;
 		}
 
 		femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
@@ -77,9 +77,9 @@
 			CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum);
 			Reduceloadx(pf,Kfs,ys); delete Kfs;
 			delete uf;
-			Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters);
-			delete old_uf; old_uf=uf->Duplicate();
-			if(sedcount>1)delete ug_sed; /*Not on first time to avoid deleting aged_ug_sed*/
+			Solverx(&uf, Kff, pf,uf_int_iter, df, femmodel->parameters);
+			delete uf_int_iter; uf_int_iter=uf->Duplicate();
+			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;
@@ -119,15 +119,15 @@
 				CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
 				Reduceloadx(pf,Kfs,ys); delete Kfs;
 				delete uf;
-				Solverx(&uf, Kff, pf,old_uf, df, femmodel->parameters);
-				delete old_uf; old_uf=uf->Duplicate();
+				Solverx(&uf, Kff, pf,uf_int_iter, df, femmodel->parameters);
+				delete uf_int_iter; uf_int_iter=uf->Duplicate();
 				if(eplcount>1) delete ug_epl; 
 				delete Kff;delete pf;
 				delete df;
 				Mergesolutionfromftogx(&ug_epl,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
 				InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
-				//femmodel->HydrologyEPLupdateDomainx();
 				ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+				femmodel->HydrologyEPLupdateDomainx();			
 				
 				if (!eplconverged){
 					if(VerboseConvergence()) _printf0_("   # EPL unstable constraints = " << num_unstable_constraints << "\n");
@@ -137,12 +137,11 @@
 					}
 				}
 				eplcount++;
-				
+
 				if(eplconverged){
 					InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,eplconverged,ConvergedEnum);
 					InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum);
 					InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug_epl);
-					//femmodel->HydrologyEPLupdateDomainx();
 					break;
 				}
 			}
@@ -151,14 +150,14 @@
 		if(!hydroconverged){
 			//compute norm(du)/norm(u)
 			dug=ug_sed->Duplicate(); _assert_(dug);
-			aged_ug_sed->Copy(dug);	
+			ug_sed_main_iter->Copy(dug);	
 			dug->AYPX(ug_sed,-1.0);
-			ndu_sed=dug->Norm(NORM_TWO); nu_sed=aged_ug_sed->Norm(NORM_TWO);
+			ndu_sed=dug->Norm(NORM_TWO); nu_sed=ug_sed_main_iter->Norm(NORM_TWO);
 			if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("Sed convergence criterion is NaN!");
 			if (!xIsNan<IssmDouble>(eps_hyd)){
 				if (!isefficientlayer){
 					if ((ndu_sed/nu_sed)<eps_hyd){
-						if(VerboseConvergence()) _printf0_(setw(50) << left << "   Converged \n");
+						if(VerboseConvergence()) _printf0_(setw(50) << left << "   Converged after, " << hydrocount << " iterations \n");
 						hydroconverged=true;
 					}
 					else{ 
@@ -168,15 +167,15 @@
 				}
 				else{
 					dug=ug_epl->Duplicate();_assert_(dug); 
-					aged_ug_epl->Copy(dug);_assert_(aged_ug_epl); 
+					ug_epl_main_iter->Copy(dug);_assert_(ug_epl_main_iter); 
 					dug->AYPX(ug_epl,-1.0);
 					ndu_epl=dug->Norm(NORM_TWO); 
-					nu_epl=aged_ug_epl->Norm(NORM_TWO);
+					nu_epl=ug_epl_main_iter->Norm(NORM_TWO);
 
 					if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("EPL convergence criterion is NaN!");
 					if (ndu_epl==0.0 && nu_epl==0.0) nu_epl=1.0e-6; /*Hacking the case where the EPL is used but empty*/
-					if ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<eps_hyd){
-						if (VerboseConvergence()) _printf0_(setw(50) << left << "   Converged \n");
+					if ((ndu_epl/nu_epl)<eps_hyd && (ndu_sed/nu_sed)<(eps_hyd*10)){
+						if (VerboseConvergence()) _printf0_(setw(50) << left << "   Converged after, " << hydrocount << " iterations \n");
 						hydroconverged=true;
 					}
 					else{ 
@@ -201,9 +200,8 @@
 	delete ug_epl;
 	delete ug_sed;
 	delete uf;
-	delete old_uf;
-	delete aged_ug_sed;
-	delete aged_ug_epl;
+	delete uf_int_iter;
+	delete ug_sed_main_iter;
+	delete ug_epl_main_iter;
 	delete dug;
-
 }
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 15351)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 15352)
@@ -852,6 +852,7 @@
 			element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
 			element->CreateKMatrix(Kff_temp,NULL);
 		}
+
 		for (i=0;i<this->loads->Size();i++){
 			load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
 			if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff_temp,NULL);
@@ -871,6 +872,7 @@
 		element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
 		element->CreateKMatrix(Kff,Kfs);
 	}
+	
 	for (i=0;i<this->loads->Size();i++){
 		load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
 		if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs);
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15351)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15352)
@@ -6131,7 +6131,6 @@
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
-
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 
 		/*Diffusivity*/
@@ -6156,7 +6155,6 @@
 						&Ke->values[0],1);
 		}
 	}
-
 	/*Clean up and return*/
 	delete gauss;
 	return Ke;
@@ -6243,11 +6241,9 @@
 
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussTria(2);
-
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
-
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis, gauss);
 
@@ -6310,7 +6306,6 @@
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
-
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis, gauss);
 
@@ -6690,7 +6685,7 @@
 
 		/*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
 		this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
-		if(eplhead[i]>=h_max && this->AnyActive()){			
+		if(eplhead[i]>=h_max && this->AnyActive()){
 			for(j=0;j<numdof;j++){
 				if(old_active[j]>0.){
 					vec_mask->SetValue(nodes[i]->Sid(),1.,INS_VAL);