Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 18658)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 18659)
@@ -1189,12 +1189,79 @@
 }/*}}}*/
 void EnthalpyAnalysis::UpdateBasalConstraints(FemModel* femmodel){/*{{{*/
+
 	/*Update basal dirichlet BCs for enthalpy: */
+	Vector<IssmDouble>* spc           = NULL;
+	IssmDouble*         serial_spc    = NULL;
+
+	spc=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes(EnthalpyAnalysisEnum));
+	/*First create a vector to figure out what elements should be constrained*/
 	for(int i=0;i<femmodel->elements->Size();i++){
 		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		UpdateBasalConstraints(element);
-	}
+		GetBasalConstraints(spc,element);
+	}
+
+	/*Assemble and serialize*/
+	spc->Assemble();
+	serial_spc=spc->ToMPISerial();
+	delete spc;
+
+	/*Then update basal constraints nodes accordingly*/
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		ApplyBasalConstraints(serial_spc,element);
+	}
+
 	femmodel->UpdateConstraintsx();
-}/*}}}*/
-void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/
+
+	/*Delete*/
+	xDelete<IssmDouble>(serial_spc);
+}/*}}}*/
+void EnthalpyAnalysis::ApplyBasalConstraints(IssmDouble* serial_spc,Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return;
+
+	/* Only update Constraints at the base of grounded ice*/
+	if(!(element->IsOnBase()) || element->IsFloating()) return;
+
+	/*Intermediary*/
+	bool        isdynamicbasalspc;
+	int         numindices;
+	int        *indices = NULL;
+	Node*       node = NULL;
+	IssmDouble	pressure;
+
+	/*Check wether dynamic basal boundary conditions are activated */
+	element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
+	if(!isdynamicbasalspc) return;
+
+	/*Get parameters and inputs: */
+	Input* pressure_input		 = element->GetInput(PressureEnum);							 _assert_(pressure_input);
+
+	/*Fetch indices of basal & surface nodes for this finite element*/
+	Penta *penta =  (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
+	penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
+
+	GaussPenta* gauss=new GaussPenta();
+	for(int i=0;i<numindices;i++){
+		gauss->GaussNode(element->GetElementType(),indices[i]);
+
+		pressure_input->GetInputValue(&pressure,gauss);
+
+		/*apply or release spc*/
+		node=element->GetNode(indices[i]);
+		if(serial_spc[node->Sid()]==1.){
+			pressure_input->GetInputValue(&pressure, gauss);
+			node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
+		}
+		else			
+			node->DofInFSet(0);
+	}
+
+	/*Free ressources:*/
+	xDelete<int>(indices);
+	delete gauss;
+}/*}}}*/
+void EnthalpyAnalysis::GetBasalConstraints(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
 
 	/* Check if ice in element */
@@ -1214,11 +1281,11 @@
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
 	if(dt==0.){
-		UpdateBasalConstraintsSteadystate(element);
+		GetBasalConstraintsSteadystate(vec_spc,element);
 	}
 	else{
-		UpdateBasalConstraintsTransient(element);
-	}
-}/*}}}*/
-void EnthalpyAnalysis::UpdateBasalConstraintsTransient(Element* element){/*{{{*/
+		GetBasalConstraintsTransient(vec_spc,element);
+	}
+}/*}}}*/
+void EnthalpyAnalysis::GetBasalConstraintsTransient(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
 
 	/* Check if ice in element */
@@ -1229,13 +1296,7 @@
 
 	/*Intermediary*/
-	bool        isdynamicbasalspc,setspc;
 	int         numindices, numindicesup, state;
 	int        *indices = NULL, *indicesup = NULL;
-	Node*       node = NULL;
 	IssmDouble	enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
-
-	/*Check wether dynamic basal boundary conditions are activated */
-	element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
-	if(!isdynamicbasalspc) return;
 
 	/*Get parameters and inputs: */
@@ -1266,22 +1327,24 @@
 		state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
 
-		setspc=false;
 		switch (state) {
 			case 0:
 				// cold, dry base: apply basal surface forcing
+				vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
 				break;
 			case 1:
 				// cold, wet base: keep at pressure melting point 
-				setspc=true;
+				vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
 				break;
 			case 2:
 				// temperate, thin refreezing base: release spc
+				vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
 				break;
 			case 3:
 				// temperate, thin melting base: set spc
-				setspc=true;
+				vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
 				break;
 			case 4:
 				// temperate, thick melting base: set grad H*n=0
+				vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
 				break;
 			default:
@@ -1289,12 +1352,4 @@
 		}
 
-		/*apply or release spc*/
-		node=element->GetNode(indices[i]);
-		if(setspc){
-			pressure_input->GetInputValue(&pressure, gauss);
-			node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
-		}
-		else			
-			node->DofInFSet(0);
 	}
 
@@ -1305,5 +1360,5 @@
 	delete gaussup;
 }/*}}}*/
-void EnthalpyAnalysis::UpdateBasalConstraintsSteadystate(Element* element){/*{{{*/
+void EnthalpyAnalysis::GetBasalConstraintsSteadystate(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
 
 	/* Check if ice in element */
@@ -1314,13 +1369,7 @@
 
 	/*Intermediary*/
-	bool        isdynamicbasalspc,setspc;
 	int         numindices, numindicesup, state;
 	int        *indices = NULL, *indicesup = NULL;
-	Node*       node = NULL;
 	IssmDouble	enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
-
-	/*Check wether dynamic basal boundary conditions are activated */
-	element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
-	if(!isdynamicbasalspc) return;
 
 	/*Get parameters and inputs: */
@@ -1349,36 +1398,28 @@
 
 		state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
-		setspc=false;
 		switch (state) {
 			case 0:
 				// cold, dry base: apply basal surface forcing
+				vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
 				break;
 			case 1:
 				// cold, wet base: keep at pressure melting point 
-				setspc=true;
+				vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
 				break;
 			case 2:
 				// temperate, thin refreezing base: release spc
+				vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
 				break;
 			case 3:
 				// temperate, thin melting base: set spc
-				setspc=true;
+				vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
 				break;
 			case 4:
 				// temperate, thick melting base: s
-				setspc=true;
+				vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
 				break;
 			default:
 				_printf0_("	unknown thermal basal state found!");
 		}
-
-		/*apply or release spc*/
-		node=element->GetNode(indices[i]);
-		if(setspc){
-			pressure_input->GetInputValue(&pressure, gauss);
-			node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
-		}
-		else			
-			node->DofInFSet(0);
 	}
 
Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 18658)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 18659)
@@ -48,7 +48,8 @@
 		static void DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element);
 		static void UpdateBasalConstraints(FemModel* femmodel);
-		static void UpdateBasalConstraints(Element* element);
-		static void UpdateBasalConstraintsTransient(Element* element);
-		static void UpdateBasalConstraintsSteadystate(Element* element);
+		static void ApplyBasalConstraints(IssmDouble* serial_spc,Element* element);
+		static void GetBasalConstraints(Vector<IssmDouble>* vec_spc,Element* element);
+		static void GetBasalConstraintsTransient(Vector<IssmDouble>* vec_spc,Element* element);
+		static void GetBasalConstraintsSteadystate(Vector<IssmDouble>* vec_spc,Element* element);
 		static int GetThermalBasalCondition(Element* element, IssmDouble enthalpy, IssmDouble enthalpy_up, IssmDouble pressure, IssmDouble pressure_up, IssmDouble watercolumn, IssmDouble meltingrate);
 		static IssmDouble GetWetIceConductivity(Element* element, IssmDouble enthalpy, IssmDouble pressure);