Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp	(revision 27052)
+++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp	(revision 27053)
@@ -264,7 +264,4 @@
 
 			parameters->AddObject(new DoubleParam(SolidearthSettingsTimeAccEnum,timeacc));
-			parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,ndeg,precomputednt));
-			parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,ndeg,precomputednt));
-			parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,ndeg,precomputednt));
 
 			if (rotation){
@@ -274,10 +271,4 @@
 				iomodel->FetchData(&love_pmtf_colinear,&dummy,&precomputednt,"md.solidearth.lovenumbers.pmtf_colinear");
 				iomodel->FetchData(&love_pmtf_ortho,&dummy,&precomputednt,"md.solidearth.lovenumbers.pmtf_ortho");
-
-				parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionColinearEnum,love_pmtf_colinear,1,precomputednt));
-				parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionOrthogonalEnum,love_pmtf_ortho,1,precomputednt));
-				parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,ndeg,precomputednt));
-				parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,ndeg,precomputednt));
-				parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,ndeg,precomputednt));
 			}
 
@@ -595,13 +586,22 @@
 				xDelete<IssmDouble>(G_viscoelastic);
 				xDelete<IssmDouble>(G_viscoelastic_local);
+				xDelete<IssmDouble>(G_viscoelastic_interpolated);
 				xDelete<IssmDouble>(U_viscoelastic);
 				xDelete<IssmDouble>(U_viscoelastic_local);
+				xDelete<IssmDouble>(U_viscoelastic_interpolated);
 				if(horiz){
 					xDelete<IssmDouble>(H_viscoelastic);
 					xDelete<IssmDouble>(H_viscoelastic_local);
+					xDelete<IssmDouble>(H_viscoelastic_interpolated);
 				}
 				if(rotation){
+					xDelete<IssmDouble>(Love_th2_interpolated);
+					xDelete<IssmDouble>(Love_tk2_interpolated);
+					if (horiz) xDelete<IssmDouble>(Love_tl2_interpolated);
 					xDelete<IssmDouble>(love_pmtf_colinear);
 					xDelete<IssmDouble>(love_pmtf_ortho);
+					xDelete<IssmDouble>(Pmtf_col_interpolated);
+					xDelete<IssmDouble>(Pmtf_ortho_interpolated);
+					xDelete<IssmDouble>(Pmtf_z_interpolated);
 
 				}
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 27052)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 27053)
@@ -6371,4 +6371,5 @@
 	bool computerotation = false;
 	bool computeviscous = false;
+	bool computesatgravi = false;
 	int  horiz;
 	bool istime=true;
@@ -6384,4 +6385,5 @@
 	IssmDouble* __restrict__ tide_love_l  = NULL;
 	IssmDouble* __restrict__ LoveRotRSL   = NULL;
+	IssmDouble* __restrict__ LoveRotSatG  = NULL;
 	IssmDouble* __restrict__ LoveRotU     = NULL;
 	IssmDouble* __restrict__ LoveRothoriz = NULL;
@@ -6395,7 +6397,9 @@
 	IssmDouble* tide_love_l  = NULL;
 	IssmDouble* LoveRotRSL   = NULL;
+	IssmDouble* LoveRotSatG  = NULL;
 	IssmDouble* LoveRotU     = NULL;
 	IssmDouble* LoveRothoriz = NULL;
 	IssmDouble*  Grot        = NULL;
+
 	IssmDouble* GUrot        = NULL;
 	IssmDouble* GNrot        = NULL;
@@ -6587,4 +6591,5 @@
 		GUrot       = xNewZeroInit<IssmDouble>(3*3*nt);
 
+
 		if (horiz){
 			GErot=xNewZeroInit<IssmDouble>(3*3*nt);
@@ -6679,4 +6684,5 @@
 		this->inputs->SetArrayInput(SealevelchangeViscousUEnum,this->lid,viscousU,3*nt);
 		this->parameters->SetParam(0,SealevelchangeViscousIndexEnum);
+
 		if(horiz){
 			viscousN=xNewZeroInit<IssmDouble>(3*nt);
@@ -6691,34 +6697,73 @@
 	#ifdef _HAVE_RESTRICT_
 	delete G;
+	delete G_gravi_precomputed;
 	if(computeelastic){
 		delete GU;
+		delete G_viscoelastic_precomputed;
+		delete U_viscoelastic_precomputed;
 		if(horiz){
 			delete GN;
 			delete GE;
+			delete H_viscoelastic_precomputed;
 		}
 		if(computerotation){
 			delete Grot;
 			delete GUrot;
+			delete tide_love_h;
+			delete tide_love_k;
+			delete tide_love_l;
+			delete LoveRotRSL;
+			delete LoveRotU;
 			if (horiz){
 				delete GNrot;
 				delete GErot;
-			}
+				delete LoveRothoriz;
+			}
+		}
+	}
+	if(computeviscous){
+		delete viscousRSL;
+		delete viscousU;
+
+		if(horiz){
+			delete viscousN;
+			delete viscousE;
 		}
 	}
 	#else
 	xDelete(G);
+	xDelete(G_gravi_precomputed);
 	if(computeelastic){
 		xDelete(GU);
+		xDelete(G_viscoelastic_precomputed);
+		xDelete(U_viscoelastic_precomputed);
 		if(horiz){
 			xDelete(GN);
 			xDelete(GE);
+			xDelete(H_viscoelastic_precomputed);
 		}
 		if(computerotation){
 			xDelete(Grot);
 			xDelete(GUrot);
+			xDelete(tide_love_h);
+			xDelete(tide_love_k);
+			xDelete(tide_love_l);
+			xDelete(LoveRotRSL);
+			xDelete(LoveRotU);
 			if (horiz){
 				xDelete(GNrot);
 				xDelete(GErot);
-			}
+				xDelete(LoveRothoriz);
+			}
+		}
+		
+	}
+	if(computeviscous){
+		xDelete(viscousRSL);
+		xDelete(viscousU);
+
+		if(horiz){
+			xDelete(viscousN);
+			xDelete(viscousE);
 		}
 	}
@@ -7262,7 +7307,11 @@
 	}
 	xDelete<IssmDouble*>(Gsubel);
+	xDelete<IssmDouble>(G_gravi_precomputed);
 	if(computeelastic){
 		xDelete<IssmDouble*>(GUsubel);
+		xDelete<IssmDouble>(G_viscoelastic_precomputed);
+		xDelete<IssmDouble>(U_viscoelastic_precomputed);
 		if(horiz){
+			xDelete<IssmDouble>(H_viscoelastic_precomputed);
 			xDelete<IssmDouble*>(GNsubel); 
 			xDelete<IssmDouble*>(GEsubel);
@@ -7278,4 +7327,5 @@
 	/*Inputs:*/
 	IssmDouble* viscousRSL=NULL;
+	IssmDouble* viscousSG=NULL;
 	IssmDouble* viscousU=NULL;
 	IssmDouble* viscousN=NULL;
@@ -7293,4 +7343,5 @@
 
 		this->inputs->GetArrayPtr(SealevelchangeViscousRSLEnum,this->lid,&viscousRSL,&dummy);
+		this->inputs->GetArrayPtr(SealevelchangeViscousSGEnum,this->lid,&viscousSG,&dummy);
 		this->inputs->GetArrayPtr(SealevelchangeViscousUEnum,this->lid,&viscousU,&dummy);
 		if(horiz){
@@ -7301,4 +7352,5 @@
 		for(int i=0;i<NUMVERTICES;i++){
 			viscousRSL[i*viscousnumsteps+newindex+offset]=(1-lincoeff)*viscousRSL[i*viscousnumsteps+newindex]+lincoeff*viscousRSL[i*viscousnumsteps+newindex+1];
+			viscousSG[i*viscousnumsteps+newindex+offset]=(1-lincoeff)*viscousSG[i*viscousnumsteps+newindex]+lincoeff*viscousSG[i*viscousnumsteps+newindex+1];
 			viscousU[i*viscousnumsteps+newindex+offset]=(1-lincoeff)*viscousU[i*viscousnumsteps+newindex]+lincoeff*viscousU[i*viscousnumsteps+newindex+1];
 			if(horiz){
@@ -7309,4 +7361,11 @@
 
 	}
+
+	/*Free ressources*/
+	xDelete<IssmDouble>(viscousRSL);
+	xDelete<IssmDouble>(viscousSG);
+	xDelete<IssmDouble>(viscousU);
+	xDelete<IssmDouble>(viscousN);
+	xDelete<IssmDouble>(viscousE);
 
 }
@@ -7393,4 +7452,6 @@
 	barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
 
+	xDelete<IssmDouble>(areae);
+
 }
 /*}}}*/
@@ -7426,4 +7487,9 @@
 	}
 
+	xDelete<IssmDouble>(G);
+	xDelete<IssmDouble>(Grot);
+	for (int l=0;l<SLGEOM_NUMLOADS;l++) xDelete<IssmDouble>(Gsub[l]);
+	xDelete<IssmDouble>(Gsub[SLGEOM_NUMLOADS]);
+
 	return;
 } /*}}}*/
@@ -7468,4 +7534,5 @@
 	IssmDouble SealevelGrd[3]={0,0,0};
 	IssmDouble RSLGrd[3]={0,0,0};
+	IssmDouble SGGrd[3]={0,0,0}; //Satellite Gravimetry
 	IssmDouble UGrd[3]={0,0,0};
 	IssmDouble NGrd[3]={0,0,0};
@@ -7478,4 +7545,5 @@
 	IssmDouble* GN=NULL;
 	IssmDouble* Grot=NULL;
+	IssmDouble* GSGrot=NULL;
 	IssmDouble* GUrot=NULL;
 	IssmDouble* GNrot=NULL;
@@ -7494,4 +7562,5 @@
 	bool percpu=false;
 	bool planethasocean=false;
+	bool SatelliteGravi=false;
 
 	this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
@@ -7501,4 +7570,5 @@
 	this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
 	this->parameters->FindParam(&planethasocean,SolidearthSettingsGrdOceanEnum);
+	this->parameters->FindParam(&SatelliteGravi,SolidearthSettingsSatelliteGraviEnum);
 
 	if(sal){
@@ -7528,4 +7598,5 @@
 				this->inputs->GetArrayPtr(SealevelchangeGrotEnum,this->lid,&Grot,&size);
 				this->inputs->GetArrayPtr(SealevelchangeGUrotEnum,this->lid,&GUrot,&size);
+				if (SatelliteGravi) this->inputs->GetArrayPtr(SealevelchangeGSatGravirotEnum,this->lid,&GSGrot,&size);
 				if (horiz){
 					this->inputs->GetArrayPtr(SealevelchangeGErotEnum,this->lid,&GErot,&size);
@@ -7535,4 +7606,5 @@
 		}
 		this->SealevelchangeGxL(&RSLGrd[0],G, Gsub, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
+		if (SatelliteGravi) this->SealevelchangeGxL(&SGGrd[0],G, Gsub, GSGrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousSGEnum,computefuture=true);
 
 		if(elastic){
@@ -7554,13 +7626,40 @@
 
 	/*Create geoid: */
-	for(int i=0;i<NUMVERTICES;i++)SealevelGrd[i]=UGrd[i]+RSLGrd[i];
+	for(int i=0;i<NUMVERTICES;i++){
+		SealevelGrd[i]=UGrd[i]+RSLGrd[i];
+		if (SatelliteGravi) SGGrd[i]+=UGrd[i];
+	}
 
 	/*Create inputs*/
 	this->AddInput(SealevelGRDEnum,SealevelGrd,P1Enum);
 	this->AddInput(BedGRDEnum,UGrd,P1Enum);
+	if (SatelliteGravi) this->AddInput(SatGraviGRDEnum,SGGrd,P1Enum);
 	if(horiz){
 		this->AddInput(BedNorthGRDEnum,NGrd,P1Enum);
 		this->AddInput(BedEastGRDEnum,EGrd,P1Enum);
 	}
+
+	/*Free ressources*/
+
+	xDelete<IssmDouble>(G);
+	xDelete<IssmDouble>(GU);
+	xDelete<IssmDouble>(GE);
+	xDelete<IssmDouble>(GN);
+	xDelete<IssmDouble>(Grot);
+	xDelete<IssmDouble>(GSGrot);
+	xDelete<IssmDouble>(GUrot);
+	xDelete<IssmDouble>(GNrot);
+	xDelete<IssmDouble>(GErot);
+
+	for (int l=0;l<SLGEOM_NUMLOADS;l++){
+		xDelete<IssmDouble>(Gsub[l]);
+		xDelete<IssmDouble>(GUsub[l]);
+		xDelete<IssmDouble>(GNsub[l]);
+		xDelete<IssmDouble>(GEsub[l]);
+	}
+	xDelete<IssmDouble>(Gsub[SLGEOM_NUMLOADS]);
+	xDelete<IssmDouble>(GUsub[SLGEOM_NUMLOADS]);
+	xDelete<IssmDouble>(GNsub[SLGEOM_NUMLOADS]);
+	xDelete<IssmDouble>(GEsub[SLGEOM_NUMLOADS]);
 
 } /*}}}*/
@@ -7577,4 +7676,10 @@
 	int viscousindex=0; //important
 	int viscousnumsteps=1; //important
+
+	IssmDouble* grdfieldinterp=NULL;
+	IssmDouble* viscoustimes=NULL;
+	IssmDouble  final_time;
+	IssmDouble  lincoeff;
+	IssmDouble  timeacc;
 
 	this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
@@ -7661,10 +7766,4 @@
 		// 3*: subtract from viscous stack the grdfield that has already been accounted for so we don't add it again at the next time step
 
-		IssmDouble* grdfieldinterp=NULL;
-		IssmDouble* viscoustimes=NULL;
-		IssmDouble  final_time;
-		IssmDouble  lincoeff;
-		IssmDouble  timeacc;
-
 		this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
 		this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
@@ -7710,6 +7809,4 @@
 		}
 
-		/*Free allocatoins:*/
-		xDelete<IssmDouble>(viscoustimes);
 	} 
 	/*}}}*/
@@ -7727,4 +7824,15 @@
 	}
 
+
+	/*Free allocatoins:*/
+	xDelete<IssmDouble>(grdfield);
+	if (computeviscous){
+		xDelete<IssmDouble>(viscoustimes);
+		xDelete<IssmDouble>(viscousfield);
+		if(computefuture){
+			xDelete<IssmDouble>(grdfieldinterp);
+		}
+	}
+
 } /*}}}*/
 
@@ -7740,4 +7848,5 @@
 
 } /*}}}*/
+
 #endif
 
Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 27052)
+++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 27053)
@@ -430,4 +430,10 @@
 	}
 
+	/*Free ressources*/
+	delete loads;
+	delete oldsealevelloads;
+	delete oceanareas;
+	delete subelementoceanareas;
+
 }
 /*}}}*/
@@ -607,4 +613,9 @@
 	#endif
 
+	xDelete<IssmDouble>(xxe);
+	xDelete<IssmDouble>(yye);
+	xDelete<IssmDouble>(zze);
+	xDelete<IssmDouble>(areae);
+
 	return;
 
@@ -672,4 +683,9 @@
 	femmodel->parameters->AddObject(new DoubleVecParam(ZzeEnum,zze,nel));
 	femmodel->parameters->AddObject(new DoubleVecParam(AreaeEnum,areae,nel));
+
+	xDelete<IssmDouble>(xxe);
+	xDelete<IssmDouble>(yye);
+	xDelete<IssmDouble>(zze);
+	xDelete<IssmDouble>(areae);
 
 	return slgeom;
@@ -864,5 +880,11 @@
 	xDelete<IssmDouble>(m2);
 	xDelete<IssmDouble>(m3);
+
+	xDelete<IssmDouble>(pmtf_col);
+	xDelete<IssmDouble>(pmtf_ortho);
+	xDelete<IssmDouble>(pmtf_z);
 	if (viscous){
+		xDelete<IssmDouble>(viscoustimes);
+																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																																			xDelete<IssmDouble>(viscouspolarmotion);
 		if (computefuture){
 			xDelete<IssmDouble>(m1interp);
@@ -935,8 +957,9 @@
 		/*free allocations:*/
 		xDelete<IssmDouble>(viscoustimes);
-	}
-
-
-}
+		if (rotation) xDelete<IssmDouble>(viscouspolarmotion);
+	}
+
+
+} /*}}}*/
 void        ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
 
@@ -974,4 +997,5 @@
 	loadcopy=vloadcopy->ToMPISerial();
 
+	delete vloadcopy;
 	return loadcopy;