Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 18974)
+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 18975)
@@ -23,19 +23,19 @@
 	/*Do we want an efficient layer*/
 	iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
 
+	/*If not return*/
 	if(!isefficientlayer) return;
+
+	/*If yes, initialize a flip flop counter*/
 	iomodel->FetchData(&eplflip_lock,HydrologydcEplflipLockEnum);
- 
 	parameters->AddObject(new IntParam(HydrologydcEplflipLockEnum,eplflip_lock));
 	
-	/*Nothing for now*/
 }/*}}}*/
 void HydrologyDCEfficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
 	bool   isefficientlayer;
 	int    hydrology_model;
 
-
 	/*Now, do we really want DC?*/
 	iomodel->Constant(&hydrology_model,HydrologyModelEnum);
 	if(hydrology_model!=HydrologydcEnum) return;
@@ -60,9 +60,8 @@
 	iomodel->FetchDataToInput(elements,SedimentHeadEnum);
 	iomodel->FetchDataToInput(elements,HydrologydcEplInitialThicknessEnum);
 	iomodel->FetchDataToInput(elements,HydrologydcEplMaxThicknessEnum);
-	iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum);
 	iomodel->FetchDataToInput(elements,HydrologydcEplThicknessEnum);
-	if(iomodel->domaintype!=Domain2DhorizontalEnum){
+		if(iomodel->domaintype!=Domain2DhorizontalEnum){
 		iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
 		iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
 	}
@@ -84,6 +83,7 @@
 	::CreateNodes(nodes,iomodel,HydrologyDCEfficientAnalysisEnum,P1Enum);
 	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
 	/*Do we really want DC?*/
@@ -97,8 +97,8 @@
 	if(!isefficientlayer) return;
 
 	IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum);
+}/*}}}*/
 
-}/*}}}*/
 void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
 	/*Nothing for now*/
 }/*}}}*/
@@ -130,16 +130,16 @@
 	_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;
@@ -243,8 +243,8 @@
 	delete gauss;
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return Ke;
+}/*}}}*/
 
-}/*}}}*/
 ElementVector* HydrologyDCEfficientAnalysis::CreatePVector(Element* element){/*{{{*/
 
 	/*Intermediaries*/
@@ -347,12 +347,15 @@
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return pe;
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	element->GetSolutionFromInputsOneDof(solution,EplHeadEnum);
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
 	bool active_element;
@@ -375,7 +378,6 @@
 	/*Fetch number of nodes for this finite element*/
 	int numnodes = basalelement->GetNumberOfNodes();
 
-
 	/*Fetch dof list and allocate solution vector*/
 	IssmDouble* sedhead     = xNew<IssmDouble>(numnodes);
 	IssmDouble* eplHeads    = xNew<IssmDouble>(numnodes);
@@ -400,6 +402,7 @@
 	xDelete<int>(doflist);
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 } /*}}}*/
+
 void HydrologyDCEfficientAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
 	return;
@@ -414,6 +417,7 @@
 	IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
 	return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity));		 
 }/*}}}*/
+
 IssmDouble HydrologyDCEfficientAnalysis::SedimentStoring(Element* element){/*{{{*/
 	IssmDouble rho_freshwater           = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
 	IssmDouble g                        = element->GetMaterialParameter(ConstantsGEnum);
@@ -534,7 +538,7 @@
 		default:
 		_error_("not Implemented Yet");
 		}
-			
+		
 		int         numnodes      = element->GetNumberOfNodes();
 		IssmDouble* thickness     = xNew<IssmDouble>(numnodes);
 		IssmDouble* eplhead       = xNew<IssmDouble>(numnodes);
@@ -601,8 +605,8 @@
 		xDelete<IssmDouble>(ice_thickness);
 		xDelete<IssmDouble>(bed);
 	}
-}
-/*}}}*/
+}/*}}}*/
+
 void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
@@ -630,6 +634,7 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
+
 void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element){
 
 	bool        active_element;
@@ -654,7 +659,7 @@
 
 	/*Intermediaries*/
 
-	int         penalty_lock;
+	int         eplflip_lock;
 	int         numnodes      =basalelement->GetNumberOfNodes();
 	IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
 	IssmDouble* old_active    =xNew<IssmDouble>(numnodes);
@@ -668,7 +673,7 @@
 	Input* active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);
 	active_element_input->GetInputValue(&active_element);
 
-	basalelement->parameters->FindParam(&penalty_lock,HydrologydcEplflipLockEnum); 
+	basalelement->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); 
 
 	basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum);	
 	basalelement-> GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum);	
@@ -690,10 +695,11 @@
 		/*If mask was already one, keep one*/
 		else if(old_active[i]>0.){
 			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 gets under colapse thickness, close the layer*/
 			if(epl_thickness[i]<colapse_thick){
 				eplzigzag_counter[basalelement->nodes[i]->Lid()] ++;
-				if(eplzigzag_counter[basalelement->nodes[i]->Lid()]<penalty_lock |penalty_lock==0){
+				/*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;
 				}
@@ -704,7 +710,6 @@
 		}
 		/*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.){
@@ -728,6 +733,7 @@
 	xDelete<IssmDouble>(residual);
 }
 /*}}}*/
+
 void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){
 	/*Constants*/
 	int      domaintype;
@@ -749,7 +755,8 @@
 	const int   numnodes = basalelement->GetNumberOfNodes();
 	IssmDouble  flag     = 0.;
 	IssmDouble* active   = xNew<IssmDouble>(numnodes);
-		
+
+	/*Pass the activity mask from elements to nodes*/
 	basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
 	bool active_element;
 	Input* 	active_element_input=basalelement->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input);		
@@ -757,11 +764,13 @@
 	
 	for(int i=0;i<numnodes;i++) flag+=active[i];
 
+	/*If any node is active all the node in the element are active*/
 	if(flag>0.){
 		for(int i=0;i<numnodes;i++){
 			active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
 		}
 	}
+	/*If the element is active all its nodes are active*/
 	else if(active_element){
 		for(int i=0;i<numnodes;i++){
 			active_vec->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL);
Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 18974)
+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 18975)
@@ -9,6 +9,7 @@
 int  HydrologyDCInefficientAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
 	return 1;
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 
 	int         hydrology_model;
@@ -16,7 +17,6 @@
 	int         transfer_flag;
 	int         penalty_lock;
 	int         hydro_maxiter;
-	/* int*        elementactive_counter =NULL; */
 	bool        isefficientlayer;
 	IssmDouble  sedimentlimit;
 	IssmDouble  penalty_factor;
@@ -55,8 +55,8 @@
 	parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
 	parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock));
 	parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter));
-	
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
 	bool   isefficientlayer;
@@ -96,21 +96,9 @@
 
 	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); */
-			
-		/* 	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){/*{{{*/
 
 	/*Fetch parameters: */
@@ -124,6 +112,7 @@
 	::CreateNodes(nodes,iomodel,HydrologyDCInefficientAnalysisEnum,P1Enum);
 	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
 	/*retrieve some parameters: */
@@ -133,6 +122,7 @@
 
 	IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum,P1Enum);
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
 
 	/*Fetch parameters: */
@@ -164,13 +154,16 @@
 void           HydrologyDCInefficientAnalysis::Core(FemModel* femmodel){/*{{{*/
 	_error_("not implemented");
 }/*}}}*/
+
 ElementVector* HydrologyDCInefficientAnalysis::CreateDVector(Element* element){/*{{{*/
 	/*Default, return NULL*/
 	return NULL;
 }/*}}}*/
+
 ElementMatrix* HydrologyDCInefficientAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
 _error_("Not implemented");
 }/*}}}*/
+
 ElementMatrix* HydrologyDCInefficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
 	/*Intermediaries*/
@@ -274,6 +267,7 @@
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return Ke;
 }/*}}}*/
+
 ElementVector* HydrologyDCInefficientAnalysis::CreatePVector(Element* element){/*{{{*/
 
 	/*Intermediaries*/
@@ -311,7 +305,6 @@
 	ElementVector* pe    = basalelement->NewElementVector();
 	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
-
 	/*Retrieve all inputs and parameters*/
 	basalelement->GetVerticesCoordinates(&xyz_list);
 	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
@@ -341,9 +334,7 @@
 
 		/*Loading term*/
 		water_input->GetInputValue(&water_load,gauss);
-
 		scalar = Jdet*gauss->weight*(water_load);
-
 		if(dt!=0.) scalar = scalar*dt;
 		for(int i=0;i<numnodes;i++){
 			pe->values[i]+=scalar*basis[i];
@@ -377,6 +368,7 @@
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return pe;
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
@@ -404,12 +396,15 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	element->GetSolutionFromInputsOneDof(solution,SedimentHeadEnum);
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
 	int        domaintype;
@@ -488,6 +483,7 @@
 	xDelete<int>(doflist);
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
+
 void HydrologyDCInefficientAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
 	return;
@@ -550,6 +546,7 @@
 	}
 	return h_max;
 }/*}}}*/
+
 void  HydrologyDCInefficientAnalysis::GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode){/*{{{*/
 	
 	int        hmax_flag;
Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 18974)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 18975)
@@ -126,14 +126,6 @@
 			}
 		
 			/*}}}*//*End of the sediment penalization loop*/
-			/*Update EPL mask*/
-			
-			/*++++++This is probably useless+++++++*/
-			/* if(isefficientlayer){ */
-			/* 	inefanalysis->ElementizeEplMask(femmodel); */
-			/* } */
-			/*+++++++++++++*/
-
 			sedconverged=false;
 			
 			/*Checking convegence on the value of the sediment head*/
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18974)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18975)
@@ -1978,10 +1978,5 @@
 	ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
 	counter=sum_counter;
 	if(VerboseSolution()) _printf0_("   Number of active nodes L2 Projection: "<< counter <<"\n");
-
-	/*+++++++this is done in the solution sequence+++++++++++++++*/
-	/*Update dof indexings*/
-	//	this->UpdateConstraintsx();
-	/*++++++++++++++++++++++++++*/
 }
 /*}}}*/
