Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18859)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18860)
@@ -1546,37 +1546,47 @@
 
 	if(!IsOnBase()) return;
 
-	/*Intermediaries*/
-	IssmDouble* xyz_list = NULL;
-	IssmDouble* xyz_list_base = NULL;
-	IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
-	IssmDouble  bed_normal[2];
-	IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
-	IssmDouble  surface=0,value=0;
-	bool grounded;
+	int approximation;
+	inputs->GetInputValue(&approximation,ApproximationEnum);
 
-	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list);
-	GetVerticesCoordinatesBase(&xyz_list_base);
+	if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
+		for(int i=0;i<NUMVERTICES;i++){
+			vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
+			vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
+		}
+	}
+	else{
+		/*Intermediaries*/
+		IssmDouble* xyz_list = NULL;
+		IssmDouble* xyz_list_base = NULL;
+		IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
+		IssmDouble  bed_normal[2];
+		IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+		IssmDouble  surface=0,value=0;
+		bool grounded;
 
-	/*Retrieve all inputs we will be needing: */
-	Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
-	Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
-	Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
-	Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
-	Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
+		/* Get node coordinates and dof list: */
+		GetVerticesCoordinates(&xyz_list);
+		GetVerticesCoordinatesBase(&xyz_list_base);
 
-	/*Create gauss point in the middle of the basal edge*/
-	Gauss* gauss=NewGaussBase(1);
-	gauss->GaussPoint(0);
+		/*Retrieve all inputs we will be needing: */
+		Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
+		Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
+		Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
+		Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
+		Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
 
-	if(!IsFloating()){ 
-		/*Check for basal force only if grounded and touching GL*/
-		if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
+		/*Create gauss point in the middle of the basal edge*/
+		Gauss* gauss=NewGaussBase(1);
+		gauss->GaussPoint(0);
+
+		if(!IsFloating()){ 
+			/*Check for basal force only if grounded and touching GL*/
+			//		if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
 			this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
 			this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
 			pressure_input->GetInputValue(&pressure, gauss);
-			base_input->GetInputValue(&base, gauss); _assert_(base<0.);
+			base_input->GetInputValue(&base, gauss); 
 
 			/*Compute Stress*/
 			IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
@@ -1600,25 +1610,22 @@
 			else                         grounded=false;
 		}
 		else{
-			grounded=true;
+			/*Check for basal elevation if floating*/
+			base_input->GetInputValue(&base, gauss);
+			bed_input->GetInputValue(&bed, gauss);
+			if(base<bed) grounded=true;
+			else         grounded=false;
 		}
+		for(int i=0;i<NUMVERTICES;i++){
+			if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
+			else         vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
+		}
+
+		/*clean up*/
+		delete gauss;
+		xDelete<IssmDouble>(xyz_list);
+		xDelete<IssmDouble>(xyz_list_base);
 	}
-	else{
-		/*Check for basal elevation if floating*/
-		base_input->GetInputValue(&base, gauss);
-		bed_input->GetInputValue(&bed, gauss);
-		if (base<bed) grounded=true;
-		else          grounded=false;
-	}
-	for(int i=0;i<NUMVERTICES;i++){
-		if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
-		else         vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
-	}
-
-	/*clean up*/
-	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(xyz_list_base);
 }
 /*}}}*/
 bool       Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18859)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18860)
@@ -687,32 +687,40 @@
 
 	if(!IsOnBase()) return;
 
-	/*Intermediaries*/
-	IssmDouble* xyz_list = NULL;
-	IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
-	IssmDouble  bed_normal[3];
-	IssmDouble  epsilon[6]; /* epsilon=[exx eyy ezz exy exz eyz];*/
-	IssmDouble  surface=0,value=0;
-	bool grounded;
+	int approximation;
+	inputs->GetInputValue(&approximation,ApproximationEnum);
+	if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
+		for(int i=0;i<NUMVERTICES;i++){
+			vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
+			vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
+		}
+	}
+	else {
+		/*Intermediaries*/
+		IssmDouble* xyz_list = NULL;
+		IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
+		IssmDouble  bed_normal[3];
+		IssmDouble  epsilon[6]; /* epsilon=[exx eyy ezz exy exz eyz];*/
+		IssmDouble  surface=0,value=0;
+		bool grounded;
 
-	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list);
+		/* Get node coordinates and dof list: */
+		GetVerticesCoordinates(&xyz_list);
 
-	/*Retrieve all inputs we will be needing: */
-	Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
-	Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
-	Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
-	Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
-	Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
-	Input* vz_input       = inputs->GetInput(VzEnum);       _assert_(vz_input);
+		/*Retrieve all inputs we will be needing: */
+		Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
+		Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
+		Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
+		Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
+		Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
+		Input* vz_input       = inputs->GetInput(VzEnum);       _assert_(vz_input);
 
-	/*Create gauss point in the middle of the basal edge*/
-	Gauss* gauss=NewGaussBase(1);
-	gauss->GaussPoint(0);
+		/*Create gauss point in the middle of the basal edge*/
+		Gauss* gauss=NewGaussBase(1);
+		gauss->GaussPoint(0);
 
-	if(!IsFloating()){ 
-		/*Check for basal force only if grounded and touching GL*/
-		if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
+		if(!IsFloating()){ 
+			/*Check for basal force only if grounded and touching GL*/
 			this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
 			this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
 			pressure_input->GetInputValue(&pressure, gauss);
@@ -744,24 +752,21 @@
 			else                         grounded=false;
 		}
 		else{
-			grounded=true;
+			/*Check for basal elevation if floating*/
+			base_input->GetInputValue(&base, gauss);
+			bed_input->GetInputValue(&bed, gauss);
+			if(base<bed) grounded=true;
+			else          grounded=false;
 		}
+		for(int i=0;i<NUMVERTICES;i++){
+			if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
+			else         vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
+		}
+
+		/*clean up*/
+		delete gauss;
+		xDelete<IssmDouble>(xyz_list);
 	}
-	else{
-		/*Check for basal elevation if floating*/
-		base_input->GetInputValue(&base, gauss);
-		bed_input->GetInputValue(&bed, gauss);
-		if (base<bed) grounded=true;
-		else          grounded=false;
-	}
-	for(int i=0;i<NUMVERTICES;i++){
-		if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
-		else         vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
-	}
-
-	/*clean up*/
-	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
 int        Penta::ObjectEnum(void){/*{{{*/
