Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17143)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 17144)
@@ -56,4 +56,5 @@
 
 	if(VerboseSolution()) _printf0_("extrapolation: call computational core:\n");
+	UpdateConstraints(femmodel->elements);
 	solutionsequence_linear(femmodel);
 
@@ -230,3 +231,42 @@
 
 }/*}}}*/
-
+void ExtrapolationAnalysis::UpdateConstraints(Elements* elements){/*{{{*/
+
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+		this->SetConstraintsOnIce(element);
+	}
+
+}/*}}}*/
+void ExtrapolationAnalysis::SetConstraintsOnIce(Element* element){/*{{{*/
+
+	int numnodes=element->GetNumberOfNodes();	
+
+	/* Intermediaries */
+	int extvar_enum;
+	IssmDouble phi,value;
+	Node* node = NULL;
+
+	/* Get further parameters */
+	element->FindParam(&extvar_enum, ExtrapolationVariableEnum);
+	
+	Input* levelset_input=element->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input);
+	Input* extvar_input=element->GetInput(extvar_enum); _assert_(extvar_input);
+	
+	GaussPenta* gauss=new GaussPenta();
+	for(int in=0;in<numnodes;in++){
+		gauss->GaussNode(element->GetElementType(),in);
+		node=element->GetNode(in);
+		levelset_input->GetInputValue(&phi,gauss);
+		extvar_input->GetInputValue(&value,gauss);
+		if(phi<=0.){
+			/* if ice, set dirichlet BC */
+			node->ApplyConstraint(1,value);
+		}
+		else {
+			/* no ice, set no spc */
+			node->DofInFSet(0); 
+		}
+	}
+	delete gauss;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h	(revision 17143)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h	(revision 17144)
@@ -30,5 +30,6 @@
 	void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 	void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-
+	void UpdateConstraints(Elements* elements);
+	void SetConstraintsOnIce(Element* element);
 };
 #endif