Changeset 26296

Timestamp:

06/04/21 16:35:40 (4 years ago)

Author:

caronlam

Message:

BUG: BarystaticLoad was considering coastline elements that were fully covered by ice as also fully grounded; CHG: lots of renaming in sealevelchange_core & analysis; CHG: it is now possible to compute GRD with Barystatic contributions & ocean mass but without sea-level loading

Location:

issm/trunk-jpl

Files:

• src/c/analyses/EsaAnalysis.cpp (modified) (1 diff)
• src/c/analyses/SealevelchangeAnalysis.cpp (modified) (18 diffs)
• src/c/classes/Elements/Element.h (modified) (2 diffs)
• src/c/classes/Elements/Penta.h (modified) (2 diffs)
• src/c/classes/Elements/Seg.h (modified) (2 diffs)
• src/c/classes/Elements/Tetra.h (modified) (2 diffs)
• src/c/classes/Elements/Tria.cpp (modified) (63 diffs)
• src/c/classes/Elements/Tria.h (modified) (2 diffs)
• src/c/cores/sealevelchange_core.cpp (modified) (17 diffs)
• src/c/shared/Enum/EnumDefinitions.h (modified) (7 diffs)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (4 diffs)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (4 diffs)
• src/m/classes/solidearthsettings.m (modified) (9 diffs)
• test/NightlyRun/test2001.m (modified) (1 diff)
• test/NightlyRun/test2002.m (modified) (6 diffs)
• test/NightlyRun/test2003.m (modified) (6 diffs)
• test/NightlyRun/test2004.m (modified) (5 diffs)
• test/NightlyRun/test2005.m (modified) (3 diffs)
• test/NightlyRun/test2006.m (modified) (3 diffs)
• test/NightlyRun/test2007.m (modified) (2 diffs)
• test/NightlyRun/test2008.m (modified) (3 diffs)
• test/NightlyRun/test2010.m (modified) (2 diffs)
• test/NightlyRun/test2011.m (modified) (6 diffs)
• test/NightlyRun/test2090.m (modified) (5 diffs)

issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp

@@ -61 +61 @@
         /*some constant parameters: */
         parameters->AddObject(iomodel->CopyConstantObject("md.esa.hemisphere",EsaHemisphereEnum));
-        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rigid",SolidearthSettingsRigidEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.selfattraction",SolidearthSettingsSelfAttractionEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum));

issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp

@@ -92 +92 @@
         int         isexternal=0;
 
-        IssmDouble* G_rigid = NULL;
-        IssmDouble* G_rigid_local = NULL;
+        IssmDouble* G_gravi = NULL;
+        IssmDouble* G_gravi_local = NULL;
         IssmDouble* G_viscoelastic = NULL;
         IssmDouble* G_viscoelastic_interpolated = NULL;
@@ -108 +108 @@
         IssmDouble  planetradius=0;
         IssmDouble  planetarea=0;
-        bool            rigid=false;
+        bool            selfattraction=false;
         bool            elastic=false;
         bool            viscous=false;
@@ -145 +145 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.abstol",SolidearthSettingsAbstolEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.maxiter",SolidearthSettingsMaxiterEnum));
-        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rigid",SolidearthSettingsRigidEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.selfattraction",SolidearthSettingsSelfAttractionEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum));
@@ -153 +153 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.polarmoi",RotationalPolarMoiEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.angularvelocity",RotationalAngularVelocityEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.grdocean",SolidearthSettingsGrdOceanEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.ocean_area_scaling",SolidearthSettingsOceanAreaScalingEnum));
-        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.sealevelloading",SolidearthSettingsSealevelLoadingEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.isgrd",SolidearthSettingsGRDEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.compute_bp_grd",SolidearthSettingsComputeBpGrdEnum));
@@ -228 +229 @@
 
         parameters->FindParam(&grdmodel,GrdModelEnum);
-        if(grdmodel!=IvinsEnum){
+        if(grdmodel==ElasticEnum){
                 /*Deal with elasticity {{{*/
-                iomodel->FetchData(&rigid,"md.solidearth.settings.rigid");
+                iomodel->FetchData(&selfattraction,"md.solidearth.settings.selfattraction");
                 iomodel->FetchData(&elastic,"md.solidearth.settings.elastic");
                 iomodel->FetchData(&viscous,"md.solidearth.settings.viscous");
@@ -236 +237 @@
                 iomodel->FetchData(&horiz,"md.solidearth.settings.horiz");
 
-                if(rigid){
+                if(selfattraction){
                         /*compute green functions for a range of angles*/
                         iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
@@ -271 +272 @@
                         // // Providing "t" will cause ensurecontiguous to be called.
                         if(viscous){
-                                IssmDouble* stacktimes=NULL;
+                                IssmDouble* viscoustimes=NULL;
                                 ntimesteps=precomputednt; 
                                 nt=reCast<int,IssmDouble>((final_time-start_time)/timeacc)+1;
 
-                                parameters->AddObject(new IntParam(StackNumStepsEnum,nt));
-                                /*Initialize stack times:*/
-                                stacktimes=xNew<IssmDouble>(nt);
+                                parameters->AddObject(new IntParam(SealevelchangeViscousNumStepsEnum,nt));
+                                /*Initialize viscous stack times:*/
+                                viscoustimes=xNew<IssmDouble>(nt);
                                 for(int t=0;t<nt;t++){
-                                        stacktimes[t]=start_time+timeacc*t;
+                                        viscoustimes[t]=start_time+timeacc*t;
                                 }
-                                parameters->AddObject(new DoubleVecParam(StackTimesEnum,stacktimes,nt));
-                                parameters->AddObject(new IntParam(StackIndexEnum,0));
-                                xDelete<IssmDouble>(stacktimes);
+                                parameters->AddObject(new DoubleVecParam(SealevelchangeViscousTimesEnum,viscoustimes,nt));
+                                parameters->AddObject(new IntParam(SealevelchangeViscousIndexEnum,0));
+                                xDelete<IssmDouble>(viscoustimes);
                         }
                         else {
@@ -300 +301 @@
 #endif
                 }
-                if(rigid){
+                if(selfattraction){
 #ifdef _HAVE_AD_
-                        G_rigid=xNew<IssmDouble>(M,"t");
+                        G_gravi=xNew<IssmDouble>(M,"t");
 #else
-                        G_rigid=xNew<IssmDouble>(M);
+                        G_gravi=xNew<IssmDouble>(M);
 #endif
                 }
@@ -310 +311 @@
                 if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
 
-                if(rigid){
+                if(selfattraction){
 
                         /*compute combined legendre + love number (elastic green function:*/
@@ -327 +328 @@
 #endif
                 }
-                if(rigid){
+                if(selfattraction){
 #ifdef _HAVE_AD_
-                        G_rigid_local=xNew<IssmDouble>(m,"t");
+                        G_gravi_local=xNew<IssmDouble>(m,"t");
 #else
-                        G_rigid_local=xNew<IssmDouble>(m);
+                        G_gravi_local=xNew<IssmDouble>(m);
 #endif
                 }
 
-                if(rigid){
+                if(selfattraction){
                         for(int i=lower_row;i<upper_row;i++){
                                 IssmDouble alpha,x;
                                 alpha= reCast<IssmDouble>(i)*degacc * M_PI / 180.0;
-                                G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
+                                G_gravi_local[i-lower_row]= .5/sin(alpha/2.0);
                         }
                 }
@@ -348 +349 @@
 
                                 for(int t=0;t<ntimesteps;t++){
-                                        G_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (love_k[(ndeg-1)*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t])*G_rigid_local[i-lower_row];
-                                        U_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (love_h[(ndeg-1)*precomputednt+t])*G_rigid_local[i-lower_row];
+                                        G_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (love_k[(ndeg-1)*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t])*G_gravi_local[i-lower_row];
+                                        U_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (love_h[(ndeg-1)*precomputednt+t])*G_gravi_local[i-lower_row];
                                         if(horiz)H_viscoelastic_local[(i-lower_row)*ntimesteps+t]= 0; 
                                 }
@@ -392 +393 @@
                         }
                 }
-                if(rigid){
+                if(selfattraction){
 
                         /*merge G_viscoelastic_local into G_viscoelastic; U_viscoelastic_local into U_viscoelastic; H_viscoelastic_local to H_viscoelastic:{{{*/
@@ -400 +401 @@
                         int  offset;
 
-                        //deal with rigid first: 
+                        //deal with selfattraction first: 
                         ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
 
@@ -407 +408 @@
 
                         /*All gather:*/
-                        ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+                        ISSM_MPI_Allgatherv(G_gravi_local, m, ISSM_MPI_DOUBLE, G_gravi, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
 
                         if(elastic){
@@ -427 +428 @@
 
                         /*Avoid singularity at 0: */
-                        G_rigid[0]=G_rigid[1];
+                        G_gravi[0]=G_gravi[1];
                         if(elastic){
                                 for(int t=0;t<ntimesteps;t++){
@@ -486 +487 @@
                         }
                         else if(elastic){
-                                nt=1; //in elastic, or if we run only rigid, we need only one step
+                                nt=1; //in elastic, or if we run only selfattraction, we need only one step
 #ifdef _HAVE_AD_
                                 G_viscoelastic_interpolated=G_viscoelastic;
@@ -499 +500 @@
 
                         /*Save our precomputed tables into parameters*/
-                        parameters->AddObject(new DoubleVecParam(SealevelchangeGRigidEnum,G_rigid,M));
+                        parameters->AddObject(new DoubleVecParam(SealevelchangeGSelfAttractionEnum,G_gravi,M));
                         if(viscous || elastic){
                                 parameters->AddObject(new DoubleVecParam(SealevelchangeGViscoElasticEnum,G_viscoelastic_interpolated,M*nt));
@@ -507 +508 @@
 
                         /*free resources: */
-                        xDelete<IssmDouble>(G_rigid);
-                        xDelete<IssmDouble>(G_rigid_local);
+                        xDelete<IssmDouble>(G_gravi);
+                        xDelete<IssmDouble>(G_gravi_local);
                         if(elastic){
                                 xDelete<IssmDouble>(love_h);

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -391 +391 @@
                 virtual void          GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0;
 
-                virtual void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom)=0;
+                virtual void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom)=0;
                 virtual void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0;
                 virtual void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom)=0;
@@ -400 +400 @@
                 virtual void       SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
                 virtual void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom)=0;
-                virtual void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
+                virtual void       SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
                 virtual void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0;
                 virtual void       SealevelchangeUpdateViscousFields()=0;

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -223 +223 @@
                 #ifdef _HAVE_SEALEVELCHANGE_
                 void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
-                void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
+                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads,  SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
@@ -232 +232 @@
                 void       SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");};
-                void       SealevelchangeDeformationConvolution(GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
+                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -178 +178 @@
 #ifdef _HAVE_SEALEVELCHANGE_
                 void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
-                void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
+                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads,  SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
@@ -187 +187 @@
                 void       SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");};
-                void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
+                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/Elements/Tetra.h

@@ -185 +185 @@
 #ifdef _HAVE_SEALEVELCHANGE_
                 void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
-                void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
+                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads,  SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
@@ -194 +194 @@
                 void       SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");};
-                void       SealevelchangeDeformationConvolution(GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
+                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
                 void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};

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

@@ -1937 +1937 @@
 
         bool istrapeze1, istrapeze2;
-        IssmDouble phi1,phi2, f11,f12,f21,f22;
+        IssmDouble phi1,phi2, d,e,f,g,h1,h2;
         int point1, point2,  i0,i1,i2,j0,j1,j2;
         IssmDouble weights1[3],weights2[3];
@@ -1951 +1951 @@
         IssmDouble xyz3[3][3]={0};
         IssmDouble xyz[8][3]={0};
-        IssmDouble f1o;
         IssmDouble w[8][NUMVERTICES]={0};
         IssmDouble areasub=0;
@@ -2019 +2018 @@
         //If just one levelset is all negative, just take the partitioning of the other, no interaction between them
         if (levelset1[0]<=0 && levelset1[1]<=0 && levelset1[2]<=0) {
-                this->GetFractionGeometry(loadweights,&phi2,&point2,&f21,&f22,&istrapeze2,levelset2);
-                this->GetBarycenterFromLevelset(platbar,plongbar, phi2, f21, f22, late, longe, point2,istrapeze2,planetradius);
+                this->GetFractionGeometry(loadweights,&phi2,&point2,&f,&g,&istrapeze2,levelset2);
+                this->GetBarycenterFromLevelset(platbar,plongbar, phi2, f, g, late, longe, point2,istrapeze2,planetradius);
                 *ploadarea=area*phi2;
                 return;
         }
         if (levelset2[0]<=0 && levelset2[1]<=0 && levelset2[2]<=0) {
-                this->GetFractionGeometry(loadweights,&phi1,&point1,&f11,&f12,&istrapeze1,levelset1);
-                this->GetBarycenterFromLevelset(platbar,plongbar, phi1, f11, f12, late, longe, point1,istrapeze1,planetradius);
+                this->GetFractionGeometry(loadweights,&phi1,&point1,&d,&e,&istrapeze1,levelset1);
+                this->GetBarycenterFromLevelset(platbar,plongbar, phi1, d, e, late, longe, point1,istrapeze1,planetradius);
                 *ploadarea=area*phi1;
                 return;
@@ -2032 +2031 @@
 
 
-        this->GetFractionGeometry(&weights1[0],&phi1,&point1,&f11,&f12,&istrapeze1,levelset1);
-        this->GetFractionGeometry(&weights2[0],&phi2,&point2,&f21,&f22,&istrapeze2,levelset2);
+        this->GetFractionGeometry(&weights1[0],&phi1,&point1,&d,&e,&istrapeze1,levelset1);
+        this->GetFractionGeometry(&weights2[0],&phi2,&point2,&f,&g,&istrapeze2,levelset2);
 
 
@@ -2039 +2038 @@
         if (istrapeze1==istrapeze2 && point1==point2 && phi1==phi2){
                 //the two levelsets are redundant: levelset1 = positivescalar * levelset2
-                this->GetBarycenterFromLevelset(platbar,plongbar, phi1, f11, f12, late, longe, point1,istrapeze1,planetradius);
+                this->GetBarycenterFromLevelset(platbar,plongbar, phi1, d, e, late, longe, point1,istrapeze1,planetradius);
                 *ploadarea=area*phi1;
                 for (int i=0;i<NUMVERTICES;i++) loadweights[i]=weights1[i];
@@ -2056 +2055 @@
         ::GetVerticesCoordinates(&xyz0[0][0],vertices,NUMVERTICES); // initial triangle
 
-        i0=point1;
-        i1=(point1+1)%3;
-        i2=(point1+2)%3;
-
-        j0=point2;
-        j1=(point2+1)%3;
-        j2=(point2+2)%3;
-
-        //f1o: fraction of segment {point_f11 -> point_f12} where the levelsets intersect (only used when f1o>=0 and f1o<=1)
-        if(point2==i0) f1o= f22*(f11-f21)/(f11*f22-f12*f21);
-        else if(point2==i1) f1o=f21*(1.0-f22-f11)/(f21*(f12-f11)-f12*f22);
-        else f1o= (f22*(1.0-f21-f11)+f21*f11)/(f22*(f12-f11) +f21*f11);
+        //Let our element be triangle ABC with:
+        i0=point1; //A
+        i1=(point1+1)%3; //B
+        i2=(point1+2)%3; //C
+
+        j0=point2; //Can be A, B or C
+        j1=(point2+1)%3; //anticlockwise point from j0
+        j2=(point2+2)%3; //clockwise point from j0
+
+        /* Below we define the relative fractional lengths of ABC where the zero-level contours of the two level sets intersect with ABC and each other. For example D is the intersection of level set 1 with side AB with fractional length d=[AD]/[AB].
+
+           levelset1 intersects ABC on D and E:
+           A------D---B
+           <--d--->             
+
+           A-----E----C
+           <--e-->
+
+           levelset2 intersects ABC on F and G:
+           j0---F------j1
+           <--f->
+
+           j0-------G--j2
+           <---g---->
+
+           levelset1 and 2 intersect on H (when that intersection exists inside the element)
+           D----H------E
+           <-h1->
+
+           F-----H----G
+           <--h2->
+        */
+
+        if (point2==i0){
+                h1= g*(d-f)/(d*g-e*f);
+                h2= e/g * h1;
+        }
+        else if (point2==i1){
+                h1=f*(1.0-g-d)/(f*(e-d)-e*g);
+                h2= 1.0-e/f * h1;
+        }
+        else if (point2==i2){
+                h1= (g*(1.0-f-d)+f*d)/(g*(e-d) +f*d); 
+                h2= (d*(f+e-1))/(g*(e-d) +f*d);
+        }
 
 
         //interpolant weights of each point. Any field F[0,1,2] provided at the original vertices [0,1,2] will be equal on point k to sum_i (F[i] * w[k][i])
-        w[0][0]=1;
-        w[1][1]=1;
-        w[2][2]=1;
-        w[3][i0]=1.0-f11; w[3][i1]=f11;
-        w[4][i0]=1.0-f12; w[4][i2]=f12;
-        w[5][j0]=1.0-f21; w[5][j1]=f21;
-        w[6][j0]=1.0-f22; w[6][j2]=f22;
-        for (int j=0;j<3;j++) w[7][j]=w[3][j]*(1.0-f1o) + w[4][j]*f1o; //we interpolate the intersection point between point_f11 and point_f12 at fraction f1o
+        w[0][0]=1; //A
+        w[1][1]=1; //B
+        w[2][2]=1; //C
+        w[3][i0]=1.0-d; w[3][i1]=d; //D
+        w[4][i0]=1.0-e; w[4][i2]=e; //E
+        w[5][j0]=1.0-f; w[5][j1]=f; //F
+        w[6][j0]=1.0-g; w[6][j2]=g; //G
+        for (int j=0;j<3;j++) w[7][j]=w[3][j]*(1.0-h1) + w[4][j]*h1; //H: we interpolate the intersection point H between D and E at fraction h1
 
         for (int k=0;k<8;k++){
@@ -2088 +2120 @@
                 //point2 can be either i0,i1 or i2. We start the search with i1 and i2 as they have less computational cost in ifs
                 if(point2==i2){ /*{{{*/
-                        if (f12>1.0-f21){ /*{{{*/
+                        if (e>1.0-f){ /*{{{*/
                                 if (!istrapeze1 && !istrapeze2){
-                                        tria1[0]=5; tria1[1]= 7; tria1[2]= 4;
+                                        tria1[0]=5; tria1[1]= 7; tria1[2]= 4; // FHE
+                                        area1=h2*(e+f-1.0)*g;
                                 }
                                 else if (!istrapeze1 && istrapeze2){
-                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 7;
-                                        tria2[0]=i0; tria2[1]= 7; tria2[2]= 5;
+                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 5; //ADF
+                                        area1=d*(1.0-f);
+                                        tria2[0]=3; tria2[1]= 7; tria2[2]= 5; //DHF
+                                        area2=d*h1*(e+f-1.0);
                                 }
                                 else if (istrapeze1 && !istrapeze2){
-                                        tria1[0]=7; tria1[1]= 6; tria1[2]= 4;
-                                        tria2[0]=4; tria2[1]= 6; tria2[2]= i2;
+                                        tria1[0]=7; tria1[1]= 6; tria1[2]= 4; //HGE
+                                        area1=g*(1.0-h2)*(e+f-1.0);
+                                        tria2[0]=4; tria2[1]= 6; tria2[2]= i2; //EGC
+                                        area2=g*(1.0-e);
                                 }
                                 else { //istrapeze1 && istrapeze2
-                                        tria1[0]=3; tria1[1]= i1; tria1[2]= 7;
-                                        tria2[0]=7; tria2[1]= i1; tria2[2]= 6;
+                                        tria1[0]=3; tria1[1]= i1; tria1[2]= 6; //DBG
+                                        area1=(1.0-d)*(1.0-g);
+                                        tria2[0]=3; tria2[1]= 6; tria2[2]= 7; //DGH
+                                        area2=g*((1.0-f)*(1.0-h2)+h2*e)+d*(1.0-e-g);
                                 }  /*}}}*/
                         }
-                        else if (f12<=1.0-f21){ /*{{{*/
+                        else if (e<=1.0-f){ /*{{{*/
                                 if (!istrapeze1 && !istrapeze2){
                                 }
                                 else if (!istrapeze1 && istrapeze2){
-                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 4;
+                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 4; //ADE
+                                        area1=d*e;
                                 }
                                 else if (istrapeze1 && !istrapeze2){
-                                        tria1[0]=5; tria1[1]= 6; tria1[2]= i2;
+                                        tria1[0]=5; tria1[1]= 6; tria1[2]= i2; //FGC
+                                        area1=f*g;
                                 }
                                 else { //istrapeze1 && istrapeze2
-                                        tria1[0]=3; tria1[1]= i1; tria1[2]= 4;
-                                        tria2[0]=4; tria2[1]= i1; tria2[2]= 5;
-                                        tria3[0]=5; tria3[1]= i1; tria3[2]= 6;
+                                        tria1[0]=3; tria1[1]= i1; tria1[2]= 5; //DBF
+                                        area1=(1.0-d)*(1.0-f);
+                                        tria2[0]=4; tria2[1]= 3; tria2[2]= 5;  //EDF
+                                        area2=d*(1.0-e-f);
+                                        tria3[0]=5; tria3[1]= i1; tria3[2]= 6; //FBG
+                                        area3=f*(1.0-g);
                                 } /*}}}*/
                         } 
                 }/*}}}*/    
                 else if(point2==i1){ /*{{{*/
-                        if (f11>1.0-f22){ /*{{{*/
+                        if (d>1.0-g){ /*{{{*/
                                 if (!istrapeze1 && !istrapeze2){
-                                        tria1[0]=6; tria1[1]= 3; tria1[2]= 7;
+                                        tria1[0]=6; tria1[1]= 3; tria1[2]= 7; //GDH
+                                        area1=(1.0-h2)*(d+g-1.0)*f;
                                 }
                                 else if (!istrapeze1 && istrapeze2){
-                                        tria1[0]=i0; tria1[1]= 6; tria1[2]= 7;
-                                        tria2[0]=i0; tria2[1]= 7; tria2[2]= 4;
+                                        tria1[0]=i0; tria1[1]= 6; tria1[2]= 4; //AGE
+                                        area1=(1.0-g)*e;
+                                        tria2[0]=6; tria2[1]= 7; tria2[2]= 4; //GHE
+                                        area2=e*(1.0-h1)*(g+d-1.0);
                                 }
                                 else if (istrapeze1 && !istrapeze2){
-                                        tria1[0]=3; tria1[1]= i1; tria1[2]= 7;
-                                        tria2[0]=7; tria2[1]= i1; tria2[2]= 5;
+                                        tria1[0]=i1; tria1[1]= 5; tria1[2]= 3; //BFD
+                                        area1=(1.0-d)*f;
+                                        tria2[0]=3; tria2[1]= 5; tria2[2]= 7; //DFH
+                                        area2=f*h2*(d+g-1.0);
                                 }
                                 else { //istrapeze1 && istrapeze2
-                                        tria1[0]=7; tria1[1]= 5; tria1[2]= 4;
-                                        tria2[0]=4; tria2[1]= 5; tria2[2]= i2;
+                                        tria1[0]=7; tria1[1]= 5; tria1[2]= 4; //HFE
+                                        area1=e*((1.0-d)*(1.0-h1)+h1*g)+f*(1.0-g-e);
+                                        tria2[0]=4; tria2[1]= 5; tria2[2]= i2;//EFC
+                                        area2=(1.0-e)*(1.0-f);
                                 }  /*}}}*/
                         }
-                        else if (f11<=1.0-f22){ /*{{{*/
+                        else if (d<=1.0-g){ /*{{{*/
                                 if (!istrapeze1 && !istrapeze2){
                                 }
                                 else if (!istrapeze1 && istrapeze2){
-                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 4;
+                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 4; //ADE
+                                        area1=d*e;
                                 }
                                 else if (istrapeze1 && !istrapeze2){
-                                        tria1[0]=6; tria1[1]= i1; tria1[2]= 5;
+                                        tria1[0]=6; tria1[1]= i1; tria1[2]= 5; //GBF
+                                        area1=f*g;
                                 }
                                 else { //istrapeze1 && istrapeze2
-                                        tria1[0]=3; tria1[1]= 6; tria1[2]= 4;
-                                        tria2[0]=4; tria2[1]= 6; tria2[2]= 5;
-                                        tria3[0]=4; tria3[1]= 5; tria3[2]= i2;
+                                        tria1[0]=3; tria1[1]= 6; tria1[2]= 4; //DGE
+                                        area1=e*(1.0-d-g);
+                                        tria2[0]=6; tria2[1]= i2; tria2[2]= 4; //GCE
+                                        area2=(1.0-g)*(1.0-e);
+                                        tria3[0]=6; tria3[1]= 5; tria3[2]= i2; //GFC
+                                        area3=g*(1.0-f);
                                 } /*}}}*/
                         }
@@ -2157 +2213 @@
                 }/*}}}*/
                 else{ /*{{{*/
-                        if (f11<=f21 && f12>=f22){  /*{{{*/
+                        if (d<=f && e>=g){  /*{{{*/
                                 if (!istrapeze1 && !istrapeze2){
-                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 7;
-                                        tria2[0]=i0; tria2[1]= 7; tria2[2]= 6;
+                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 7; //ADH
+                                        area1=h1*d*e;
+                                        tria2[0]=i0; tria2[1]= 7; tria2[2]= 6; //AHG
+                                        area2=(1.0-h2)*f*g;
                                 }
                                 else if (!istrapeze1 && istrapeze2){
-                                        tria1[0]=6; tria1[1]= 7; tria1[2]= 4;
+                                        tria1[0]=6; tria1[1]= 7; tria1[2]= 4; //GHE
+                                        area1=(e-g)*d*(1.0-h1);
                                 }
                                 else if (istrapeze1 && !istrapeze2){
-                                        tria1[0]=3; tria1[1]= 5; tria1[2]= 7;
+                                        tria1[0]=3; tria1[1]= 5; tria1[2]= 7; //DFH
+                                        area1=(f-d)*g*h2;
                                 }
                                 else { //istrapeze1 && istrapeze2
-                                        tria1[0]=7; tria1[1]= 5; tria1[2]= i1;
-                                        tria2[0]=7; tria2[1]= i1; tria2[2]= 4;
-                                        tria3[0]=4; tria3[1]= i1; tria3[2]= i2;
+                                        tria1[0]=5; tria1[1]= i1; tria1[2]= 7; //FBH
+                                        area1=g*h2*(1.0-f);
+                                        tria2[0]=7; tria2[1]= i1; tria2[2]= i2;//HBC
+                                        area2=1.0+d*(h1-e*h1-1.0)+g*h2*(d-1.0);
+                                        tria3[0]=7; tria3[1]= i2; tria3[2]= 4; //HCE
+                                        area3=d*(1.0-h1)*(1.0-e);
                                 } /*}}}*/
                         }
-                        else if (f11>=f21 && f12<=f22){ /*{{{*/
+                        else if (d>=f && e<=g){ /*{{{*/
                                 if (!istrapeze1 && !istrapeze2){
-                                        tria1[0]=i0; tria1[1]= 5; tria1[2]= 7;
-                                        tria2[0]=i0; tria2[1]= 4; tria2[2]= 7;
+                                        tria1[0]=i0; tria1[1]= 5; tria1[2]= 7; //AFH
+                                        area1=h2*f*g;
+                                        tria2[0]=i0; tria2[1]= 7; tria2[2]= 4; //AHE
+                                        area2=(1.0-h1)*d*e;
                                 }else if (!istrapeze1 && istrapeze2){
-                                        tria1[0]=5; tria1[1]= 3; tria1[2]= 7;
+                                        tria1[0]=5; tria1[1]= 3; tria1[2]= 7; //FDH
+                                        area1=(d-f)*e*h1;
                                 }else if (istrapeze1 && !istrapeze2){
-                                        tria1[0]=4; tria1[1]= 7; tria1[2]= 6;
+                                        tria1[0]=4; tria1[1]= 7; tria1[2]= 6; //EHG
+                                        area1=(g-e)*f*(1.0-h2);
                                 }else { //istrapeze1 && istrapeze2
-                                        tria1[0]=3; tria1[1]= i1; tria1[2]= 7;
-                                        tria2[0]=7; tria2[1]= i1; tria2[2]= 6;
-                                        tria3[0]=6; tria3[1]= i1; tria3[2]= i2;
+                                        tria1[0]=3; tria1[1]= i1; tria1[2]= 7; //DBH
+                                        area1=e*h1*(1.0-d);
+                                        tria2[0]=7; tria2[1]= i1; tria2[2]= 6; //HCG
+                                        area2=f*(1.0-h2)*(1.0-g);
+                                        tria3[0]=6; tria3[1]= i1; tria3[2]= i2; //HBC
+                                        area3=1.0+f*(h2-g*h2-1.0)+e*h1*(f-1.0);
                                 }  /*}}}*/
                         }
-                        else if (f11<=f21 && f12<=f22){ /*{{{*/
+                        else if (d<=f && e<=g){ /*{{{*/
                                 if (!istrapeze1 && !istrapeze2){
-                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 4;
+                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 4;//ADE
+                                        area1=d*e;
                                 }else if (!istrapeze1 && istrapeze2){
                                 }else if (istrapeze1 && !istrapeze2){
-                                        tria1[0]=3; tria1[1]= 5; tria1[2]= 4;
-                                        tria2[0]=4; tria2[1]= 5; tria2[2]= 6;
+                                        tria1[0]=3; tria1[1]= 5; tria1[2]= 6; //DFG
+                                        area1=g*(f-d);
+                                        tria2[0]=3; tria2[1]= 6; tria2[2]= 4; //DGE
+                                        area2=d*(g-e);
                                 }else { //istrapeze1 && istrapeze2
-                                        tria1[0]=6; tria1[1]= 5; tria1[2]= i1;
-                                        tria2[0]=6; tria2[1]= i1; tria2[2]= i2;
+                                        tria1[0]=5; tria1[1]= i2; tria1[2]= 6; //FCG
+                                        area1=f*(1.0-g);
+                                        tria2[0]=5; tria2[1]= i1; tria2[2]= i2;//FBC
+                                        area2=1.0-f;
                                 }  /*}}}*/
                         }
-                        else if (f11>=f21 && f12>=f22){ /*{{{*/
+                        else if (d>=f && e>=g){ /*{{{*/
                                 if (!istrapeze1 && !istrapeze2){
-                                        tria1[0]=i0; tria1[1]= 5; tria1[2]= 6;
+                                        tria1[0]=i0; tria1[1]= 5; tria1[2]= 6; //AFG
+                                        area1=f*g;
                                 }else if (!istrapeze1 && istrapeze2){
-                                        tria1[0]=5; tria1[1]= 3; tria1[2]= 6;
-                                        tria2[0]=6; tria2[1]= 3; tria2[2]= 4;
+                                        tria1[0]=5; tria1[1]= 4; tria1[2]= 6; //FEG
+                                        area1=f*(e-g);
+                                        tria2[0]=5; tria2[1]= 3; tria2[2]= 4; //FDE
+                                        area2=e*(d-f);
                                 }else if (istrapeze1 && !istrapeze2){
                                 }else { //istrapeze1 && istrapeze2
-                                        tria1[0]=4; tria1[1]= 3; tria1[2]= i1;
-                                        tria2[0]=i1; tria2[1]= i2; tria2[2]= 4;
+                                        tria1[0]=3; tria1[1]= i1; tria1[2]= i2; //DBC
+                                        area1=1.0-d;
+                                        tria2[0]=3; tria2[1]= i2; tria2[2]= 4;//DCE
+                                        area2=d*(1.0-e);
                                 }  /*}}}*/
                         } 
@@ -2221 +2301 @@
                         }
                 }
-                area1= GetTriangleAreaSpherical(xyz1);
+                area1*=area; //dimensionalize the fractional area from [0:1] to [0:area]
         }
         if(tria2[0]>-1){ 
@@ -2230 +2310 @@
                         }
                 }
-                area2= GetTriangleAreaSpherical(xyz2);
+                area2*=area;
         }
         if(tria3[0]>-1){ 
@@ -2239 +2319 @@
                         }
                 }
-                area3= GetTriangleAreaSpherical(xyz3);
+                area3*=area;
         }
 
@@ -6100 +6180 @@
         IssmDouble xyz_list[NUMVERTICES][3];
 
-        /*stacks:*/
-        IssmDouble* stackRSL = NULL;
-        IssmDouble* stackU = NULL;
-        IssmDouble* stackN = NULL;
-        IssmDouble* stackE = NULL;
+        /*viscous stacks:*/
+        IssmDouble* viscousRSL = NULL;
+        IssmDouble* viscousU = NULL;
+        IssmDouble* viscousN = NULL;
+        IssmDouble* viscousE = NULL;
 
         #ifdef _HAVE_RESTRICT_
@@ -6112 +6192 @@
         IssmDouble* __restrict__ GE=NULL;
         IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
-        IssmDouble* __restrict__ G_rigid_precomputed=NULL;
+        IssmDouble* __restrict__ G_gravi_precomputed=NULL;
         IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
         IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
@@ -6121 +6201 @@
         IssmDouble* GE=NULL;
         IssmDouble* G_viscoelastic_precomputed=NULL;
-        IssmDouble* G_rigid_precomputed=NULL;
+        IssmDouble* G_gravi_precomputed=NULL;
         IssmDouble* U_viscoelastic_precomputed=NULL;
         IssmDouble* H_viscoelastic_precomputed=NULL;
@@ -6128 +6208 @@
         /*viscoelastic green function:*/
         int index;
-        int         M;
+        int M;
+        IssmDouble doubleindex,lincoef;
 
         /*Computational flags:*/
-        bool computerigid = false;
+        bool computeselfattraction = false;
         bool computeelastic = false;
         bool computerotation = false;
@@ -6140 +6221 @@
         IssmDouble start_time,final_time;
         int  nt,precomputednt;
+        int grd, grdmodel;
 
         /*Rotational:*/
@@ -6155 +6237 @@
         IssmDouble polenudge;
         /*}}}*/
+
         /*Recover parameters:{{{ */
         rho_earth=FindParam(MaterialsEarthDensityEnum);
-        this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
+        this->parameters->FindParam(&computeselfattraction,SolidearthSettingsSelfAttractionEnum);
         this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
         this->parameters->FindParam(&computerotation,SolidearthSettingsRotationEnum);
@@ -6166 +6249 @@
         this->parameters->FindParam(&NewtonG,ConstantsNewtonGravityEnum);
         this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
+        this->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); 
+        this->parameters->FindParam(&grdmodel,GrdModelEnum);
+
+        /*early return:*/
+        if (!grd || grdmodel!=ElasticEnum) return; //Love numbers won't be found in this case, return before loading them
+        if(!computeselfattraction)return;
 
         if(computerotation){
@@ -6177 +6266 @@
         /*}}}*/
 
-        /*early return:*/
-        if(!computerigid)return;
-
         /*Recover precomputed green function kernels:{{{*/
-        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGRigidEnum)); _assert_(parameter);
-        parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
+        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGSelfAttractionEnum)); _assert_(parameter);
+        parameter->GetParameterValueByPointer((IssmDouble**)&G_gravi_precomputed,&M);
 
         if(computeelastic){
@@ -6215 +6301 @@
         }
         else{
-                nt=1; //in elastic, or if we run only rigid, we need only one step
+                nt=1; //in elastic, or if we run only selfattraction, we need only one step
         }
 
@@ -6245 +6331 @@
                         delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
                         alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
-                        index=reCast<int,IssmDouble>( alpha/M_PI*reCast<IssmDouble,int>(M-1) );
+                        doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
+                        index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
                         _assert_(index>=0 && index<M);
+
+
+                        lincoef=doubleindex-index; //where between index and index+1 is alpha
+                        if (index==M-1){ //avoids out of bound case
+                                index-=1;
+                                lincoef=1;
+                        }
+                        
 
                         if(horiz){
@@ -6266 +6361 @@
 
                                 /*Rigid earth gravitational perturbation: */
-                                G[timeindex]+=G_rigid_precomputed[index];
+                                G[timeindex]+=G_gravi_precomputed[index+0]*(1.-lincoef) 
+                                             +G_gravi_precomputed[index+1]*lincoef; //linear interpolation
 
                                 /*Find a way to interpolate precomputed Gkernels to our solution time stepping:*/
                                 if(computeelastic){
-                                        G[timeindex]+=G_viscoelastic_precomputed[index*nt+t];
+                                        G[timeindex]+=G_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
+                                                     +G_viscoelastic_precomputed[(index+1)*nt+t]*lincoef; //linear interpolation
                                 }
                                 G[timeindex]=G[timeindex]*ratioe;
@@ -6276 +6373 @@
                                 /*Elastic components:*/
                                 if(computeelastic){
-                                        GU[timeindex] =  ratioe * U_viscoelastic_precomputed[index*nt+t];
+                                        GU[timeindex] =  ratioe *(U_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
+                                                                 +U_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
                                         if(horiz){
-                                                GN[timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*N_azim;
-                                                GE[timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*E_azim;
+                                                GN[timeindex] = N_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
+                                                                               +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
+                                                GE[timeindex] = E_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
+                                                                               +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
                                         }
                                 }
@@ -6376 +6476 @@
         }
         /*}}}*/
-        /*Initialize stacks: {{{*/
+        /*Initialize viscous stacks: {{{*/
         if(computeviscous){
-                stackRSL=xNewZeroInit<IssmDouble>(3*nt);
-                stackU=xNewZeroInit<IssmDouble>(3*nt);
-
-                this->inputs->SetArrayInput(StackRSLEnum,this->lid,stackRSL,3*nt);
-                this->inputs->SetArrayInput(StackUEnum,this->lid,stackU,3*nt);
-                this->parameters->SetParam(0,StackIndexEnum);
+                viscousRSL=xNewZeroInit<IssmDouble>(3*nt);
+                viscousU=xNewZeroInit<IssmDouble>(3*nt);
+
+                this->inputs->SetArrayInput(SealevelchangeViscousRSLEnum,this->lid,viscousRSL,3*nt);
+                this->inputs->SetArrayInput(SealevelchangeViscousUEnum,this->lid,viscousU,3*nt);
+                this->parameters->SetParam(0,SealevelchangeViscousIndexEnum);
                 if(horiz){
-                        stackN=xNewZeroInit<IssmDouble>(3*nt);
-                        stackE=xNewZeroInit<IssmDouble>(3*nt);
-                        this->inputs->SetArrayInput(StackNEnum,this->lid,stackRSL,3*nt);
-                        this->inputs->SetArrayInput(StackEEnum,this->lid,stackU,3*nt);
+                        viscousN=xNewZeroInit<IssmDouble>(3*nt);
+                        viscousE=xNewZeroInit<IssmDouble>(3*nt);
+                        this->inputs->SetArrayInput(SealevelchangeViscousNEnum,this->lid,viscousRSL,3*nt);
+                        this->inputs->SetArrayInput(SealevelchangeViscousEEnum,this->lid,viscousU,3*nt);
                 }
         }
@@ -6472 +6572 @@
         slgeom->late[this->lid]=late;
 
-        /*compute fractional areas and load weights for ocean:*/
+        /*compute areas and load weights for ocean and flag elements only partially in the ocean:*/
         if(isoceanonly){ 
                 slgeom->LoadArea[SLGEOM_OCEAN][this->lid]=area;
@@ -6582 +6682 @@
         /*Deal with ice loads if we are on grounded ice:*/
         if(isice && !isoceanonly && computeice){
-                if(isiceonly){
+                if(isiceonly && !isocean){
                         slgeom->LoadArea[SLGEOM_ICE][this->lid]=area;
                         for(int i=0;i<NUMVERTICES;i++) slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid]=1.0/3.0;
@@ -6635 +6735 @@
         IssmDouble loadweightsocean[3]; //to keep memory of these loads, no need to recompute for bottom pressure.
         IssmDouble xyz_list[NUMVERTICES][3];
-        IssmDouble latbar=0;
-        IssmDouble longbar=0;
+        IssmDouble latbar=slgeom->late[this->lid];
+        IssmDouble longbar=slgeom->longe[this->lid];
         IssmDouble constant;
         IssmDouble nanconstant=NAN;
@@ -6645 +6745 @@
 
         #ifdef _ISSM_DEBUG_
-        this->AddInput(SealevelBarystaticIceLatbarEnum,&nanconstant,P0Enum); 
-        this->AddInput(SealevelBarystaticIceLongbarEnum,&nanconstant,P0Enum); 
-        this->AddInput(SealevelBarystaticHydroLatbarEnum,&nanconstant,P0Enum); 
-        this->AddInput(SealevelBarystaticHydroLongbarEnum,&nanconstant,P0Enum); 
-        this->AddInput(SealevelBarystaticOceanLatbarEnum,&nanconstant,P0Enum); 
-        this->AddInput(SealevelBarystaticOceanLongbarEnum,&nanconstant,P0Enum); 
+        this->AddInput(SealevelBarystaticIceLatbarEnum,&latbar,P0Enum); 
+        this->AddInput(SealevelBarystaticIceLongbarEnum,&longbar,P0Enum); 
+        this->AddInput(SealevelBarystaticHydroLatbarEnum,&latbar,P0Enum); 
+        this->AddInput(SealevelBarystaticHydroLongbarEnum,&longbar,P0Enum); 
+        this->AddInput(SealevelBarystaticOceanLatbarEnum,&latbar,P0Enum); 
+        this->AddInput(SealevelBarystaticOceanLongbarEnum,&longbar,P0Enum); 
         #endif
         
@@ -6725 +6825 @@
 }
 /*}}}*/
-void       Tria::SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){ /*{{{*/
+void       Tria::SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){ /*{{{*/
 
         /*Declarations:{{{*/
@@ -6740 +6840 @@
         #ifdef _HAVE_RESTRICT_
         IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
-        IssmDouble* __restrict__ G_rigid_precomputed=NULL;
+        IssmDouble* __restrict__ G_gravi_precomputed=NULL;
         IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
         IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
@@ -6750 +6850 @@
         #else
         IssmDouble* G_viscoelastic_precomputed=NULL;
-        IssmDouble* G_rigid_precomputed=NULL;
+        IssmDouble* G_gravi_precomputed=NULL;
         IssmDouble* U_viscoelastic_precomputed=NULL;
         IssmDouble* H_viscoelastic_precomputed=NULL;
@@ -6759 +6859 @@
         #endif
         
-        /*elastic green function:*/
+        /*viscoelastic green function:*/
         int index;
-        int         M;
+        int M;
+        IssmDouble doubleindex,lincoef;
 
         /*Computational flags:*/
-        bool computerigid = false;
+        bool computeselfattraction = false;
         bool computeelastic = false;
         bool computeviscous = false;
         int  horiz;
+        int grd, grdmodel;
 
         bool istime=true;
@@ -6777 +6879 @@
         /*Recover parameters:{{{ */
         rho_earth=FindParam(MaterialsEarthDensityEnum);
-        this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
+        this->parameters->FindParam(&computeselfattraction,SolidearthSettingsSelfAttractionEnum);
         this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
         this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
@@ -6784 +6886 @@
         this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
         this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
+        this->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); 
+        this->parameters->FindParam(&grdmodel,GrdModelEnum);
         /*}}}*/
 
         /*early return:*/
-        if(!computerigid)return;
+        if (!grd || grdmodel!=ElasticEnum) return; //Love numbers won't be found in this case, return before loading them
+        if(!computeselfattraction)return;
 
         /*Recover precomputed green function kernels:{{{*/
-        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGRigidEnum)); _assert_(parameter);
-        parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
+        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGSelfAttractionEnum)); _assert_(parameter);
+        parameter->GetParameterValueByPointer((IssmDouble**)&G_gravi_precomputed,&M);
 
         if(computeelastic){
@@ -6826 +6931 @@
         }
         else{
-                nt=1; //in elastic, or if we run only rigid, we need only one step
+                nt=1; //in elastic, or if we run only selfattraction, we need only one step
         }
         Gsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
@@ -6884 +6989 @@
                                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
                                 alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0)));
-                                index=reCast<int,IssmDouble>( alpha/M_PI*reCast<IssmDouble,int>(M-1) );
-
+                                doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
+                                index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
+
+                                lincoef=doubleindex-index; //where between index and index+1 is alpha
+                                if (index==M-1){ //avoids out of bound case
+                                        index-=1;
+                                        lincoef=1;
+                                }
 
                                 for(int t=0;t<nt;t++){
@@ -6891 +7002 @@
 
                                         /*Rigid earth gravitational perturbation: */
-                                        Gsubel[l][timeindex]+=G_rigid_precomputed[index];
+                                        Gsubel[l][timeindex]+=G_gravi_precomputed[index+0]*(1.-lincoef) 
+                                                             +G_gravi_precomputed[index+1]*lincoef; //linear interpolation
 
                                         if(computeelastic){
-                                                Gsubel[l][timeindex]+=G_viscoelastic_precomputed[index*nt+t];
+                                                Gsubel[l][timeindex]+=G_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
+                                                                     +G_viscoelastic_precomputed[(index+1)*nt+t]*lincoef; //linear interpolation
                                         }
                                         Gsubel[l][timeindex]*=ratioe;
@@ -6900 +7013 @@
                                         /*Elastic components:*/
                                         if(computeelastic){
-                                                GUsubel[l][timeindex] =  ratioe * U_viscoelastic_precomputed[index*nt+t];
-
+                                                GUsubel[l][timeindex] =  ratioe *(U_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
+                                                                                 +U_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
                                                 if(horiz){
-                                                        GNsubel[l][timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*N_azim;
-                                                        GEsubel[l][timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*E_azim;
+                                                        GNsubel[l][timeindex] = N_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
+                                                                                              +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
+                                                        GEsubel[l][timeindex] = E_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
+                                                                                              +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
                                                 }
                                         }
@@ -6951 +7066 @@
         
         /*Inputs:*/
-        IssmDouble* stackRSL=NULL;
-        IssmDouble* stackU=NULL;
-        IssmDouble* stackN=NULL;
-        IssmDouble* stackE=NULL;
-        IssmDouble* stacktimes=NULL;
-        int         stacknumsteps;
-        int         stackindex=0; 
+        IssmDouble* viscousRSL=NULL;
+        IssmDouble* viscousU=NULL;
+        IssmDouble* viscousN=NULL;
+        IssmDouble* viscousE=NULL;
+        IssmDouble* viscoustimes=NULL;
+        int         viscousnumsteps;
+        int         viscousindex=0; 
         int         newindex=0; 
         int         dummy;
@@ -6970 +7085 @@
                 this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
                 
-                this->parameters->FindParam(&stacknumsteps,StackNumStepsEnum);
-                this->parameters->FindParam(&stacktimes,NULL,StackTimesEnum);
-                this->parameters->FindParam(&stackindex,StackIndexEnum);
+                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
+                this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
+                this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
                 this->parameters->FindParam(&currenttime,TimeEnum);
 
-                this->inputs->GetArrayPtr(StackRSLEnum,this->lid,&stackRSL,&dummy);
-                this->inputs->GetArrayPtr(StackUEnum,this->lid,&stackU,&dummy);
+                this->inputs->GetArrayPtr(SealevelchangeViscousRSLEnum,this->lid,&viscousRSL,&dummy);
+                this->inputs->GetArrayPtr(SealevelchangeViscousUEnum,this->lid,&viscousU,&dummy);
                 if(horiz){
-                        this->inputs->GetArrayPtr(StackNEnum,this->lid,&stackN,&dummy);
-                        this->inputs->GetArrayPtr(StackEEnum,this->lid,&stackE,&dummy);
+                        this->inputs->GetArrayPtr(SealevelchangeViscousNEnum,this->lid,&viscousN,&dummy);
+                        this->inputs->GetArrayPtr(SealevelchangeViscousEEnum,this->lid,&viscousE,&dummy);
                 }
 
@@ -6986 +7101 @@
                 int offset=1;
                 lincoeff=0;
-                newindex=stacknumsteps-2;
-
-                for(int t=stackindex;t<stacknumsteps;t++){
-                        if (stacktimes[t]>currenttime){
+                newindex=viscousnumsteps-2;
+
+                for(int t=viscousindex;t<viscousnumsteps;t++){
+                        if (viscoustimes[t]>currenttime){
                                 newindex=t-1;
-                                lincoeff=(currenttime-stacktimes[newindex])/(stacktimes[t]-stacktimes[newindex]);
+                                lincoeff=(currenttime-viscoustimes[newindex])/(viscoustimes[t]-viscoustimes[newindex]);
                                 foundtime=true;
                                 offset=0;
@@ -6999 +7114 @@
 
                 if(!foundtime) lincoeff=1;
-                stacktimes[newindex]=currenttime;
+                viscoustimes[newindex]=currenttime;
                 for(int i=0;i<NUMVERTICES;i++){
-                        stackRSL[i*stacknumsteps+newindex+offset]=(1-lincoeff)*stackRSL[i*stacknumsteps+newindex]+lincoeff*stackRSL[i*stacknumsteps+newindex+1];
-                        stackU[i*stacknumsteps+newindex+offset]=(1-lincoeff)*stackU[i*stacknumsteps+newindex]+lincoeff*stackU[i*stacknumsteps+newindex+1];
+                        viscousRSL[i*viscousnumsteps+newindex+offset]=(1-lincoeff)*viscousRSL[i*viscousnumsteps+newindex]+lincoeff*viscousRSL[i*viscousnumsteps+newindex+1];
+                        viscousU[i*viscousnumsteps+newindex+offset]=(1-lincoeff)*viscousU[i*viscousnumsteps+newindex]+lincoeff*viscousU[i*viscousnumsteps+newindex+1];
                         if(horiz){
-                                stackN[i*stacknumsteps+newindex+offset]=(1-lincoeff)*stackN[i*stacknumsteps+newindex]+lincoeff*stackN[i*stacknumsteps+newindex+1];
-                                stackE[i*stacknumsteps+newindex+offset]=(1-lincoeff)*stackE[i*stacknumsteps+newindex]+lincoeff*stackE[i*stacknumsteps+newindex+1];
-                        }
-                }
-                stackindex=newindex+offset;
+                                viscousN[i*viscousnumsteps+newindex+offset]=(1-lincoeff)*viscousN[i*viscousnumsteps+newindex]+lincoeff*viscousN[i*viscousnumsteps+newindex+1];
+                                viscousE[i*viscousnumsteps+newindex+offset]=(1-lincoeff)*viscousE[i*viscousnumsteps+newindex]+lincoeff*viscousE[i*viscousnumsteps+newindex+1];
+                        }
+                }
+                viscousindex=newindex+offset;
 
                 
-                this->parameters->SetParam(stackindex,StackIndexEnum);
-                this->parameters->SetParam(stacktimes,stacknumsteps,StackTimesEnum);
+                this->parameters->SetParam(viscousindex,SealevelchangeViscousIndexEnum);
+                this->parameters->SetParam(viscoustimes,viscousnumsteps,SealevelchangeViscousTimesEnum);
 
                 /*free allocations:*/
-                xDelete<IssmDouble>(stacktimes);
+                xDelete<IssmDouble>(viscoustimes);
         }
 
@@ -7099 +7214 @@
 }
 /*}}}*/
-void       Tria::SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/
+void       Tria::SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* polarmotionvector,SealevelGeometry* slgeom){ /*{{{*/
 
         /*sal green function:*/
         IssmDouble* G=NULL;
         IssmDouble* Gsub[SLGEOM_NUMLOADS];
-        IssmDouble SealevelGRD[NUMVERTICES]={0};
-        IssmDouble oceanaverage=0;
-        IssmDouble oceanarea=0;
-        IssmDouble rho_water;
         bool computefuture=false;
         
@@ -7120 +7231 @@
         IssmDouble Grotm3[3];
                 
-        this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
+        this->parameters->FindParam(&sal,SolidearthSettingsSelfAttractionEnum);
         this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
         this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
         this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
-        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
 
         if(sal){
@@ -7132 +7242 @@
                 this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
 
-                this->SealevelchangeGxL(sealevelpercpu, G, Gsub, loads, slgeom, nel,percpu=true,StackRSLEnum,computefuture=false);
+                this->SealevelchangeGxL(sealevelpercpu, G, Gsub, loads, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
         }
 
@@ -7142 +7252 @@
                 for(int i=0;i<NUMVERTICES;i++){
                         if(slgeom->lids[this->vertices[i]->lid]==this->lid){
-                                sealevelpercpu[this->vertices[i]->lid]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2];
+                                sealevelpercpu[this->vertices[i]->lid]+=Grotm1[i]*polarmotionvector[0]+Grotm2[i]*polarmotionvector[1]+Grotm3[i]*polarmotionvector[2];
                         }
                 }
@@ -7151 +7261 @@
 
         /*sal green function:*/
-        IssmDouble* G=NULL;
-        IssmDouble* Gsub[SLGEOM_NUMLOADS];
         IssmDouble oceanaverage=0;
         IssmDouble oceanarea=0;
@@ -7188 +7296 @@
         return;
 } /*}}}*/
-void       Tria::SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/
-
-        IssmDouble Sealevel[3]={0,0,0};
-        IssmDouble SealevelRSL[3]={0,0,0};
-        IssmDouble SealevelU[3]={0,0,0};
-        IssmDouble SealevelN[3]={0,0,0};
-        IssmDouble SealevelE[3]={0,0,0};
+void       Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* polarmotionvector,SealevelGeometry* slgeom){ /*{{{*/
+
+        IssmDouble SealevelGrd[3]={0,0,0};
+        IssmDouble RSLGrd[3]={0,0,0};
+        IssmDouble UGrd[3]={0,0,0};
+        IssmDouble NGrd[3]={0,0,0};
+        IssmDouble EGrd[3]={0,0,0};
         int nel,nbar;
         bool sal = false;
@@ -7224 +7332 @@
         bool elastic=false;
         bool percpu=false;
+        bool planethasocean=false;
 
         this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
-        this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
+        this->parameters->FindParam(&sal,SolidearthSettingsSelfAttractionEnum);
         this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
         this->parameters->FindParam(&elastic,SolidearthSettingsElasticEnum);
         this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
-
+        this->parameters->FindParam(&planethasocean,SolidearthSettingsGrdOceanEnum);
         
                 
@@ -7257 +7366 @@
                         }
                 }
-                this->SealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel,percpu=false,StackRSLEnum,computefuture=true);
+                this->SealevelchangeGxL(&RSLGrd[0],G, Gsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
 
                 if(elastic){
-                        this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel,percpu=false,StackUEnum,computefuture=true);
+                        this->SealevelchangeGxL(&UGrd[0],GU, GUsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
                         if(horiz ){
-                                this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel,percpu=false,StackNEnum,computefuture=true);
-                                this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel,percpu=false,StackEEnum,computefuture=true);
+                                this->SealevelchangeGxL(&NGrd[0],GN, GNsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
+                                this->SealevelchangeGxL(&EGrd[0],GE, GEsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
                         }
                 }
@@ -7275 +7384 @@
                 for(int i=0;i<NUMVERTICES;i++){
                         if(slgeom->lids[this->vertices[i]->lid]==this->lid){
-                                SealevelRSL[i]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2];
+                                RSLGrd[i]+=Grotm1[i]*polarmotionvector[0]+Grotm2[i]*polarmotionvector[1]+Grotm3[i]*polarmotionvector[2];
                         }
                 }
@@ -7286 +7395 @@
                         for(int i=0;i<NUMVERTICES;i++){
                                 if(slgeom->lids[this->vertices[i]->lid]==this->lid){
-                                        SealevelU[i]+=GUrotm1[i]*rotationvector[0]+GUrotm2[i]*rotationvector[1]+GUrotm3[i]*rotationvector[2];
+                                        UGrd[i]+=GUrotm1[i]*polarmotionvector[0]+GUrotm2[i]*polarmotionvector[1]+GUrotm3[i]*polarmotionvector[2];
                                 }
                         }
@@ -7299 +7408 @@
                                 for(int i=0;i<NUMVERTICES;i++){
                                         if(slgeom->lids[this->vertices[i]->lid]==this->lid){
-                                                SealevelN[i]+=GNrotm1[i]*rotationvector[0]+GNrotm2[i]*rotationvector[1]+GNrotm3[i]*rotationvector[2];
-                                                SealevelE[i]+=GErotm1[i]*rotationvector[0]+GErotm2[i]*rotationvector[1]+GErotm3[i]*rotationvector[2];
+                                                NGrd[i]+=GNrotm1[i]*polarmotionvector[0]+GNrotm2[i]*polarmotionvector[1]+GNrotm3[i]*polarmotionvector[2];
+                                                EGrd[i]+=GErotm1[i]*polarmotionvector[0]+GErotm2[i]*polarmotionvector[1]+GErotm3[i]*polarmotionvector[2];
                                         }
                                 }
@@ -7307 +7416 @@
         }
 
+        if (planethasocean){ //We must also output the RSL on vertices to compute the ocean mass conservation
+                for(int i=0;i<NUMVERTICES;i++){
+                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
+                                sealevelpercpu[this->vertices[i]->lid]=RSLGrd[i];
+                        }
+                }
+        }
+
         /*Create geoid: */
-        for(int i=0;i<NUMVERTICES;i++)Sealevel[i]=SealevelU[i]+SealevelRSL[i];
+        for(int i=0;i<NUMVERTICES;i++)SealevelGrd[i]=UGrd[i]+RSLGrd[i];
         
         /*Create inputs*/
-        this->AddInput(SealevelGRDEnum,Sealevel,P1Enum);
-        this->AddInput(BedGRDEnum,SealevelU,P1Enum);
+        this->AddInput(SealevelGRDEnum,SealevelGrd,P1Enum);
+        this->AddInput(BedGRDEnum,UGrd,P1Enum);
         if(horiz){
-                this->AddInput(BedNorthGRDEnum,SealevelN,P1Enum);
-                this->AddInput(BedEastGRDEnum,SealevelE,P1Enum);
+                this->AddInput(BedNorthGRDEnum,NGrd,P1Enum);
+                this->AddInput(BedEastGRDEnum,EGrd,P1Enum);
         }
 
 
 } /*}}}*/
-void       Tria::SealevelchangeGxL(IssmDouble* sealevelout, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu, int stackenum, bool computefuture) { /*{{{*/
-
-        IssmDouble* sealevel=NULL;
+void       Tria::SealevelchangeGxL(IssmDouble* grdfieldout, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/
+
+        IssmDouble* grdfield=NULL;
         int i,e,l,t,nbar;
-        bool viscous=false;
-        IssmDouble* stack=NULL;
+        bool computeviscous=false;
+        IssmDouble* viscousfield=NULL;
         int nt=1; //important
-        int stackindex=0; //important
-        int stacknumsteps=1; //important
-
-        this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
-        if(viscous){
-                this->parameters->FindParam(&stacknumsteps,StackNumStepsEnum);
-                if(computefuture) nt=stacknumsteps;
+        int viscousindex=0; //important
+        int viscousnumsteps=1; //important
+
+        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
+        if(computeviscous){
+                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
+                if(computefuture) nt=viscousnumsteps;
                 else nt=1;
 
                 //allocate
-                sealevel=xNewZeroInit<IssmDouble>(3*nt);
-        }
-        else sealevel=xNewZeroInit<IssmDouble>(3*nt);
+                grdfield=xNewZeroInit<IssmDouble>(3*nt);
+        }
+        else grdfield=xNewZeroInit<IssmDouble>(3*nt);
 
         if(loads->sealevelloads){
@@ -7347 +7464 @@
                         for (e=0;e<nel;e++){
                                 for(t=0;t<nt;t++){
-                                        int index=i*nel*stacknumsteps+e*stacknumsteps+t;
-                                        sealevel[i*nt+t]+=G[index]*(loads->sealevelloads[e]+loads->loads[e]);
+                                        int index=i*nel*viscousnumsteps+e*viscousnumsteps+t;
+                                        grdfield[i*nt+t]+=G[index]*(loads->sealevelloads[e]+loads->loads[e]);
                                 }
                         }
@@ -7355 +7472 @@
                                 for (e=0;e<nbar;e++){
                                         for(t=0;t<nt;t++){
-                                                int index=i*nbar*stacknumsteps+e*stacknumsteps+t;
-                                                sealevel[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
+                                                int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
+                                                grdfield[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
                                         }
                                 }
@@ -7362 +7479 @@
                                         for (e=0;e<nbar;e++){
                                                 for(t=0;t<nt;t++){
-                                                        int index=i*nbar*stacknumsteps+e*stacknumsteps+t;
-                                                        sealevel[i*nt+t]+=Gsub[l][index]*(loads->subsealevelloads[e]);
+                                                        int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
+                                                        grdfield[i*nt+t]+=Gsub[l][index]*(loads->subsealevelloads[e]);
                                                 }
                                         }
@@ -7376 +7493 @@
                         for (e=0;e<nel;e++){
                                 for(t=0;t<nt;t++){
-                                        int index=i*nel*stacknumsteps+e*stacknumsteps+t;
-                                        sealevel[i*nt+t]+=G[index]*(loads->loads[e]);
+                                        int index=i*nel*viscousnumsteps+e*viscousnumsteps+t;
+                                        grdfield[i*nt+t]+=G[index]*(loads->loads[e]);
                                 }
                         }
@@ -7384 +7501 @@
                                 for (e=0;e<nbar;e++){
                                         for(t=0;t<nt;t++){
-                                                int index=i*nbar*stacknumsteps+e*stacknumsteps+t;
-                                                sealevel[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
+                                                int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
+                                                grdfield[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
                                         }
                                 }
@@ -7392 +7509 @@
         }
         
-        if(viscous){
-                IssmDouble* sealevelinterp=NULL;
-                IssmDouble* stacktimes=NULL;
+        if(computeviscous){
+                IssmDouble* grdfieldinterp=NULL;
+                IssmDouble* viscoustimes=NULL;
                 IssmDouble  final_time;
                 IssmDouble  lincoeff;
                 IssmDouble  timeacc;
 
-                this->parameters->FindParam(&stackindex,StackIndexEnum);
-                this->parameters->FindParam(&stacktimes,NULL,StackTimesEnum);
+                this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
+                this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
                 this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
                 this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
-                this->inputs->GetArrayPtr(stackenum,this->lid,&stack,NULL);
+                this->inputs->GetArrayPtr(viscousenum,this->lid,&viscousfield,NULL);
                 if(computefuture){
-                        sealevelinterp=xNew<IssmDouble>(3*nt); 
-                        if(stacktimes[stackindex]<final_time){
-                                lincoeff=(stacktimes[stackindex+1]-stacktimes[stackindex])/timeacc;
-                                for(int t=stackindex;t<nt;t++){
-                                        if(t==stackindex){
+                        grdfieldinterp=xNew<IssmDouble>(3*nt); 
+                        if(viscoustimes[viscousindex]<final_time){
+                                lincoeff=(viscoustimes[viscousindex+1]-viscoustimes[viscousindex])/timeacc;
+                                for(int t=viscousindex;t<nt;t++){
+                                        if(t==viscousindex){
                                                 for(int i=0;i<NUMVERTICES;i++){
                                                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
-                                                        sealevelinterp[i*nt+t]=  sealevel[i*nt+0];
+                                                        grdfieldinterp[i*nt+t]=  grdfield[i*nt+0];
                                                 }
                                         }
@@ -7418 +7535 @@
                                                 for(int i=0;i<NUMVERTICES;i++){
                                                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
-                                                        sealevelinterp[i*nt+t]=  (1-lincoeff)*sealevel[i*nt+(t-stackindex-1)]+lincoeff*sealevel[i*nt+(t-stackindex)];
+                                                        grdfieldinterp[i*nt+t]=  (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]+lincoeff*grdfield[i*nt+(t-viscousindex)];
                                                 }
                                         }
@@ -7425 +7542 @@
                 }
 
-                /*update sealevel at present time using stack at present time: */
+                /*update grdfield at present time using viscous stack at present time: */
                 for(int i=0;i<NUMVERTICES;i++){
                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
-                        sealevel[i*nt+0]+=stack[i*stacknumsteps+stackindex]; 
-                }
-
-                if(computefuture){ /*update stack with future deformation from present load: */
-
-                        for(int t=nt-1;t>=stackindex;t--){
+                        grdfield[i*nt+0]+=viscousfield[i*viscousnumsteps+viscousindex]; 
+                }
+
+                if(computefuture){ /*update viscous stack with future deformation from present load: */
+
+                        for(int t=nt-1;t>=viscousindex;t--){
                                 for(int i=0;i<NUMVERTICES;i++){
                                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
-                                        stack[i*stacknumsteps+t]+=sealevelinterp[i*nt+t]-sealevelinterp[i*nt+stackindex]-stack[i*stacknumsteps+stackindex];
+                                        //offset viscousfield to remove all deformation that has already been added
+                                        viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*nt+t]-grdfieldinterp[i*nt+viscousindex]-viscousfield[i*viscousnumsteps+viscousindex];
                                 }
                         }
-                        /*Re-add stack now that we updated:*/
-                        this->inputs->SetArrayInput(stackenum,this->lid,stack,3*stacknumsteps);
+                        /*Re-add viscous stack now that we updated:*/
+                        this->inputs->SetArrayInput(viscousenum,this->lid,viscousfield,3*viscousnumsteps);
                 }
 
                 /*Free allocatoins:*/
-                xDelete<IssmDouble>(stacktimes);
+                xDelete<IssmDouble>(viscoustimes);
         }
 
@@ -7451 +7569 @@
                 for(i=0;i<NUMVERTICES;i++){
                         if(slgeom->lids[this->vertices[i]->lid]==this->lid){
-                                sealevelout[this->vertices[i]->lid]=sealevel[i*nt+0];
+                                grdfieldout[this->vertices[i]->lid]=grdfield[i*nt+0];
                         }
                 }
         }
         else{
-                for(i=0;i<NUMVERTICES;i++) sealevelout[i]=sealevel[i*nt+0];
+                for(i=0;i<NUMVERTICES;i++) grdfieldout[i]=grdfield[i*nt+0];
         }
 

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

@@ -170 +170 @@
                 void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
                 void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
-                void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom);
+                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom);
                 void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae);
                 void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
@@ -176 +176 @@
                 void       SealevelchangeConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom);
                 void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom);
-                void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom);
+                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom);
                 void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom);
                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads,  SealevelGeometry* slgeom);

issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp

@@ -24 +24 @@
 void TransferSealevel(FemModel* femmodel,int forcingenum);
 bool slcconvergence(Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
-IssmDouble  SealevelloadsOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble oceanarea);
-void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom);
+IssmDouble  SealevelloadsOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble totaloceanarea);
+void PolarMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom);
 void ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom);
 void ivins_deformation_core(FemModel* femmodel);
@@ -221 +221 @@
         BarystaticContributions* barycontrib=NULL;
         GenericParam<BarystaticContributions*>* barycontribparam=NULL;
-        IssmDouble rotationaxismotionvector[3]={0};
+        IssmDouble polarmotionvector[3]={0};
         
         GrdLoads*              loads=NULL;
         Vector<IssmDouble>*    oldsealevelloads=NULL;
         Vector<IssmDouble>*    oceanareas=NULL;
-        IssmDouble             oceanarea;
+        IssmDouble             totaloceanarea;
         Vector<IssmDouble>*    subelementoceanareas=NULL;
         IssmDouble             oceanaverage;
@@ -244 +244 @@
         int  grd=0;
         int  grdmodel; 
-        int  computesealevel=0;
+        int  sealevelloading=0;
         bool viscous=false;
+        bool rotation=false;
+        bool planethasocean=false;
         IssmDouble*           sealevelpercpu=NULL;
 
@@ -255 +257 @@
         /*retrieve parameters:{{{*/
         femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); 
+        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
         femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
         femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
-        femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum);
+        femmodel->parameters->FindParam(&step,StepEnum);
+        femmodel->parameters->FindParam(&time,TimeEnum);
+        femmodel->parameters->FindParam(&sealevelloading,SolidearthSettingsSealevelLoadingEnum);
         femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum);
-        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
         femmodel->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
+        femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
+        femmodel->parameters->FindParam(&planethasocean,SolidearthSettingsGrdOceanEnum);
         /*}}}*/
 
@@ -314 +320 @@
         loads->BroadcastLoads();
 
-        //compute rotation axis motion:
-        RotationAxisMotion(&rotationaxismotionvector[0],femmodel,loads,slgeom);
-
-        /*skip computation of sea level if requested, which means sea level loads should be zeroed */
-        if(!computesealevel){
+        //compute polar motion:
+        PolarMotion(&polarmotionvector[0],femmodel,loads,slgeom);
+
+        /*skip computation of sea level equation if requested, which means sea level loads should be zeroed */
+        if(!sealevelloading){
                 loads->sealevelloads=xNewZeroInit<IssmDouble>(nel);
                 loads->subsealevelloads=xNewZeroInit<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
+
                 goto deformation;
         }
@@ -332 +339 @@
                 for(Object* & object : femmodel->elements->objects){
                         Element* element = xDynamicCast<Element*>(object);
-                        element->SealevelchangeConvolution(sealevelpercpu, loads , rotationaxismotionvector,slgeom);
+                        element->SealevelchangeConvolution(sealevelpercpu, loads , polarmotionvector,slgeom);
                 }
                 
@@ -347 +354 @@
                         oceanareas->Assemble(); 
                         subelementoceanareas->Assemble();
-                        oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
-                        if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
+                        oceanareas->Sum(&totaloceanarea); _assert_(totaloceanarea>0.);
+                        if(scaleoceanarea) totaloceanarea=3.619e+14; // use true ocean area, m^2
                 }
         
                 //Conserve ocean mass: 
-                oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, oceanarea);
-                ConserveOceanMass(femmodel,loads,barycontrib->Total()/oceanarea - oceanaverage,slgeom);
+                oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, totaloceanarea);
+                ConserveOceanMass(femmodel,loads,barycontrib->Total()/totaloceanarea - oceanaverage,slgeom);
 
                 //broadcast sea level loads 
                 loads->BroadcastSealevelLoads();
 
-                //compute rotation axis motion:
-                RotationAxisMotion(&rotationaxismotionvector[0],femmodel,loads, slgeom);
+                //compute polar motion:
+                PolarMotion(&polarmotionvector[0],femmodel,loads, slgeom);
 
                 //convergence?
@@ -376 +383 @@
         for(Object* & object : femmodel->elements->objects){
                 Element* element = xDynamicCast<Element*>(object);
-                element->SealevelchangeDeformationConvolution(loads, rotationaxismotionvector,slgeom);
+                element->SealevelchangeDeformationConvolution(sealevelpercpu, loads, polarmotionvector,slgeom);
         }
 
         if(VerboseSolution()) _printf0_("         updating GRD fields\n");
 
-        /*Update bedrock motion and geoid:*/
-        if(computesealevel){
-                femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/oceanarea- oceanaverage/rho_water); //given that we converged, no need to recompute ocean average
+        if(planethasocean){
+                if (!sealevelloading){ //we haven't done so before, so we need to compute the ocean average and area
+                        loads->sealevelloads=NULL; //needed to trigger the calculation of areas
+                        /*retrieve sea level average  and ocean area:*/
+                        for(Object* & object : femmodel->elements->objects){
+                                Element* element = xDynamicCast<Element*>(object);
+                                element->SealevelchangeOceanAverage(loads, oceanareas, subelementoceanareas, sealevelpercpu, slgeom);
+                        }
+                        loads->sealevelloads=xNewZeroInit<IssmDouble>(nel);
+                        loads->AssembleSealevelLoads();
+                        /*compute ocean areas:*/
+                        oceanareas->Assemble(); 
+                        subelementoceanareas->Assemble();
+                        oceanareas->Sum(&totaloceanarea); _assert_(totaloceanarea>0.);
+                        if(scaleoceanarea) totaloceanarea=3.619e+14; // use true ocean area, m^2
+                                
+                        //Conserve ocean mass
+                        //Note that here we create sea-level loads but they will not generate GRD as we have already run all the convolutions
+                        oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, totaloceanarea);
+                        ConserveOceanMass(femmodel,loads,barycontrib->Total()/totaloceanarea - oceanaverage,slgeom);
+                }
+
+                femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/totaloceanarea- oceanaverage/rho_water);
 
                 //cumulate barystatic contributions and save to results: 
                 barycontrib->Cumulate(femmodel->parameters);
-                barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
+                barycontrib->Save(femmodel->results,femmodel->parameters,totaloceanarea);
                 barycontrib->Reset();
         }
 
+        if (rotation) {
+                femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,polarmotionvector[0],step,time));
+                femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,polarmotionvector[1],step,time));
+                femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorZZEnum,polarmotionvector[2],step,time));
+        }
+
+        /*Update bedrock motion and geoid:*/
         femmodel->inputs->AXPY(1,SealevelGRDEnum,SealevelEnum);
         femmodel->inputs->AXPY(1,BedGRDEnum,BedEnum);
@@ -409 +443 @@
         /*parameters: */
         int  step;
-        int computesealevel=0;
         bool isocean=false;
         IssmDouble time;
@@ -417 +450 @@
         IssmDouble cumgmslc=0;
         IssmDouble gmtslc=0;
-
-        /*early return if we are not computing sea level, but rather deformation: */
-        femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum);
-        if (!computesealevel)return;
 
         /*early return if we have no ocean transport:*/
@@ -610 +639 @@
         femmodel->parameters->FindParam(&areae,&nel,AreaeEnum);
         
-        /*initialize SealevelMasks structure: */
+        /*initialize SealevelloadMasks structure: */
         slgeom=new SealevelGeometry(femmodel->elements->Size(),femmodel->vertices->Size());
         
@@ -639 +668 @@
         for(Object* & object : femmodel->elements->objects){
                 Element*   element=xDynamicCast<Element*>(object);
-                element->SealevelchangeGeometryFractionKernel(slgeom);
+                element->SealevelchangeGeometrySubElementKernel(slgeom);
         }
 
@@ -699 +728 @@
 
 } /*}}}*/
-IssmDouble  SealevelloadsOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* suboceanareas, IssmDouble oceanarea){ /*{{{*/
+IssmDouble  SealevelloadsOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* suboceanareas, IssmDouble totaloceanarea){ /*{{{*/
 
         IssmDouble sealevelloadsaverage;        
@@ -715 +744 @@
         delete vsubsealevelloadsvolume; 
         
-        return (sealevelloadsaverage+subsealevelloadsaverage)/oceanarea;
+        return (sealevelloadsaverage+subsealevelloadsaverage)/totaloceanarea;
 } /*}}}*/
-void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom){ /*{{{*/
+void PolarMotion(IssmDouble* polarmotionvector, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom){ /*{{{*/
 
         IssmDouble  moi_list[3]={0,0,0};
@@ -746 +775 @@
 
                 element->SealevelchangeMomentOfInertiaCentroid(&moi_list[0],loads,slgeom);
-
                 element->SealevelchangeMomentOfInertiaSubElement(&moi_list_sub[0],loads, slgeom);
 
@@ -754 +782 @@
         }
 
+        for (int i=0;i<3;i++) moi_list[i]=0;
 
         ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
@@ -770 +799 @@
 
         /*Assign output pointers:*/
-        m[0]=m1;
-        m[1]=m2;
-        m[2]=m3;
+        polarmotionvector[0]=m1;
+        polarmotionvector[1]=m2;
+        polarmotionvector[2]=m3;
 } /*}}}*/
 void        ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -360 +360 @@
         SolidearthSettingsViscousEnum,
         SealevelchangeGeometryDoneEnum,
+        SealevelchangeViscousNumStepsEnum,
+        SealevelchangeViscousTimesEnum,
+        SealevelchangeViscousIndexEnum,
         RotationalEquatorialMoiEnum,
         TidalLoveHEnum,
@@ -370 +373 @@
         LoveTimeFreqEnum,
         LoveIsTimeEnum,
-        SealevelchangeGRigidEnum,
+        SealevelchangeGSelfAttractionEnum,
         SealevelchangeGViscoElasticEnum,
-        SolidearthSettingsComputesealevelchangeEnum,
+        SolidearthSettingsSealevelLoadingEnum,
         SolidearthSettingsGRDEnum,
         SolidearthSettingsRunFrequencyEnum,
@@ -379 +382 @@
         SolidearthSettingsHorizEnum,
         SolidearthSettingsMaxiterEnum,
+        SolidearthSettingsGrdOceanEnum,
         SolidearthSettingsOceanAreaScalingEnum,
         RotationalPolarMoiEnum,
         SolidearthSettingsReltolEnum,
         SealevelchangeRequestedOutputsEnum,
-        SolidearthSettingsRigidEnum,
+        SolidearthSettingsSelfAttractionEnum,
         SolidearthSettingsRotationEnum,
         SolidearthSettingsMaxSHCoeffEnum,
@@ -397 +401 @@
         SettingsSolverResidueThresholdEnum,
         SettingsWaitonlockEnum,
-        StackNumStepsEnum,
-        StackTimesEnum,
-        StackIndexEnum,
         SmbAIceEnum,
         SmbAIdxEnum,
@@ -832 +833 @@
         SealevelchangeGEsubelHydroEnum,
         SealevelchangeGNsubelHydroEnum,
+        SealevelchangeViscousRSLEnum,
+        SealevelchangeViscousUEnum,
+        SealevelchangeViscousNEnum,
+        SealevelchangeViscousEEnum,
         SedimentHeadEnum,
         SedimentHeadOldEnum,
@@ -953 +958 @@
         SolidearthExternalDisplacementUpRateEnum,
         SolidearthExternalGeoidRateEnum,
-        StackRSLEnum,
-        StackUEnum,
-        StackNEnum,
-        StackEEnum,
         StrainRateeffectiveEnum,
         StrainRateparallelEnum,
@@ -1430 +1431 @@
         SealevelInertiaTensorYZEnum,
         SealevelInertiaTensorZZEnum,
+        SealevelchangePolarMotionEnum,
         SealevelNmotionEnum,
         SealevelUmotionEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -378 +378 @@
                 case LoveTimeFreqEnum : return "LoveTimeFreq";
                 case LoveIsTimeEnum : return "LoveIsTime";
-                case SealevelchangeGRigidEnum : return "SealevelchangeGRigid";
+                case SealevelchangeGSelfAttractionEnum : return "SealevelchangeGSelfAttraction";
                 case SealevelchangeGViscoElasticEnum : return "SealevelchangeGViscoElastic";
-                case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange";
+                case SolidearthSettingsSealevelLoadingEnum : return "SolidearthSettingsSealevelLoading";
                 case SolidearthSettingsGRDEnum : return "SolidearthSettingsGRD";
                 case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency";
@@ -391 +391 @@
                 case SolidearthSettingsReltolEnum : return "SolidearthSettingsReltol";
                 case SealevelchangeRequestedOutputsEnum : return "SealevelchangeRequestedOutputs";
-                case SolidearthSettingsRigidEnum : return "SolidearthSettingsRigid";
+                case SolidearthSettingsSelfAttractionEnum : return "SolidearthSettingsSelfAttraction";
                 case SolidearthSettingsRotationEnum : return "SolidearthSettingsRotation";
                 case SealevelchangeRunCountEnum : return "SealevelchangeRunCount";
@@ -404 +404 @@
                 case SettingsSolverResidueThresholdEnum : return "SettingsSolverResidueThreshold";
                 case SettingsWaitonlockEnum : return "SettingsWaitonlock";
-                case StackNumStepsEnum : return "StackNumSteps";
-                case StackTimesEnum : return "StackTimes";
-                case StackIndexEnum : return "StackIndex";
+                case SealevelchangeViscousNumStepsEnum : return "SealevelchangeViscousNumSteps";
+                case SealevelchangeViscousTimesEnum : return "SealevelchangeViscousTimes";
+                case SealevelchangeViscousIndexEnum : return "SealevelchangeViscousIndex";
                 case SmbAIceEnum : return "SmbAIce";
                 case SmbAIdxEnum : return "SmbAIdx";
@@ -957 +957 @@
                 case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate";
                 case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate";
-                case StackRSLEnum : return "StackRSL";
-                case StackUEnum : return "StackU";
-                case StackNEnum : return "StackN";
-                case StackEEnum : return "StackE";
+                case SealevelchangeViscousRSLEnum : return "SealevelchangeViscousRSL";
+                case SealevelchangeViscousUEnum : return "SealevelchangeViscousU";
+                case SealevelchangeViscousNEnum : return "SealevelchangeViscousN";
+                case SealevelchangeViscousEEnum : return "SealevelchangeViscousE";
                 case StrainRateeffectiveEnum : return "StrainRateeffective";
                 case StrainRateparallelEnum : return "StrainRateparallel";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -387 +387 @@
    if(stage==4){
               if (strcmp(name,"LoveIsTime")==0) return LoveIsTimeEnum;
-              else if (strcmp(name,"SealevelchangeGRigid")==0) return SealevelchangeGRigidEnum;
+              else if (strcmp(name,"SealevelchangeGSelfAttraction")==0) return SealevelchangeGSelfAttractionEnum;
               else if (strcmp(name,"SealevelchangeGViscoElastic")==0) return SealevelchangeGViscoElasticEnum;
-              else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum;
+              else if (strcmp(name,"SolidearthSettingsSealevelLoading")==0) return SolidearthSettingsSealevelLoadingEnum;
               else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum;
               else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
@@ -400 +400 @@
               else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum;
               else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum;
-              else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum;
+              else if (strcmp(name,"SolidearthSettingsSelfAttraction")==0) return SolidearthSettingsSelfAttractionEnum;
               else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum;
               else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum;
@@ -413 +413 @@
               else if (strcmp(name,"SettingsSolverResidueThreshold")==0) return SettingsSolverResidueThresholdEnum;
               else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
-              else if (strcmp(name,"StackNumSteps")==0) return StackNumStepsEnum;
-              else if (strcmp(name,"StackTimes")==0) return StackTimesEnum;
-              else if (strcmp(name,"StackIndex")==0) return StackIndexEnum;
+              else if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum;
+              else if (strcmp(name,"SealevelchangeViscousTimes")==0) return SealevelchangeViscousTimesEnum;
+              else if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum;
               else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum;
               else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum;
@@ -978 +978 @@
               else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum;
               else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
-              else if (strcmp(name,"StackRSL")==0) return StackRSLEnum;
-              else if (strcmp(name,"StackU")==0) return StackUEnum;
-              else if (strcmp(name,"StackN")==0) return StackNEnum;
-              else if (strcmp(name,"StackE")==0) return StackEEnum;
+              else if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum;
+              else if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum;
+              else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum;
+              else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum;
               else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
               else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;

issm/trunk-jpl/src/m/classes/solidearthsettings.m

@@ -9 +9 @@
                 abstol                 = 0;
                 maxiter                = 0;
-                rigid                  = 0;
-                elastic                = 0;
-                viscous                = 0;
-                rotation               = 0;
+                selfattraction         = 1;
+                elastic                = 1;
+                viscous                = 1;
+                rotation               = 1;
+                grdocean               = 1;
                 ocean_area_scaling     = 0;
                 runfrequency           = 1; %how many time steps we skip before we run grd_core
-                computesealevelchange  = 1; %will sea-level be coputed? 
+                sealevelloading        = 1; %will sea-level loads be computed? 
                 isgrd                  = 0; %will GRD patterns be computed? 
                 compute_bp_grd         = 0; %will GRD patterns for bottom pressures be computed? 
@@ -21 +22 @@
                 timeacc                = 0; %time step accurary required to compute Green tables
                 horiz                  = 0; %compute horizontal deformation
-                grdmodel               = 0; %grd model (0 by default, 1 for elastic, 2 for Ivins)
+                grdmodel               = 0; %grd model (0 by default, 1 for (visco)-elastic, 2 for Ivins)
                 cross_section_shape    = 0; %cross section only used when grd model is Ivins
         end
+        methods (Static)
+                function self = loadobj(self) % {{{
+                        % This function is directly called by matlab when a model object is
+                        % loaded. If the input is a struct it is an old version of this class and
+                        % old fields must be recovered (make sure they are in the deprecated
+                        % model properties)
+
+                        if isstruct(self)
+                                % 2021, Jun 4
+                                if isfield(self,'rigid')
+                                        self.selfattraction = self.rigid;
+                                end
+                                if isfield(self,'computesealevelchange')
+                                        self.sealevelloading = self.computesealevelchange;
+                                end
+                                self = structtoobj(solidearthsettings(),self);
+
+                        end
+                end% }}}
+        end
         methods
+
                 function self = solidearthsettings(varargin) % {{{
                         switch nargin
@@ -43 +65 @@
 
                 %computational flags: 
-                self.rigid=1;
+                self.selfattraction=1;
                 self.elastic=1;
                 self.viscous=1;
                 self.rotation=1;
+                self.grdocean=1;
                 self.ocean_area_scaling=0;
                 self.compute_bp_grd=0;
                 self.isgrd=0;
-                self.computesealevelchange=1;
+                self.sealevelloading=1;
 
                 %numerical discretization accuracy
@@ -83 +106 @@
                         md = checkfield(md,'fieldname','solidearth.settings.grdmodel','>=',0,'<=',2);
                         md = checkfield(md,'fieldname','solidearth.settings.cross_section_shape','numel',[1],'values',[1,2]);
-
-                        %checks on computational flags
-                        if self.elastic==1 & self.rigid==0,
-                                error('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set');
+                        if self.elastic==1 & self.selfattraction==0,
+                                error('solidearthsettings checkconsistency error message: need selfattraction on if elastic flag is set');
                         end
                         if self.viscous==1 & self.elastic==0,
@@ -103 +124 @@
                                         end
                                 end
+                                if self.sealevelloading==1 & self.grdocean==0
+                                        error('solidearthsettings checkconsistency error message: need grdocean on if sealevelloading flag is set');
+                                end
                         end
 
@@ -116 +140 @@
                         fielddisplay(self,'abstol','sea level change absolute convergence criterion, NaN: not applied');
                         fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
+                        fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module [default: 1]'); 
                         fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); 
-                        fielddisplay(self,'computesealevelchange','compute sealevel change (default 1)');
+                        fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default 1)');
                         fielddisplay(self,'isgrd','compute GRD patterns (default 1)');
                         fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default 1)');
                         fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)');
-                        fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
-                        fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
-                        fielddisplay(self,'viscous','viscous earth graviational potential perturbation');
-                        fielddisplay(self,'rotation','earth rotational potential perturbation');
+                        fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field');
+                        fielddisplay(self,'elastic','enables elastic deformation from surface loading');
+                        fielddisplay(self,'viscous','enables viscous deformation from surface loading');
+                        fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields');
                         fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
                         fielddisplay(self,'timeacc','time accuracy (default 1 yr) for numerical discretization of the Green''s functions');
-                        fielddisplay(self,'grdmodel','type of deformation model, 1 for elastic, 2 for visco-elastic from Ivins');
+                        fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)');
                         fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
                 end % }}}
@@ -134 +159 @@
                         WriteData(fid,prefix,'object',self,'fieldname','abstol','name','md.solidearth.settings.abstol','format','Double');
                         WriteData(fid,prefix,'object',self,'fieldname','maxiter','name','md.solidearth.settings.maxiter','format','Integer');
-                        WriteData(fid,prefix,'object',self,'fieldname','rigid','name','md.solidearth.settings.rigid','format','Boolean');
+                        WriteData(fid,prefix,'object',self,'fieldname','selfattraction','name','md.solidearth.settings.selfattraction','format','Boolean');
                         WriteData(fid,prefix,'object',self,'fieldname','elastic','name','md.solidearth.settings.elastic','format','Boolean');
                         WriteData(fid,prefix,'object',self,'fieldname','viscous','name','md.solidearth.settings.viscous','format','Boolean');
                         WriteData(fid,prefix,'object',self,'fieldname','rotation','name','md.solidearth.settings.rotation','format','Boolean');
+                        WriteData(fid,prefix,'object',self,'fieldname','grdocean','name','md.solidearth.settings.grdocean','format','Boolean');
                         WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','name','md.solidearth.settings.ocean_area_scaling','format','Boolean');
                         WriteData(fid,prefix,'object',self,'fieldname','runfrequency','name','md.solidearth.settings.runfrequency','format','Integer');
@@ -143 +169 @@
                         WriteData(fid,prefix,'object',self,'fieldname','timeacc','name','md.solidearth.settings.timeacc','format','Double','scale',md.constants.yts);
                         WriteData(fid,prefix,'object',self,'fieldname','horiz','name','md.solidearth.settings.horiz','format','Integer');
-                        WriteData(fid,prefix,'object',self,'fieldname','computesealevelchange','name','md.solidearth.settings.computesealevelchange','format','Integer');
+                        WriteData(fid,prefix,'object',self,'fieldname','sealevelloading','name','md.solidearth.settings.sealevelloading','format','Integer');
                         WriteData(fid,prefix,'object',self,'fieldname','isgrd','name','md.solidearth.settings.isgrd','format','Integer');
                         WriteData(fid,prefix,'object',self,'fieldname','compute_bp_grd','name','md.solidearth.settings.compute_bp_grd','format','Integer');
@@ -154 +180 @@
                         writejsdouble(fid,[modelname '.solidearth.settings.reltol'],self.reltol);
                         writejsdouble(fid,[modelname '.solidearth.settings.abstol'],self.abstol);
-                        writejsdouble(fid,[modelname '.solidearth.settings.rigid'],self.rigid);
+                        writejsdouble(fid,[modelname '.solidearth.settings.selfattraction'],self.selfattraction);
                         writejsdouble(fid,[modelname '.solidearth.settings.elastic'],self.elastic);
                         writejsdouble(fid,[modelname '.solidearth.settings.viscous'],self.viscous);
                         writejsdouble(fid,[modelname '.solidearth.settings.rotation'],self.rotation);
+                        writejsdouble(fid,[modelname '.solidearth.settings.grdocean'],self.rotation);
                         writejsdouble(fid,[modelname '.solidearth.settings.ocean_area_scaling'],self.ocean_area_scaling);
                         writejsdouble(fid,[modelname '.solidearth.settings.run_frequency'],self.run_frequency);

issm/trunk-jpl/test/NightlyRun/test2001.m

@@ -17 +17 @@
 md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
 md.solidearth.settings.cross_section_shape=1;    % for square-edged x-section 
-md.solidearth.settings.computesealevelchange=0;  %do not compute sea level, only deformation
+md.solidearth.settings.grdocean=0;  %do not compute sea level, only deformation
+md.solidearth.settings.sealevelloading=0;  %do not compute sea level, only deformation
 md.solidearth.requested_outputs={'Sealevel','BedGRD'};
 

issm/trunk-jpl/test/NightlyRun/test2002.m

@@ -3 +3 @@
 %mesh earth:
 md=model;
-md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500.); %700 km resolution mesh
+load ../Data/SlcTestMesh.mat;
+md.mesh=SlcMesh; %700 km resolution mesh
 
-%Geometry for the bed, arbitrary thickness of 1000: 
+%Geometry for the bed, arbitrary thickness of 100: 
 md.geometry.bed=zeros(md.mesh.numberofvertices,1);
 md.geometry.base=md.geometry.bed;
@@ -70 +71 @@
 md.solidearth.settings.reltol=NaN;
 md.solidearth.settings.abstol=1e-3;
-md.solidearth.settings.computesealevelchange=1;
+md.solidearth.settings.sealevelloading=1;
 md.solidearth.settings.isgrd=1;
 md.solidearth.settings.ocean_area_scaling=0;
@@ -87 +88 @@
 
 %eustatic run:
-md.solidearth.settings.rigid=0;
+md.solidearth.settings.selfattraction=0;
 md.solidearth.settings.elastic=0;
 md.solidearth.settings.rotation=0;
@@ -96 +97 @@
 Beustatic=md.results.TransientSolution.Bed;
 
-%eustatic + rigid run:
-md.solidearth.settings.rigid=1;
+%eustatic + selfattraction run:
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=0;
 md.solidearth.settings.rotation=0;
 md.solidearth.settings.viscous=0;
 md=solve(md,'tr');
-Srigid=md.results.TransientSolution.Sealevel;
-Brigid=md.results.TransientSolution.Bed;
+Sselfattraction=md.results.TransientSolution.Sealevel;
+Bselfattraction=md.results.TransientSolution.Bed;
 
-%eustatic + rigid + elastic run:
-md.solidearth.settings.rigid=1;
+%eustatic + selfattraction + elastic run:
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
 md.solidearth.settings.rotation=0;
@@ -114 +115 @@
 Belastic=md.results.TransientSolution.Bed;
 
-%eustatic + rigid + elastic + rotation run:
-md.solidearth.settings.rigid=1;
+%eustatic + selfattraction + elastic + rotation run:
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
 md.solidearth.settings.rotation=1;
@@ -124 +125 @@
 
 %Fields and tolerances to track changes
-field_names={'Seustatic','Srigid','Selastic','Srotation','Beustatic','Brigid','Belastic','Brotation'};
+field_names={'Seustatic','Sselfattraction','Selastic','Srotation','Beustatic','Bselfattraction','Belastic','Brotation'};
 field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13};
-field_values={Seustatic,Srigid,Selastic,Srotation,Beustatic,Brigid,Belastic,Brotation};
+field_values={Seustatic,Sselfattraction,Selastic,Srotation,Beustatic,Bselfattraction,Belastic,Brotation};

issm/trunk-jpl/test/NightlyRun/test2003.m

@@ -3 +3 @@
 %mesh earth:
 md=model;
-md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1000.); %1000 km resolution mesh
+load ../Data/SlcTestMesh.mat;
+md.mesh=SlcMesh; %700 km resolution mesh
 
 %Geometry for the bed, arbitrary thickness of 1000: 
@@ -62 +63 @@
 md.solidearth.settings.reltol=NaN;
 md.solidearth.settings.abstol=1e-3;
-md.solidearth.settings.computesealevelchange=0;
+md.solidearth.settings.sealevelloading=0;
+md.solidearth.settings.grdocean=0;
 md.solidearth.settings.isgrd=1;
 md.solidearth.settings.ocean_area_scaling=0;
@@ -68 +70 @@
 md.solidearth.settings.horiz=1;
 md.solidearth.requested_outputs={'Sealevel','Bed', 'BedEast', 'BedNorth'};
-
 
 %Physics: 
@@ -81 +82 @@
 md.timestepping.final_time=1;
 
-%eustatic + rigid + elastic run:
-md.solidearth.settings.rigid=1;
+%eustatic + selfattraction + elastic run:
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
 md.solidearth.settings.rotation=0;
@@ -94 +95 @@
 BNnoRotation=md.results.TransientSolution.BedNorth;
 
-%eustatic + rigid + elastic + rotation run:
-md.solidearth.settings.rigid=1;
+%eustatic + selfattraction + elastic + rotation run:
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
 md.solidearth.settings.rotation=1;
@@ -108 +109 @@
 
 %Fields and tolerances to track changes
-field_names     ={'noRotation','Rotation'};
+field_names     ={'Sealevel', 'Uplift', 'NorthDisplacement', 'EastDisplacement'};
 field_tolerances={1e-13,1e-13,1e-13,1e-13};
 field_values={SRotation-SnoRotation,BURotation-BUnoRotation,BNRotation-BNnoRotation,BERotation-BEnoRotation};

issm/trunk-jpl/test/NightlyRun/test2004.m

@@ -395 +395 @@
 md.solidearth.settings.reltol=NaN;
 md.solidearth.settings.abstol=1e-3;
-md.solidearth.settings.computesealevelchange=1;
+md.solidearth.settings.sealevelloading=1;
 md.solidearth.settings.isgrd=1;
 md.solidearth.settings.ocean_area_scaling=0;
@@ -423 +423 @@
 
 %eustatic run:
-md.solidearth.settings.rigid=0;
+md.solidearth.settings.selfattraction=0;
 md.solidearth.settings.elastic=0;
 md.solidearth.settings.rotation=0;
@@ -432 +432 @@
 Seustatic=md.results.TransientSolution.Sealevel;
 
-%eustatic + rigid run:
-md.solidearth.settings.rigid=1;
+%eustatic + selfattraction run:
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=0;
 md.solidearth.settings.rotation=0;
 md=solve(md,'Transient');
-Srigid=md.results.TransientSolution.Sealevel;
-
-%eustatic + rigid + elastic run:
-md.solidearth.settings.rigid=1;
+Sselfattraction=md.results.TransientSolution.Sealevel;
+
+%eustatic + selfattraction + elastic run:
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
 md.solidearth.settings.rotation=0;
@@ -446 +446 @@
 Selastic=md.results.TransientSolution.Sealevel;
 
-%eustatic + rigid + elastic + rotation run:
-md.solidearth.settings.rigid=1;
+%eustatic + selfattraction + elastic + rotation run:
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
 md.solidearth.settings.rotation=1;
@@ -457 +457 @@
 field_names     ={'Eustatic','Rigid','Elastic','Rotation'};
 field_tolerances={1e-13,1e-13,1e-13,1e-13};
-field_values={Seustatic,Srigid,Selastic,Srotation};
+field_values={Seustatic,Sselfattraction,Selastic,Srotation};

issm/trunk-jpl/test/NightlyRun/test2005.m

@@ -1 @@
-%Test Name: EarthSlc
+%Test Name: Earthslc_time
 
 %mesh earth:
 md=model;
-md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh
+load ../Data/SlcTestMesh.mat;
+md.mesh=SlcMesh; %700 km resolution mesh
 
-%parameterize solidearth solution:
+%Geometry for the bed, arbitrary thickness of 100: 
+md.geometry.bed=zeros(md.mesh.numberofvertices,1);
+md.geometry.base=md.geometry.bed;
+md.geometry.thickness=100*ones(md.mesh.numberofvertices,1);
+md.geometry.surface=md.geometry.bed+md.geometry.thickness;
+
 %solidearth loading:  {{{
-md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
-md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
-md.dsl.global_average_thermosteric_sea_level_change=[1;0];
-md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
-md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
+md.masstransport.spcthickness=[md.geometry.thickness;0];
+md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
 %antarctica
-late=sum(md.mesh.lat(md.mesh.elements),2)/3;
-longe=sum(md.mesh.long(md.mesh.elements),2)/3;
+xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
+ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3;
+ze=md.mesh.z(md.mesh.elements)*[1;1;1]/3;
+re=sqrt(xe.^2+ye.^2+ze.^2);
+
+late=asind(ze./re);
+longe=atan2d(ye,xe);
 pos=find(late < -80);
-md.solidearth.surfaceload.icethicknesschange(pos)=-100;
+md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
+posant=pos;
 %greenland
-pos=find(late>70 & late<80 & longe>-60 & longe<-30);
-md.solidearth.surfaceload.icethicknesschange(pos)=-100;
+pos=find(late>60 & late<90 & longe>-75 & longe<-15);
+md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
+posgre=pos;
 
 %elastic loading from love numbers:
-md.solidearth.lovenumbers=lovenumbers('maxdeg',100);
+md.solidearth.lovenumbers=lovenumbers('maxdeg',1000);
 
 %}}}
 %mask:  {{{
 mask=gmtmask(md.mesh.lat,md.mesh.long);
+oceanmask=-ones(md.mesh.numberofvertices,1);
+pos=find(mask==0); oceanmask(pos)=1;
+
 icemask=ones(md.mesh.numberofvertices,1);
-pos=find(mask==0);
-icemask(pos)=-1;
-pos=find(sum(mask(md.mesh.elements),2)<3);
-icemask(md.mesh.elements(pos,:))=-1;
+icemask(md.mesh.elements(posant))=-1;
+icemask(md.mesh.elements(posgre))=-1;
+
 md.mask.ice_levelset=icemask;
-md.mask.ocean_levelset=-icemask;
-
-%make sure that the elements that have loads are fully grounded:
-pos=find(md.solidearth.surfaceload.icethicknesschange);
-md.mask.ocean_levelset(md.mesh.elements(pos,:))=1;
-
-%make sure wherever there is an ice load, that the mask is set to ice:
-%pos=find(md.solidearth.surfaceload.icethicknesschange); % TODO: Do we need to do this twice?
-md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
+md.mask.ocean_levelset=oceanmask;
 % }}}
 
-md.solidearth.settings.ocean_area_scaling=0;
+%time stepping: 
+md.timestepping.start_time=0;
+md.timestepping.time_step=1;
+md.timestepping.final_time=10;
 
-%Geometry for the bed, arbitrary: 
-md.geometry.bed=-ones(md.mesh.numberofvertices,1);
+%masstransport:
+md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
+md.initialization.str=0;
 
 %Materials: 
@@ -56 +68 @@
 
 %Solution parameters
+md.cluster.np=3;
 md.solidearth.settings.reltol=NaN;
 md.solidearth.settings.abstol=1e-3;
-md.solidearth.settings.computesealevelchange=1;
+md.solidearth.settings.sealevelloading=1;
+md.solidearth.settings.isgrd=1;
+md.solidearth.settings.ocean_area_scaling=0;
+md.solidearth.settings.grdmodel=1;
 
-% max number of iteration reverted back to 10 (i.e., the original default value)
-md.solidearth.settings.maxiter=10;
-
-%eustatic + rigid + elastic + rotation run:
-md.solidearth.settings.rigid=1;
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
 md.solidearth.settings.rotation=1;
+md.solidearth.settings.viscous=0;
 
-%transient settings: 
-md.timestepping.start_time=0;
-md.timestepping.final_time=10;
-md.timestepping.time_step=1;
-md.transient.isslc=1;
-md.transient.issmb=0;
-md.transient.ismasstransport=0;
+%Physics: 
+md.transient.issmb=0; 
 md.transient.isstressbalance=0;
 md.transient.isthermal=0;
-dh=md.solidearth.surfaceload.icethicknesschange;
-deltathickness=zeros(md.mesh.numberofelements+1,10);
-for i=1:10,
-        deltathickness(1:end-1,i)=dh*i;
+md.transient.ismasstransport=1;
+md.transient.isslc=1;
+md.solidearth.requested_outputs={'Sealevel'};
+
+dh=md.masstransport.spcthickness;
+deltathickness=zeros(md.mesh.numberofvertices+1,10);
+for i=0:10,
+        deltathickness(1:end-1,i+1)=md.geometry.thickness+dh(1:end-1)*i;
 end
-deltathickness(end,:)=0:1:9;
-md.solidearth.surfaceload.icethicknesschange=deltathickness;
+deltathickness(end,:)=0:1:10;
+md.masstransport.spcthickness=deltathickness;
 %hack: 
 md.geometry.surface=zeros(md.mesh.numberofvertices,1);
@@ -89 +101 @@
 md.geometry.base=-ones(md.mesh.numberofvertices,1);
 md.geometry.bed=md.geometry.base;
+
 
 %run transient solution: 

issm/trunk-jpl/test/NightlyRun/test2006.m

@@ -3 +3 @@
 %mesh earth:
 md=model;
-md.cluster=generic('name',oshostname(),'np',5);
-md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh
-
-%parameterize solidearth solution:
+load ../Data/SlcTestMesh.mat;
+md.mesh=SlcMesh; %700 km resolution mesh
+
+%Geometry for the bed, arbitrary thickness of 100: 
+md.geometry.bed=zeros(md.mesh.numberofvertices,1);
+md.geometry.base=md.geometry.bed;
+md.geometry.thickness=100*ones(md.mesh.numberofvertices,1);
+md.geometry.surface=md.geometry.bed+md.geometry.thickness;
+
 %solidearth loading:  {{{
-md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
-md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
-md.dsl.global_average_thermosteric_sea_level_change=[1;0];
-md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
-md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
+md.masstransport.spcthickness=[md.geometry.thickness;0];
+md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
 %antarctica
-late=sum(md.mesh.lat(md.mesh.elements),2)/3;
-longe=sum(md.mesh.long(md.mesh.elements),2)/3;
+xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
+ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3;
+ze=md.mesh.z(md.mesh.elements)*[1;1;1]/3;
+re=sqrt(xe.^2+ye.^2+ze.^2);
+
+late=asind(ze./re);
+longe=atan2d(ye,xe);
 pos=find(late < -80);
-md.solidearth.surfaceload.icethicknesschange(pos)=-100;
+md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
+posant=pos;
 %greenland
-pos=find(late>70 & late<80 & longe>-60 & longe<-30);
-md.solidearth.surfaceload.icethicknesschange(pos)=-100;
+pos=find(late>60 & late<90 & longe>-75 & longe<-15);
+md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
+posgre=pos;
 
 %elastic loading from love numbers:
@@ -28 +37 @@
 %mask:  {{{
 mask=gmtmask(md.mesh.lat,md.mesh.long);
+oceanmask=-ones(md.mesh.numberofvertices,1);
+pos=find(mask==0); oceanmask(pos)=1;
+
 icemask=ones(md.mesh.numberofvertices,1);
-pos=find(mask==0);
-icemask(pos)=-1;
-pos=find(sum(mask(md.mesh.elements),2)<3);
-icemask(md.mesh.elements(pos,:))=-1;
+icemask(md.mesh.elements(posant))=-1;
+icemask(md.mesh.elements(posgre))=-1;
+
 md.mask.ice_levelset=icemask;
-md.mask.ocean_levelset=-icemask;
-
-%make sure that the elements that have loads are fully grounded:
-pos=find(md.solidearth.surfaceload.icethicknesschange);
-md.mask.ocean_levelset(md.mesh.elements(pos,:))=1;
-
-%make sure wherever there is an ice load, that the mask is set to ice:
-pos=find(md.solidearth.surfaceload.icethicknesschange);
-md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
+md.mask.ocean_levelset=oceanmask;
 % }}}
 
-md.solidearth.settings.ocean_area_scaling=0;
-
-%Geometry for the bed, arbitrary: 
-md.geometry.bed=-ones(md.mesh.numberofvertices,1);
+%time stepping: 
+md.timestepping.start_time=0;
+md.timestepping.time_step=1;
+md.timestepping.final_time=10;
+
+%masstransport:
+md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
+md.initialization.str=0;
 
 %Materials: 
@@ -57 +68 @@
 
 %Solution parameters
+md.cluster.np=3;
 md.solidearth.settings.reltol=NaN;
 md.solidearth.settings.abstol=1e-3;
-md.solidearth.settings.computesealevelchange=1;
-
-% max number of iteration reverted back to 10 (i.e., the original default value)
-md.solidearth.settings.maxiter=10;
-
-%eustatic + rigid + elastic + rotation run:
-md.solidearth.settings.rigid=1;
+md.solidearth.settings.sealevelloading=1;
+md.solidearth.settings.isgrd=1;
+md.solidearth.settings.ocean_area_scaling=0;
+md.solidearth.settings.grdmodel=1;
+
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
 md.solidearth.settings.rotation=1;
-
-%transient settings: 
-md.timestepping.start_time=0;
-md.timestepping.final_time=10;
-md.timestepping.time_step=1;
-md.transient.isslc=1;
-md.transient.issmb=0;
-md.transient.ismasstransport=0;
+md.solidearth.settings.viscous=0;
+
+%Physics: 
+md.transient.issmb=0; 
 md.transient.isstressbalance=0;
 md.transient.isthermal=0;
-dh=md.solidearth.surfaceload.icethicknesschange;
-deltathickness=zeros(md.mesh.numberofelements+1,10);
-for i=1:10,
-        deltathickness(1:end-1,i)=dh*i;
+md.transient.ismasstransport=1;
+md.transient.isslc=1;
+md.solidearth.requested_outputs={'Sealevel'};
+
+dh=md.masstransport.spcthickness;
+deltathickness=zeros(md.mesh.numberofvertices+1,10);
+for i=0:10,
+        deltathickness(1:end-1,i+1)=md.geometry.thickness+dh(1:end-1)*i;
 end
-deltathickness(end,:)=0:1:9;
-md.solidearth.surfaceload.icethicknesschange=deltathickness;
+deltathickness(end,:)=0:1:10;
+md.masstransport.spcthickness=deltathickness;
 
 %hack: 

issm/trunk-jpl/test/NightlyRun/test2007.m

@@ -1 @@
-%Test Name: offline solid earth solution
+%Test Name: External_OfflineSolidearthSolution
 
 %mesh earth:
 md=model;
-md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh
+load ../Data/SlcTestMesh.mat;
+md.mesh=SlcMesh; %700 km resolution mesh
 
-%Geometry for the bed, arbitrary thickness of 1000: 
+%Geometry for the bed, arbitrary
 md.geometry.bed=-ones(md.mesh.numberofvertices,1);
 md.geometry.base=md.geometry.bed;
@@ -46 +47 @@
 %Solution parameters
 md.cluster.np=3;
+md.solidearth.settings.isgrd=0;
 md.solidearth.settings.horiz=1;
 

issm/trunk-jpl/test/NightlyRun/test2008.m

@@ -1 @@
-%Test Name: External:AdditionalSolidearthSolution
+%Test Name: External_AdditionalSolidearthSolution
 
 %mesh earth:
 md=model;
-md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh
+load ../Data/SlcTestMesh.mat;
+md.mesh=SlcMesh; %700 km resolution mesh
 
 %Geometry for the bed, arbitrary thickness of 100: 
@@ -70 +71 @@
 md.solidearth.settings.reltol=NaN;
 md.solidearth.settings.abstol=1e-3;
-md.solidearth.settings.computesealevelchange=1;
+md.solidearth.settings.sealevelloading=1;
 md.solidearth.settings.isgrd=1;
 md.solidearth.settings.ocean_area_scaling=0;
@@ -86 +87 @@
 md.solidearth.settings.maxiter=10;
 
-%eustatic + rigid + elastic:
-md.solidearth.settings.rigid=1;
+%eustatic + selfattraction + elastic:
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
 md.solidearth.settings.rotation=0;

issm/trunk-jpl/test/NightlyRun/test2010.m

@@ -3 +3 @@
 %mesh earth:
 md=model;
-rad_e = 6.371012*10^3; % mean radius of Earth, km
-md.mesh=gmshplanet('radius',rad_e,'resolution',1000.0);  % km resolution
+load ../Data/SlcTestMesh.mat;
+md.mesh=SlcMesh; %700 km resolution mesh
+
+%Geometry for the bed, arbitrary thickness of 1000: 
+md.geometry.bed=-ones(md.mesh.numberofvertices,1);
+md.geometry.base=md.geometry.bed;
+md.geometry.thickness=100*ones(md.mesh.numberofvertices,1);
+md.geometry.surface=md.geometry.bed+md.geometry.thickness;
+
 
 %parameterize slc solution:
-%slc loading:  {{{
-late=sum(md.mesh.lat(md.mesh.elements),2)/3;
-longe=sum(md.mesh.long(md.mesh.elements),2)/3;
+%solidearth loading:  {{{
+md.masstransport.spcthickness=[md.geometry.thickness;0];
+md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
 
-md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
-pos=find(late <-75 & longe >0);
-md.solidearth.surfaceload.icethicknesschange(pos(6:7))=-1;
 
-md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
-md.dsl.global_average_thermosteric_sea_level_change=[0;0];
-md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
-md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
+xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
+ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3;
+ze=md.mesh.z(md.mesh.elements)*[1;1;1]/3;
+re=sqrt(xe.^2+ye.^2+ze.^2);
 
-md.solidearth.settings.ocean_area_scaling = 1;
+late=asind(ze./re);
+longe=atan2d(ye,xe);
+%greenland
+pos=find(late>60 & late<90 & longe>-75 & longe<-15);
+md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
+posice=pos;
 
 %elastic loading from love numbers:
-md.solidearth.lovenumbers=lovenumbers('maxdeg',1000);
+md.solidearth.lovenumbers=lovenumbers('maxdeg',100);
 
 %}}}
@@ -29 +38 @@
 mask=gmtmask(md.mesh.lat,md.mesh.long);
 icemask=ones(md.mesh.numberofvertices,1);
-pos=find(mask==0);
-icemask(pos)=-1;
-pos=find(sum(mask(md.mesh.elements),2)<3);
-icemask(md.mesh.elements(pos,:))=-1;
+icemask(md.mesh.elements(posice,:))=-0.5;
+
+oceanmask=-ones(md.mesh.numberofvertices,1);
+pos=find(mask==0); oceanmask(pos)=1; icemask(~pos)=1;
+
 md.mask.ice_levelset=icemask;
-md.mask.ocean_levelset=-icemask;
+md.mask.ocean_levelset=oceanmask;
 
-%make sure that the elements that have loads are fully grounded:
-pos=find(md.solidearth.surfaceload.icethicknesschange);
-md.mask.ocean_levelset(md.mesh.elements(pos,:))=1;
+% use model representation of ocen area (not the true area)
+md.solidearth.settings.ocean_area_scaling = 0;
 
-%make sure wherever there is an ice load, that the mask is set to ice:
-md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
-pos=find(md.solidearth.surfaceload.icethicknesschange);
-md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
-% }}}
+%materials
+md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
+md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
+md.initialization.str=0;
 
-%geometry {{{
-di=md.materials.rho_ice/md.materials.rho_water;
-md.geometry.thickness=ones(md.mesh.numberofvertices,1);
-md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
-md.geometry.base=md.geometry.surface-md.geometry.thickness;
-md.geometry.bed=md.geometry.base;
-% }}}
-%materials {{{
-md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
-md.materials.rheology_B=paterson(md.initialization.temperature);
-md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
-% }}}
-%Miscellaneous {{{
+md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+
+%Miscellaneous
 md.miscellaneous.name='test2010';
-% }}}
-%Solution parameters {{{
+
+%Solution parameters
 md.solidearth.settings.reltol=NaN;
 md.solidearth.settings.abstol=1e-3;
-md.solidearth.settings.computesealevelchange=1;
-% }}}
+md.solidearth.settings.sealevelloading=0;
+md.solidearth.settings.grdocean=1;
+md.solidearth.settings.isgrd=1;
+md.solidearth.settings.ocean_area_scaling=0;
+md.solidearth.settings.grdmodel=1;
+md.solidearth.settings.horiz=1;
+md.solidearth.requested_outputs={'Sealevel','SealevelBarystaticIceArea','SealevelBarystaticIceLoad','SealevelBarystaticIceMask','SealevelBarystaticIceLatbar' 'SealevelBarystaticIceLongbar'};
 
-%eustatic + rigid + elastic run:
-md.solidearth.settings.rigid=1;
+%Physics: 
+md.transient.issmb=0; 
+md.transient.isstressbalance=0;
+md.transient.isthermal=0;
+md.transient.ismasstransport=1;
+md.transient.isslc=1;
+ 
+md.timestepping.start_time=0;
+md.timestepping.time_step=1;
+md.timestepping.final_time=1;
+
+%eustatic + selfattraction + elastic + rotation run:
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
 md.solidearth.settings.rotation=1;
+md.solidearth.settings.viscous=0;
 md.cluster=generic('name',oshostname(),'np',3);
+%md.verbose=verbose('111111111');
+md=solve(md,'Transient');
 
+moi_p=md.solidearth.rotational.polarmoi;
+moi_e=md.solidearth.rotational.equatorialmoi;
+tide_love_k2=md.solidearth.lovenumbers.tk(3);
+load_love_k2=md.solidearth.lovenumbers.k(3);
+tide_love_k2secular=md.solidearth.lovenumbers.tk2secular;
 % uncomment following 2 lines for
-md=solve(md,'Sealevelchange');
-eus=md.results.SealevelchangeSolution.Bslc;
-slc=md.results.SealevelchangeSolution.Sealevel;
-moixz=md.results.SealevelchangeSolution.SealevelInertiaTensorXZ;
-moiyz=md.results.SealevelchangeSolution.SealevelInertiaTensorYZ;
-moizz=md.results.SealevelchangeSolution.SealevelInertiaTensorZZ;
+eus=md.results.TransientSolution.Bslc;
+slc=md.results.TransientSolution.Sealevel;
+moixz=md.results.TransientSolution.SealevelInertiaTensorXZ / (1/(1-tide_love_k2/tide_love_k2secular) * (1+load_love_k2)/(moi_p-moi_e) );
+moiyz=md.results.TransientSolution.SealevelInertiaTensorYZ / (1/(1-tide_love_k2/tide_love_k2secular) * (1+load_love_k2)/(moi_p-moi_e) );
+moizz=md.results.TransientSolution.SealevelInertiaTensorZZ / ( -(1+load_love_k2)/moi_p);
+
+areaice=md.results.TransientSolution.SealevelBarystaticIceArea;
+loadice=md.results.TransientSolution.SealevelBarystaticIceLoad;
 
 % analytical moi => just checking FOR ICE only!!! {{{
 % ...have to mute ** slc induced MOI in Tria.cpp ** prior to the comparison
-%rad_e = rad_e*1e3; % now in meters
-%areas=GetAreasSphericalTria(md.mesh.elements,md.mesh.lat,md.mesh.long,rad_e);
-%lat=late*pi/180; lon=longe*pi/180;
-%moi_xz = sum(-md.materials.rho_freshwater.*md.solidearth.surfaceload.icethicknesschange.*areas.*rad_e^2.*sin(lat).*cos(lat).*cos(lon));
-%moi_yz = sum(-md.materials.rho_freshwater.*md.solidearth.surfaceload.icethicknesschange.*areas.*rad_e^2.*sin(lat).*cos(lat).*sin(lon));
+rad_e = md.solidearth.planetradius;
+
+lat=md.results.TransientSolution.SealevelBarystaticIceLatbar*pi/180;
+lon=md.results.TransientSolution.SealevelBarystaticIceLongbar*pi/180;
+moi_xz = sum(-loadice.*areaice.*rad_e^2.*sin(lat).*cos(lat).*cos(lon));
+moi_yz = sum(-loadice.*areaice.*rad_e^2.*sin(lat).*cos(lat).*sin(lon));
+moi_zz = sum(-loadice.*areaice.*rad_e^2.*(1.0-sin(lat).^2));
+theoretical_value_check=[moixz/moi_xz moiyz/moi_yz moizz/moizz]
 % }}}
 
 %Fields and tolerances to track changes
 field_names     ={'eus','slc','moixz','moiyz','moizz'};
-field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13};
+field_tolerances={1e-13,1e-13,1e-13,1e-13};
 field_values={eus,slc,moixz,moiyz,moizz};
 

issm/trunk-jpl/test/NightlyRun/test2011.m

@@ -1 @@
-%Test Name: EarthSlc
+%Test Name: SlcBarystatic
 
 %mesh earth:
 md=model;
-md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500.); %700 km resolution mesh
+load ../Data/SlcTestMesh.mat;
+md.mesh=SlcMesh; %700 km resolution mesh
 
 %Geometry for the bed, arbitrary thickness of 1000: 
@@ -70 +71 @@
 
 %Miscellaneous
-md.miscellaneous.name='test2002';
+md.miscellaneous.name='test2011';
 
 %Solution parameters
@@ -76 +77 @@
 md.solidearth.settings.reltol=NaN;
 md.solidearth.settings.abstol=1e-3;
-md.solidearth.settings.computesealevelchange=1;
+md.solidearth.settings.sealevelloading=1;
 md.solidearth.settings.isgrd=1;
 md.solidearth.settings.ocean_area_scaling=0;
@@ -94 +95 @@
 
 %eustatic + rigid + elastic + rotation run:
-md.solidearth.settings.rigid=1;
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
-md.solidearth.settings.rotation=1;
+md.solidearth.settings.rotation=0;
+md.solidearth.settings.viscous=0;
 md=solve(md,'tr');
 
@@ -106 +108 @@
 weights=md.results.TransientSolution(1).SealevelBarystaticIceWeights; 
 dH=md.results.TransientSolution(1).DeltaIceThickness(md.mesh.elements); 
-dHavg=sum(dH.*weights);
+dHavg=sum(dH.*weights,2);
 oceanarea=sum(md.results.TransientSolution(1).SealevelBarystaticOceanArea);
 bslc2=-sum(dHavg.*area)*md.materials.rho_ice/md.materials.rho_water/oceanarea;
@@ -113 +115 @@
 field_names={'BarystaticIce','BarystaticIce2','BarystaticIceDiff'};
 field_tolerances={1e-13,1e-13,1e-13};
-field_values={bslc,bslc2,bslc2-bslc2};
+field_values={bslc,bslc2,bslc2-bslc};

issm/trunk-jpl/test/NightlyRun/test2090.m

@@ -3 +3 @@
 %mesh earth:
 md=model;
-md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh
+load ../Data/SlcTestMesh.mat;
+md.mesh=SlcMesh; %700 km resolution mesh
 
 %Geometry for the bed, arbitrary thickness of 1000: 
@@ -73 +74 @@
 time1=md.timestepping.start_time:md.timestepping.time_step:md.timestepping.final_time;
 md.masstransport.spcthickness=repmat(md.masstransport.spcthickness, [1 length(time1)]);
-md.masstransport.spcthickness(1:end-1,2:end)=0;
+md.masstransport.spcthickness(1:end-1,3:end)=0;
 md.masstransport.spcthickness(end,:)=time1;
 
@@ -92 +93 @@
 
 %Solution parameters
-md.cluster.np=1;
+md.cluster.np=3;
 md.solidearth.settings.reltol=NaN;
 md.solidearth.settings.abstol=1e-3;
-md.solidearth.settings.computesealevelchange=1;
+md.solidearth.settings.sealevelloading=1;
 md.solidearth.settings.isgrd=1;
 md.solidearth.settings.ocean_area_scaling=0;
 md.solidearth.settings.grdmodel=1;
 md.solidearth.settings.timeacc=md.timestepping.time_step;
-md.solidearth.settings.degacc=.01;
+md.solidearth.settings.degacc=.1;
 
 %Physics:
@@ -113 +114 @@
 md.solidearth.settings.maxiter=10;
 
-%eustatic + rigid + elastic + rotation run:
-md.solidearth.settings.rigid=1;
+%eustatic + selfattraction + elastic + rotation run:
+md.solidearth.settings.selfattraction=1;
 md.solidearth.settings.elastic=1;
 md.solidearth.settings.viscous=1;
@@ -127 +128 @@
 end
 
+hold on
+plot(md.timestepping.start_time:md.timestepping.time_step:(md.timestepping.final_time-md.timestepping.time_step), [B(916,:)]);
+
 %Fields and tolerances to track changes
 field_names={'Sealevel','Bed'};