Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 18718)
+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 18719)
@@ -30,6 +30,7 @@
 	bool   isefficientlayer;
 	int    hydrology_model;
 
+
 	/*Now, do we really want DC?*/
 	iomodel->Constant(&hydrology_model,HydrologyModelEnum);
 	if(hydrology_model!=HydrologydcEnum) return;
@@ -60,8 +61,8 @@
 		iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
 		iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
 	}
+}/*}}}*/
 
-}/*}}}*/
 void HydrologyDCEfficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
 	/*Now, do we really want DC?*/
@@ -97,19 +98,44 @@
 	/*Nothing for now*/
 }/*}}}*/
 
+void HydrologyDCEfficientAnalysis::InitZigZagCounter(FemModel* femmodel){/*{{{*/
+	int*   eplzigzag_counter =NULL;
+	
+	eplzigzag_counter=xNewZeroInit<int>(femmodel->nodes->Size());
+	femmodel->parameters->AddObject(new IntVecParam(EplZigZagCounterEnum,eplzigzag_counter,femmodel->nodes->Size()));
+	xDelete<int>(eplzigzag_counter);
+}/*}}}*/
+
+void HydrologyDCEfficientAnalysis::ResetCounter(FemModel* femmodel){/*{{{*/
+
+	int*     eplzigzag_counter=NULL;
+	Element* element=NULL;
+
+	femmodel->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
+	for(int i=0;i<femmodel->nodes->Size();i++){
+
+		eplzigzag_counter[i]=0;
+	}
+	femmodel->parameters->SetParam(eplzigzag_counter,femmodel->nodes->Size(),EplZigZagCounterEnum);
+	xDelete<int>(eplzigzag_counter);
+}/*}}}*/
+
 /*Finite Element Analysis*/
 void HydrologyDCEfficientAnalysis::Core(FemModel* femmodel){/*{{{*/
 	_error_("not implemented");
 }/*}}}*/
-ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/
+
+ ElementVector* HydrologyDCEfficientAnalysis::CreateDVector(Element* element){/*{{{*/
 	/*Default, return NULL*/
 	return NULL;
 }/*}}}*/
-ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
+
+ ElementMatrix* HydrologyDCEfficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
 _error_("Not implemented");
 }/*}}}*/
-ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
+ ElementMatrix* HydrologyDCEfficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
+
 	/*Intermediaries*/
 	bool     active_element;
 	int      domaintype;
@@ -611,14 +637,14 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
-void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask,Element* element){
+void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element){
 
 	bool        active_element;
 	int         i,j;
 	int         domaintype;
 	IssmDouble  h_max;
 	IssmDouble  sedheadmin;
-	Element*   basalelement=NULL;
+	Element*    basalelement=NULL;
 
 	/*Get basal element*/
 	element->FindParam(&domaintype,DomainTypeEnum);
@@ -635,19 +661,22 @@
 
 	/*Intermediaries*/
 
+	int         penalty_lock;
 	int         numnodes      =basalelement->GetNumberOfNodes();
 	IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
 	IssmDouble* old_active    =xNew<IssmDouble>(numnodes);
 	IssmDouble* sedhead       =xNew<IssmDouble>(numnodes);
 	IssmDouble* eplhead       =xNew<IssmDouble>(numnodes);
 	IssmDouble* residual      =xNew<IssmDouble>(numnodes);
-
+	
 	IssmDouble init_thick    =basalelement->GetMaterialParameter(HydrologydcEplInitialThicknessEnum);
 	IssmDouble colapse_thick =basalelement->GetMaterialParameter(HydrologydcEplColapseThicknessEnum);
 
 	Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
 	active_element_input->GetInputValue(&active_element);
 
+	basalelement->parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum); 
+
 	basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum);	
 	basalelement-> GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum);	
 	basalelement-> GetInputListOnVertices(&sedhead[0],SedimentHeadEnum);
@@ -661,6 +690,7 @@
 	for(i=0;i<numnodes;i++){
 		/*Activate EPL if residual is >0 */
 		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);
 		}
 
@@ -669,27 +699,34 @@
 			vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
 			/*If epl thickness gets under 10-3 initial thickness, close the layer*/
 			if(epl_thickness[i]<colapse_thick){
-				vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
-				epl_thickness[i]=init_thick;
+				eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
+				if(eplzigzag_counter[basalelement->nodes[i]->Sid()]<penalty_lock |penalty_lock==0){
+					vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL);
+					epl_thickness[i]=init_thick;
+				}
+				else{
+					vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
+				}
 			}
 		}
 		/*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++){
 				if(old_active[j]>0.){
 					vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
 				}
 				/*Increase of the domain is on the downstream node in term of sediment head*/
 				if(sedhead[j] == sedheadmin){
+					if(old_active[i]==0) eplzigzag_counter[basalelement->nodes[j]->Lid()] ++;
 					vec_mask->SetValue(basalelement->nodes[j]->Sid(),1.,INS_VAL);
 				}
 			}
 		}
-		
 	}
+	basalelement->AddInput(HydrologydcEplThicknessEnum,epl_thickness,basalelement->GetElementType());
+
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	xDelete<IssmDouble>(epl_thickness);
 	xDelete<IssmDouble>(old_active);
@@ -722,7 +759,13 @@
 	IssmDouble* active   = xNew<IssmDouble>(numnodes);
 		
 	basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
+
+	bool active_element;
 	
+	Input* 	active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);		
+	active_element_input->GetInputValue(&active_element);
+
+	
 	for(int i=0;i<numnodes;i++) flag+=active[i];
 
 	if(flag>0.){
@@ -730,6 +773,12 @@
 			active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
 		}
 	}
+	else if(active_element){
+		/*Checking Stuff*/
+		for(int i=0;i<numnodes;i++){
+			active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
+		}		
+	}
 	else{
 		/*Do not do anything: at least one node was active for this element but this element is not solved for*/
 	}
Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h	(revision 18718)
+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h	(revision 18719)
@@ -19,7 +19,9 @@
 		void CreateNodes(Nodes* nodes,IoModel* iomodel);
 		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
 		void CreateLoads(Loads* loads, IoModel* iomodel);
-
+		void InitZigZagCounter(FemModel* femmodel);
+		void ResetCounter(FemModel* femmodel);
+			
 		/*Finite element Analysis*/
 		void           Core(FemModel* femmodel);
 		ElementVector* CreateDVector(Element* element);
@@ -37,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,Element* element);
+		void HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, 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/analyses/HydrologyDCInefficientAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h	(revision 18718)
+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h	(revision 18719)
@@ -40,6 +40,5 @@
 		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 ElementizeEplMask(FemModel* femmodel);
-		void ResetCounter(FemModel* femmodel);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 18718)
+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 18719)
@@ -16,7 +16,7 @@
 	int         transfer_flag;
 	int         penalty_lock;
 	int         hydro_maxiter;
-	int*        elementactive_counter =NULL;
+	/* int*        elementactive_counter =NULL; */
 	bool        isefficientlayer;
 	IssmDouble  sedimentlimit;
 	IssmDouble  penalty_factor;
@@ -55,12 +55,6 @@
 	parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
 	parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock));
 	parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter));
-
-	if(isefficientlayer==1){
-		elementactive_counter=xNewZeroInit<int>(iomodel->numberofelements);
-		parameters->AddObject(new IntVecParam(ElementActiveCounterEnum,elementactive_counter,iomodel->numberofelements));
-		xDelete<int>(elementactive_counter);
-	}
 	
 }/*}}}*/
 void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
@@ -103,18 +97,18 @@
 	if(isefficientlayer){
 		iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum);
 
-		for(int i=0;i<elements->Size();i++){
-			Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-			Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input);
+		/* for(int i=0;i<elements->Size();i++){ */
+		/* 	Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i)); */
+		/* 	Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input); */
 			
-			if(node_mask_input->Max()>0.) {
-				element_active = true;
-			}
-			else{
-				element_active = false;
-			}
-			element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active));
-		}
+		/* 	if(node_mask_input->Max()>0.) { */
+		/* 		element_active = true; */
+		/* 	} */
+		/* 	else{ */
+		/* 		element_active = false; */
+		/* 	} */
+		/* 	element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active)); */
+		/* } */
 	}
 }/*}}}*/
 void HydrologyDCInefficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
@@ -683,50 +677,18 @@
 void HydrologyDCInefficientAnalysis::ElementizeEplMask(FemModel* femmodel){/*{{{*/
 
 	bool     element_active;
-	bool     old_active;
-	int      penalty_lock;
-	int*     elementactive_counter=NULL;
 	Element* element=NULL;
 
-	femmodel->parameters->FindParam(&elementactive_counter,NULL,ElementActiveCounterEnum);
 	for(int i=0;i<femmodel->elements->Size();i++){
 		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		element->parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum);
+			Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input);
 		
-		Input* node_mask_input = element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(node_mask_input);
-		Input* old_active_input = element->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(old_active_input);
-		old_active_input->GetInputValue(&old_active);
-		
 		if(node_mask_input->Max()>0.){
 			element_active = true;
 		}
 		else{
 			element_active = false;
 		}
-		if(old_active!=element_active){
-			elementactive_counter[i]=elementactive_counter[i]+1;
-		}
-		if(elementactive_counter[i]>penalty_lock){
-			element_active = true;
-		}
 		element->AddInput(new BoolInput(HydrologydcMaskEplactiveEltEnum,element_active));
 	}
-	femmodel->parameters->SetParam(elementactive_counter,femmodel->elements->Size(),ElementActiveCounterEnum);
-	xDelete<int>(elementactive_counter);
 }/*}}}*/
-
-void HydrologyDCInefficientAnalysis::ResetCounter(FemModel* femmodel){/*{{{*/
-
-	int*     elementactive_counter=NULL;
-	Element* element=NULL;
-	int N;
-
-	femmodel->parameters->FindParam(&elementactive_counter,&N,ElementActiveCounterEnum);
-	for(int i=0;i<femmodel->elements->Size();i++){
-		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		
-		elementactive_counter[i]=0;
-	}		
-	femmodel->parameters->SetParam(elementactive_counter,femmodel->elements->Size(),ElementActiveCounterEnum);
-	xDelete<int>(elementactive_counter);
-}/*}}}*/
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18718)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18719)
@@ -108,7 +108,7 @@
 	EplHeadOldEnum,
 	EplHeadSlopeXEnum,
 	EplHeadSlopeYEnum,
-	ElementActiveCounterEnum,
+	EplZigZagCounterEnum,
 	HydrologydcMaxIterEnum,
 	HydrologydcRelTolEnum,
 	HydrologydcSpcsedimentHeadEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18718)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18719)
@@ -116,7 +116,7 @@
 		case EplHeadOldEnum : return "EplHeadOld";
 		case EplHeadSlopeXEnum : return "EplHeadSlopeX";
 		case EplHeadSlopeYEnum : return "EplHeadSlopeY";
-		case ElementActiveCounterEnum : return "ElementActiveCounter";
+		case EplZigZagCounterEnum : return "EplZigZagCounter";
 		case HydrologydcMaxIterEnum : return "HydrologydcMaxIter";
 		case HydrologydcRelTolEnum : return "HydrologydcRelTol";
 		case HydrologydcSpcsedimentHeadEnum : return "HydrologydcSpcsedimentHead";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18718)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18719)
@@ -116,7 +116,7 @@
 	      else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
 	      else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
 	      else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
-	      else if (strcmp(name,"ElementActiveCounter")==0) return ElementActiveCounterEnum;
+	      else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum;
 	      else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum;
 	      else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum;
 	      else if (strcmp(name,"HydrologydcSpcsedimentHead")==0) return HydrologydcSpcsedimentHeadEnum;
Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 18718)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 18719)
@@ -68,6 +68,7 @@
 		GetSolutionFromInputsx(&ug_epl,femmodel);
 		/*Initialize the element mask*/
 		inefanalysis->ElementizeEplMask(femmodel);
+		effanalysis->InitZigZagCounter(femmodel);
 	}
 	/*The real computation starts here, outermost loop is on the two layer system*/
 	for(;;){
@@ -221,7 +222,7 @@
 					if(VerboseSolution()) _printf0_("eplconverged...\n");
 					InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum);
 					InputUpdateFromSolutionx(femmodel,ug_epl);
-					inefanalysis->ResetCounter(femmodel);
+					effanalysis->ResetCounter(femmodel);
 					break;
 				}
 			}
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18718)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18719)
@@ -1808,23 +1808,27 @@
 #endif
 void FemModel::HydrologyEPLupdateDomainx(void){ /*{{{*/
 
-	Vector<IssmDouble>* mask          = NULL;
-	IssmDouble*         serial_mask   = NULL;
-	Vector<IssmDouble>* active        = NULL;
-	IssmDouble*         serial_active = NULL;
-	
+	Vector<IssmDouble>* mask							= NULL;
+	IssmDouble*         serial_mask				= NULL;
+	Vector<IssmDouble>* active						= NULL;
+	IssmDouble*         serial_active			= NULL;
+	int*                eplzigzag_counter =	NULL;
+		
 	HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
 	/*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/
-	mask=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
+	mask=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,element);
+		effanalysis->HydrologyEPLGetMask(mask,eplzigzag_counter,element);
 	}
-
+	printf("POINTER = %p\n",eplzigzag_counter);
+	this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum);
 	/*Assemble and serialize*/
 	mask->Assemble();
 	serial_mask=mask->ToMPISerial();
+	xDelete<int>(eplzigzag_counter);
 	delete mask;
 
 	/*Update Mask*/
