Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26237)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26238)
@@ -6142,4 +6142,5 @@
 	IssmDouble dY21cos_dlon,dY21sin_dlon;
 	IssmDouble LoveRotRSL,LoveRotU,LoveRothoriz;
+	IssmDouble polenudge;
 	/*}}}*/
 	/*Recover parameters:{{{ */
@@ -6277,5 +6278,19 @@
 
 		for(int i=0;i<3;i++){
-			lati=latitude[i];
+
+			//Avoid singularity of the poles
+			if (latitude[i]==M_PI/2.){
+				//North pole: nudge south
+				polenudge=-1.e12;
+			}
+			else if (latitude[i]==-M_PI/2.){
+				//South pole: nudge north
+				polenudge=1.e12;
+			}
+			else {
+				polenudge=0.0;
+			}
+
+			lati=latitude[i]+polenudge;
 			longi=longitude[i];
 
@@ -6285,28 +6300,29 @@
 			Y20   = -(1.0/6.0 - 0.5*cos(2.0*lati));
 
-			Grotm1[i]=  LoveRotRSL*Y21cos;
-			Grotm2[i]=  LoveRotRSL*Y21sin;
-			Grotm3[i]=  LoveRotRSL*Y20;
+			Grotm1[i]= LoveRotRSL*Y21cos;
+			Grotm2[i]= LoveRotRSL*Y21sin;
+			Grotm3[i]= LoveRotRSL*Y20;
 
 			if (computeelastic){
-				GUrotm1[i]=  LoveRotU*Y21cos;
-				GUrotm2[i]=  LoveRotU*Y21sin;
-				GUrotm3[i]=  LoveRotU*Y20;
+				GUrotm1[i]= LoveRotU*Y21cos;
+				GUrotm2[i]= LoveRotU*Y21sin;
+				GUrotm3[i]= LoveRotU*Y20;
 				if (horiz){
 				//bed_N = Love_l * d(Y)/dlat ;
 				dY21cos_dlat=-cos(2.*lati)*cos(longi);
 				dY21sin_dlat=-cos(2.*lati)*sin(longi);
-				dY20_dlat=-sin(2.0*lati);
-				GNrotm1[i]=  LoveRothoriz*dY21cos_dlat;
-				GNrotm2[i]=  LoveRothoriz*dY21sin_dlat;
-				GNrotm3[i]=  LoveRothoriz*dY20_dlat;
+				dY20_dlat=   -sin(2.*lati);
+				GNrotm1[i]= LoveRothoriz*dY21cos_dlat;
+				GNrotm2[i]= LoveRothoriz*dY21sin_dlat;
+				GNrotm3[i]= LoveRothoriz*dY20_dlat;
 
 				//bed_E = Love_l * 1/cos(lat) * d(Y)/dlon ;
-				dY21cos_dlon=-Y21sin;
-				dY21sin_dlon=Y21cos;
+				dY21cos_dlon=-Y21sin/cos(lati);
+				dY21sin_dlon=Y21cos/cos(lati);
 				//dY20_dlon=0.;
-				GNrotm1[i]=  LoveRothoriz*dY21cos_dlon/cos(lati);
-				GNrotm2[i]=  LoveRothoriz*dY21sin_dlon/cos(lati);
-				GNrotm3[i]=  0.0;
+
+				GErotm1[i]= LoveRothoriz*dY21cos_dlon;
+				GErotm2[i]= LoveRothoriz*dY21sin_dlon;
+				GErotm3[i]= 0.0;
 				}
 			}