Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17632)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17633)
@@ -670,5 +670,5 @@
 		else shelf=false;
 	}
-	else if(migration_style==NoneEnum || migration_style==AgressiveMigrationEnum || migration_style==SoftMigrationEnum){ //Floating if all nodes are floating
+	else if(migration_style==NoneEnum || migration_style==AgressiveMigrationEnum || migration_style==SoftMigrationEnum || migration_style==ContactEnum){ //Floating if all nodes are floating
 		if(this->inputs->Min(MaskGroundediceLevelsetEnum) > 0.) shelf=false;
 		else shelf=true;
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17632)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17633)
@@ -1323,18 +1323,60 @@
 	/*Intermediaries*/
 	IssmDouble* xyz_list = NULL;
-
-	/*Figure out wether element is grounded or floating*/
-	bool grounded = true;
+	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;
+
+	IssmDouble* basis = xNew<IssmDouble>(NUMVERTICES);
+
 
 	/* Get node coordinates and dof list: */
-	::GetVerticesCoordinates(xyz_list,vertices,NUMVERTICES);
+	GetVerticesCoordinates(&xyz_list);
 
 	/*Retrieve all inputs we will be needing: */
-	Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
-	Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
-	Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
-
-	_error_("STOP");
-
+	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);
+
+	Gauss* gauss=NewGaussBase(1);
+	gauss->GaussPoint(0);
+
+	if(!IsFloating()){ 
+
+		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.);
+
+		/*Compute Stress*/
+		IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
+		IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
+		IssmDouble sigma_xy=2.*viscosity*epsilon[2];
+
+		/*Get normal vector to the bed */
+		NormalBase(&bed_normal[0],xyz_list);
+
+		/*basalforce*/
+		sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
+
+		/*Compute water pressure*/
+		IssmDouble rho_ice   = matpar->GetRhoIce();
+		IssmDouble rho_water = matpar->GetRhoWater();
+		IssmDouble gravity   = matpar->GetG();
+		water_pressure=gravity*rho_water*base;
+
+		/*Compare basal stress to water pressure and determine whether it should ground*/
+		if (sigma_nn<water_pressure) grounded=true;
+		else                         grounded=false;
+	}
+	else{
+		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);
@@ -1343,6 +1385,6 @@
 
 	/*clean up*/
+	delete gauss;
 	xDelete<IssmDouble>(xyz_list);
-
 }
 /*}}}*/