Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Tria.cpp	(revision 21374)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Tria.cpp	(revision 21375)
@@ -1066,5 +1066,5 @@
 
 	}
-	else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum){
+	else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
 		/*Check that not all nodes are grounded or floating*/
 		if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
@@ -3622,4 +3622,5 @@
 	IssmDouble late,longe,re;
 	IssmDouble lati,longi,ri;
+	bool notfullygrounded=false;
 
 	/*elastic green function:*/
@@ -3642,4 +3643,15 @@
 		*peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
 		return;
+	}
+
+	/*early return if we are fully floating: */
+	if (this->inputs->Max(MaskGroundediceLevelsetEnum)<=0){
+		*peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
+		return;
+	}
+	
+	/*If we are here, we are on ice that is fully grounded or half-way to floating: */
+	if ((this->inputs->Min(MaskGroundediceLevelsetEnum))<0){
+		notfullygrounded=true; //used later on.
 	}
 
@@ -3703,14 +3715,44 @@
 	/*}}}*/
 
-	/*Compute area of element:*/
+	/*Compute area of element. Scale it by grounded fraction if not fully grounded: */
 	area=GetAreaSpherical();
-
+	if(notfullygrounded){
+		IssmDouble phi=0;
+		IssmDouble xyz_list[NUMVERTICES][3];
+		::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+		phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias
+		area*=phi;
+	}
+	
 	/*Compute ice thickness change: */
 	Input*	deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); 
 	if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
-	deltathickness_input->GetInputAverage(&I);
-
-	/*Are we on an ice shelf? if so, does not contribute to sea level rise!: */
-	if (this->IsFloating()) I=0;
+
+	/*If we are fully grounded, take the average over the element: */
+	if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
+	else{ 
+		IssmDouble total_weight=0;
+		bool mainlyfloating = true;
+		int         point1;
+		IssmDouble  fraction1,fraction2;
+		
+		/*Recover portion of element that is grounded*/
+		this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+		Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
+
+		/* Start  looping on the number of gaussian points and average over these gaussian points: */
+		total_weight=0;
+		I=0;
+		for(int ig=gauss->begin();ig<gauss->end();ig++){
+			IssmDouble Ig=0;
+			gauss->GaussPoint(ig);
+			deltathickness_input->GetInputValue(&Ig);
+			I+=Ig*gauss->weight;
+			total_weight+=gauss->weight;
+		}
+		I=I/total_weight;
+		delete gauss;
+	}
 
 	/*Compute eustatic compoent:*/