Changeset 26238

Timestamp:

05/04/21 14:56:31 (4 years ago)

Author:

caronlam

Message:

BUG: typo in GErotm1,2,3 and pole singularity in function SealevelchangeGeometryInitial

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -6142 +6142 @@
         IssmDouble dY21cos_dlon,dY21sin_dlon;
         IssmDouble LoveRotRSL,LoveRotU,LoveRothoriz;
+        IssmDouble polenudge;
         /*}}}*/
         /*Recover parameters:{{{ */
@@ -6277 +6278 @@
 
                 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 +6300 @@
                         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;
                                 }
                         }