Index: issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 18903)
+++ issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 18975)
@@ -24,10 +24,11 @@
 	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){/*{{{*/
@@ -35,5 +36,4 @@
 	bool   isefficientlayer;
 	int    hydrology_model;
-
 
 	/*Now, do we really want DC?*/
@@ -61,7 +61,6 @@
 	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);
@@ -85,4 +84,5 @@
 	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
 }/*}}}*/
+
 void HydrologyDCEfficientAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
@@ -98,6 +98,6 @@
 
 	IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum);
-
-}/*}}}*/
+}/*}}}*/
+
 void HydrologyDCEfficientAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
 	/*Nothing for now*/
@@ -131,14 +131,14 @@
 }/*}}}*/
 
- 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*/
@@ -244,6 +244,6 @@
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return Ke;
-
-}/*}}}*/
+}/*}}}*/
+
 ElementVector* HydrologyDCEfficientAnalysis::CreatePVector(Element* element){/*{{{*/
 
@@ -348,10 +348,13 @@
 	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){/*{{{*/
 
@@ -375,5 +378,4 @@
 	/*Fetch number of nodes for this finite element*/
 	int numnodes = basalelement->GetNumberOfNodes();
-
 
 	/*Fetch dof list and allocate solution vector*/
@@ -401,4 +403,5 @@
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 } /*}}}*/
+
 void HydrologyDCEfficientAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
@@ -415,4 +418,5 @@
 	return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity));		 
 }/*}}}*/
+
 IssmDouble HydrologyDCEfficientAnalysis::SedimentStoring(Element* element){/*{{{*/
 	IssmDouble rho_freshwater           = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
@@ -535,5 +539,5 @@
 		_error_("not Implemented Yet");
 		}
-			
+		
 		int         numnodes      = element->GetNumberOfNodes();
 		IssmDouble* thickness     = xNew<IssmDouble>(numnodes);
@@ -602,6 +606,6 @@
 		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. 
@@ -631,4 +635,5 @@
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
+
 void  HydrologyDCEfficientAnalysis::HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask, int* eplzigzag_counter, Element* element){
 
@@ -655,5 +660,5 @@
 	/*Intermediaries*/
 
-	int         penalty_lock;
+	int         eplflip_lock;
 	int         numnodes      =basalelement->GetNumberOfNodes();
 	IssmDouble* epl_thickness =xNew<IssmDouble>(numnodes);
@@ -669,5 +674,5 @@
 	active_element_input->GetInputValue(&active_element);
 
-	basalelement->parameters->FindParam(&penalty_lock,HydrologydcEplflipLockEnum); 
+	basalelement->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); 
 
 	basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum);	
@@ -691,8 +696,9 @@
 		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;
@@ -705,5 +711,4 @@
 		/*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++){
@@ -729,4 +734,5 @@
 }
 /*}}}*/
+
 void HydrologyDCEfficientAnalysis::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec, Element* element){
 	/*Constants*/
@@ -750,5 +756,6 @@
 	IssmDouble  flag     = 0.;
 	IssmDouble* active   = xNew<IssmDouble>(numnodes);
-		
+
+	/*Pass the activity mask from elements to nodes*/
 	basalelement->GetInputListOnVertices(&active[0],HydrologydcMaskEplactiveNodeEnum);
 	bool active_element;
@@ -758,4 +765,5 @@
 	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++){
@@ -763,4 +771,5 @@
 		}
 	}
+	/*If the element is active all its nodes are active*/
 	else if(active_element){
 		for(int i=0;i<numnodes;i++){
