Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16787)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16788)
@@ -9519,6 +9519,7 @@
 void Penta::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){
 
-	if (!IsOnBed())return;
-
+	if (!IsOnBed()){
+		return;
+	}
 	Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
 	tria->HydrologyEPLGetActive(active_vec);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16787)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16788)
@@ -6989,10 +6989,9 @@
 		
 		for(int i=0;i<numdof;i++){
-			/*Keeping thickness to 1 if EPL is not active*/
-			if(activeEpl[i]==0.0){
-				thickness[i]=init_thick;
-			}
-			else{
-
+			/*Keeping thickness to initial value if EPL is not active*/
+			/* if(activeEpl[i]==0.0){ */
+			/* 	thickness[i]=init_thick; */
+			/* } */
+			/* else{ */
 				/*Compute first the effective pressure in the EPL*/
 				EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
@@ -7003,5 +7002,5 @@
 				/*And proceed to the real thing*/
 				thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
-			}
+				//}
 		}
 		this->inputs->AddInput(new TriaInput(HydrologydcEplThicknessEnum,thickness,P1Enum));
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 16787)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 16788)
@@ -1423,8 +1423,47 @@
 void FemModel::UpdateConstraintsL2ProjectionEPLx(void){ /*{{{*/
 
-	for(int i=0;i<elements->Size();i++){
+	Vector<IssmDouble>* active        = NULL;
+	IssmDouble*         serial_active = NULL;
+
+	/*update node activity. If one element is connected to mask=1, all nodes are active*/
+	active=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
+	for (int i=0;i<elements->Size();i++){
 		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-		element->UpdateConstraintsL2ProjectionEPL();
-	}
+		element->HydrologyEPLGetActive(active);
+	}
+
+	/*Assemble and serialize*/
+	active->Assemble();
+	serial_active=active->ToMPISerial();
+	delete active;
+
+	/*Update node activation accordingly*/
+	int counter =0;
+	for (int i=0;i<nodes->Size();i++){
+		Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
+		if(node->InAnalysis(L2ProjectionEPLAnalysisEnum)){
+			if(serial_active[node->Sid()]==1.){
+				node->Activate();
+				counter++;
+			}
+			else{
+				node->Deactivate();
+			}
+		}
+	}
+	xDelete<IssmDouble>(serial_active);
+	int sum_counter;
+	ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	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");
+
+	/*Update dof indexings*/
+	this->UpdateConstraintsx();
+
+	/* for(int i=0;i<elements->Size();i++){ */
+	/* 	Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i)); */
+	/* 	element->UpdateConstraintsL2ProjectionEPL(); */
+	/* } */
 
 }
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 16787)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 16788)
@@ -648,4 +648,10 @@
 		return;
 	}
+	if(!element->IsOnBed()){
+		unstable=0;
+		active=0;
+		*punstable=unstable;
+		return;
+	}
 
 	/*Get sediment water head h*/
@@ -654,8 +660,10 @@
 	parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum);
 
-	if (h>h_max)
+	if (h>h_max){
 	 new_active=1;
-	else
+	}
+	else{
 	 new_active=0;
+	}
 
 	if(this->active==new_active){