Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 19090)
+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 19091)
@@ -645,7 +645,7 @@
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
 
-void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element){
+void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, int* eplzigzag_counter, Element* element){
 
 	bool        active_element;
 	int         i,j;
@@ -671,6 +671,8 @@
 
 	int         eplflip_lock;
 	int         numnodes      =basalelement->GetNumberOfNodes();
+	IssmDouble  new_active;
+	IssmDouble  recurent;
 	IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
 	IssmDouble* old_active    =xNew<IssmDouble>(numnodes);
 	IssmDouble* sedhead       =xNew<IssmDouble>(numnodes);
@@ -696,44 +698,55 @@
 	for(i=1;i<numnodes;i++) if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];
 
 	for(i=0;i<numnodes;i++){
-		/*Activate EPL if residual is >0 */
+
+		/*Dealing with the initial value to define zigzaging, preceding iteration the first time we see the node and current evolving vector after*/
+		recurence->GetValue(&recurent,basalelement->nodes[i]->Sid());
+		if (recurent>0)vec_mask->GetValue(&old_active[i],basalelement->nodes[i]->Sid());
+		new_active=old_active[i];
+		recurence->SetValue(basalelement->nodes[i]->Sid(),1.,ADD_VAL);
+
+		/*Now starting to look at the activations*/
 		if(residual[i]>0.){
-			if(old_active[i]==0) eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
 			vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
+			new_active=1.0;
 		}
-
 		/*If mask was already one, keep one or colapse*/
 		else if(old_active[i]>0.){
-			
 			vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
+			new_active=1.0;
 			/*If epl thickness gets under colapse thickness, close the layer*/
 			if(epl_thickness[i]<colapse_thick){
-				eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
-				/*Avoid flipfloping between open and closed states*/
-				if(eplzigzag_counter[basalelement->nodes[i]->Lid()]<eplflip_lock |eplflip_lock==0){
-					vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
-					epl_thickness[i]=init_thick;
-				}
+				vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
+				new_active=0.0;
 			}
 		}
+
+		if(basalelement->nodes[i]->Sid()==42){
+			printf("old active is %e, and new active is %e, with counter= %i\n",old_active[i],new_active,eplzigzag_counter[basalelement->nodes[i]->Lid()]);
+		}
 		/*Increase of the efficient system is needed if the epl head reach the maximum value (sediment max value for now)*/
 		GetHydrologyDCInefficientHmax(&h_max,basalelement,basalelement->nodes[i]);
 		if(eplhead[i]>=h_max && active_element){
 			for(j=0;j<numnodes;j++){
 				/*Increase of the domain is on the downstream node in term of sediment head*/
 				if(sedhead[j] == sedheadmin){
-					if(old_active[j]==0) {
-						eplzigzag_counter[basalelement->nodes[j]->Lid()] ++;
-						if(eplzigzag_counter[basalelement->nodes[i]->Lid()]<eplflip_lock |eplflip_lock==0){
-							vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
-						}
-						else{
-							vec_mask->SetValue(basalelement->nodes[j]->Sid(),0.,INS_VAL);
-						}
-					}
+					vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
 				}
 			}
 		}
+		/*Now checking for nodes zigzaging between open and close*/
+		if (old_active[i] != new_active){
+			eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
+			/*If node changed too much of state, fix it to it's last known state*/
+			if(eplzigzag_counter[basalelement->nodes[i]->Lid()]>eplflip_lock & eplflip_lock!=0){
+				vec_mask->SetValue(basalelement->nodes[i]->Sid(),old_active[i],INS_VAL);
+				new_active=old_active[i];
+			}
+		}
+			/*If node is now closed bring its thickness back to initial*/
+		if (new_active==0.){
+			epl_thickness[i]=init_thick;
+		}
 	}
 	basalelement->AddInput(HydrologydcEplThicknessEnum,epl_thickness,basalelement->GetElementType());
 
Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h	(revision 19090)
+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h	(revision 19091)
@@ -39,7 +39,7 @@
 		IssmDouble SedimentStoring(Element* element);
 		IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
 		IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* sed_head_input, Input* epl_head_input, Input* thick_input, Input* base_input);
-		void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element);
+		void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, Vector<IssmDouble>* recurence, int* eplzigzag_counter, Element* element);
 		void HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element);
 		void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode);
 		void ComputeEPLThickness(FemModel* femmodel);
Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 19090)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 19091)
@@ -91,7 +91,7 @@
 		femmodel->UpdateConstraintsx();
 		femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum);
 		
-		/*Reset constraint on the ZigZag Lock, this thing doesn't work, it have to disapear*/
+		/*Reset constraint on the ZigZag Lock*/
 		ResetConstraintsx(femmodel);
 
 		/*{{{*//*Treating the sediment*/
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 19090)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 19091)
@@ -1912,6 +1912,7 @@
 void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/
 
 	Vector<IssmDouble>* mask							= NULL;
+	Vector<IssmDouble>* recurence  				= NULL;
 	IssmDouble*         serial_mask				= NULL;
 	Vector<IssmDouble>* active						= NULL;
 	IssmDouble*         serial_active			= NULL;
@@ -1921,11 +1922,11 @@
 	HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
 	/*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/
 	mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
+	recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
 	this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); 
-
 	for (int i=0;i<elements->Size();i++){
 		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-		effanalysis->HydrologyEPLGetMask(mask,eplzigzag_counter,element);
+		effanalysis->HydrologyEPLGetMask(mask,recurence,eplzigzag_counter,element);
 	}
 	this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum);
 	/*Assemble and serialize*/
@@ -1933,6 +1934,7 @@
 	serial_mask=mask->ToMPISerial();
 	xDelete<int>(eplzigzag_counter);
 	delete mask;
+	delete recurence;
 
 	/*Update Mask*/
 	InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
