Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp	(revision 27130)
+++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp	(revision 27131)
@@ -124,4 +124,7 @@
 	IssmDouble  planetradius=0;
 	IssmDouble  planetarea=0;
+	IssmDouble  constant=0;
+	IssmDouble  rho_earth;
+	int		isgrd=0;
 	bool		selfattraction=false;
 	bool		elastic=false;
@@ -137,4 +140,6 @@
 	int     numoutputs;
 	char**  requestedoutputs = NULL;
+	int* recvcounts = NULL;
+	int* displs=NULL;
 
 	/*transition vectors: */
@@ -158,4 +163,5 @@
 	if(isexternal) parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.external.nature",SolidearthExternalNatureEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.runfrequency",SolidearthSettingsRunFrequencyEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.degacc",SolidearthSettingsDegreeAccuracyEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.reltol",SolidearthSettingsReltolEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.abstol",SolidearthSettingsAbstolEnum));
@@ -181,7 +187,7 @@
 	parameters->AddObject(new DoubleParam(CumBslcHydroEnum,0.0));
 	parameters->AddObject(new DoubleParam(CumGmtslcEnum,0.0));
-
 	/*compute planet area and plug into parameters:*/
 	iomodel->FetchData(&planetradius,"md.solidearth.planetradius");
+	iomodel->FetchData(&rho_earth,"md.materials.earth_density");
 	planetarea=4*M_PI*planetradius*planetradius;
 	parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea));
@@ -245,5 +251,6 @@
 
 	parameters->FindParam(&grdmodel,GrdModelEnum);
-	if(grdmodel==ElasticEnum){
+	parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum);
+	if(grdmodel==ElasticEnum && isgrd){
 		/*Deal with elasticity {{{*/
 		iomodel->FetchData(&selfattraction,"md.solidearth.settings.selfattraction");
@@ -259,4 +266,7 @@
 		}
 
+		//default values
+		nt=1;
+		ntimesteps=1;
 		/*love numbers: */
 		if(viscous || elastic){
@@ -273,7 +283,7 @@
 
 			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));
+			//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){
@@ -286,14 +296,11 @@
 				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));
+				//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));
 			}
 
 			parameters->AddObject(new DoubleMatParam(LoveTimeFreqEnum,love_timefreq,precomputednt,1));
 			parameters->AddObject(new BoolParam(LoveIsTimeEnum,love_istime));
-
-			/*Free allocations:*/
-			xDelete<IssmDouble>(love_timefreq);
 
 			// AD performance is sensitive to calls to ensurecontiguous.
@@ -320,15 +327,8 @@
 				}
 			}
-			else {
-				ntimesteps=1;
-				nt=1;
-			}
-
 #ifdef _HAVE_AD_
-			G_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
 			U_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
 			if(horiz)H_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
 #else
-			G_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
 			U_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
 			if(horiz)H_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
@@ -336,99 +336,97 @@
 		}
 		if(selfattraction){
-#ifdef _HAVE_AD_
-			G_gravi=xNew<IssmDouble>(M,"t");
-#else
-			G_gravi=xNew<IssmDouble>(M);
-#endif
-		}
-
-		if(rotation) parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
-		if(selfattraction){
-
 			/*compute combined legendre + love number (elastic green function):*/
 			m=DetermineLocalSize(M,IssmComm::GetComm());
 			GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
+#ifdef _HAVE_AD_
+			G_gravi=xNew<IssmDouble>(M,"t");
+			G_gravi_local=xNew<IssmDouble>(m,"t");
+			G_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
+			G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
+#else
+			G_gravi=xNew<IssmDouble>(M);
+			G_gravi_local=xNew<IssmDouble>(m);
+			G_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
+			G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
+#endif
 		}
 		if(viscous | elastic){
 #ifdef _HAVE_AD_
-			G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
 			U_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
 			if(horiz)H_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
 #else
-			G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
 			U_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
 			if(horiz)H_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
 #endif
 		}
-		if(selfattraction){
-#ifdef _HAVE_AD_
-			G_gravi_local=xNew<IssmDouble>(m,"t");
-#else
-			G_gravi_local=xNew<IssmDouble>(m);
-#endif
-		}
-
+
+		if(rotation) parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
+		constant=3/rho_earth/planetarea;
 		if(selfattraction){
 			for(int i=lower_row;i<upper_row;i++){
 				IssmDouble alpha,x;
 				alpha= reCast<IssmDouble>(i)*degacc * M_PI / 180.0;
-				G_gravi_local[i-lower_row]= .5/sin(alpha/2.0);
-			}
-		}
-		if(viscous | elastic){
-			for(int i=lower_row;i<upper_row;i++){
-				IssmDouble alpha,x;
-				alpha= reCast<IssmDouble>(i)*degacc * M_PI / 180.0;
-
-				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_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; 
-				}
-
-				IssmDouble Pn = 0.; 
-				IssmDouble Pn1 = 0.; 
-				IssmDouble Pn2 = 0.; 
-				IssmDouble Pn_p = 0.; 
-				IssmDouble Pn_p1 = 0.; 
-				IssmDouble Pn_p2 = 0.; 
-
-				for (int n=0;n<ndeg;n++) {
-
-					/*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
-					if(n==0){
-						Pn=1; 
-						Pn_p=0; 
+				G_gravi_local[i-lower_row]= constant*.5/sin(alpha/2.0);
+			}
+			if(viscous | elastic){
+				for(int i=lower_row;i<upper_row;i++){
+					IssmDouble alpha,x;
+					alpha= reCast<IssmDouble>(i)*degacc * M_PI / 180.0;
+
+					for(int t=0;t<ntimesteps;t++){
+						G_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (1.0+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; 
 					}
-					else if(n==1){ 
-						Pn = cos(alpha); 
-						Pn_p = 1; 
+
+					IssmDouble Pn = 0.; 
+					IssmDouble Pn1 = 0.; 
+					IssmDouble Pn2 = 0.; 
+					IssmDouble Pn_p = 0.; 
+					IssmDouble Pn_p1 = 0.; 
+					IssmDouble Pn_p2 = 0.; 
+
+					for (int n=0;n<ndeg;n++) {
+
+						/*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
+						if(n==0){
+							Pn=1; 
+							Pn_p=0; 
+						}
+						else if(n==1){ 
+							Pn = cos(alpha); 
+							Pn_p = 1; 
+						}
+						else{
+							Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
+							Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
+						}
+						Pn2=Pn1; Pn1=Pn;
+						Pn_p2=Pn_p1; Pn_p1=Pn_p;
+
+						for(int t=0;t<ntimesteps;t++){
+							IssmDouble deltalove_G;
+							IssmDouble deltalove_U;
+
+							deltalove_G = (love_k[n*precomputednt+t]-love_k[(ndeg-1)*precomputednt+t]-love_h[n*precomputednt+t]+love_h[(ndeg-1)*precomputednt+t]);
+							deltalove_U = (love_h[n*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t]);
+
+							G_viscoelastic_local[(i-lower_row)*ntimesteps+t] += constant*deltalove_G*Pn;		                // gravitational potential 
+							U_viscoelastic_local[(i-lower_row)*ntimesteps+t] += constant*deltalove_U*Pn;		                // vertical (up) displacement 
+							if(horiz)H_viscoelastic_local[(i-lower_row)*ntimesteps+t] += constant*sin(alpha)*love_l[n*precomputednt+t]*Pn_p;		// horizontal displacements 
+						}
 					}
-					else{
-						Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
-						Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
+				}
+			}
+			else { //just copy G_gravi into G_viscoelastic
+				for(int i=lower_row;i<upper_row;i++){
+					for(int t=0;t<ntimesteps;t++){
+						G_viscoelastic_local[(i-lower_row)*ntimesteps+t]= G_gravi_local[i-lower_row];
 					}
-					Pn2=Pn1; Pn1=Pn;
-					Pn_p2=Pn_p1; Pn_p1=Pn_p;
-
-					for(int t=0;t<ntimesteps;t++){
-						IssmDouble deltalove_G;
-						IssmDouble deltalove_U;
-
-						deltalove_G = (love_k[n*precomputednt+t]-love_k[(ndeg-1)*precomputednt+t]-love_h[n*precomputednt+t]+love_h[(ndeg-1)*precomputednt+t]);
-						deltalove_U = (love_h[n*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t]);
-
-						G_viscoelastic_local[(i-lower_row)*ntimesteps+t] += deltalove_G*Pn;		                // gravitational potential 
-						U_viscoelastic_local[(i-lower_row)*ntimesteps+t] += deltalove_U*Pn;		                // vertical (up) displacement 
-						if(horiz)H_viscoelastic_local[(i-lower_row)*ntimesteps+t] += sin(alpha)*love_l[n*precomputednt+t]*Pn_p;		// horizontal displacements 
-					}
-				}
-			}
-		}
-		if(selfattraction){
-
+				}
+			}
 			/*merge G_viscoelastic_local into G_viscoelastic; U_viscoelastic_local into U_viscoelastic; H_viscoelastic_local to H_viscoelastic:{{{*/
-			int* recvcounts=xNew<int>(IssmComm::GetSize());
-			int* displs=xNew<int>(IssmComm::GetSize());
+			recvcounts=xNew<int>(IssmComm::GetSize());
+			displs=xNew<int>(IssmComm::GetSize());
 			int  rc;
 			int  offset;
@@ -436,18 +434,15 @@
 			//deal with selfattraction first: 
 			ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
-
 			/*displs: */
 			ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
-
 			/*All gather:*/
 			ISSM_MPI_Allgatherv(G_gravi_local, m, ISSM_MPI_DOUBLE, G_gravi, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
 
+			rc=m*ntimesteps;
+			offset=lower_row*ntimesteps;
+			ISSM_MPI_Allgather(&rc,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
+			ISSM_MPI_Allgather(&offset,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
+			ISSM_MPI_Allgatherv(G_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, G_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
 			if(elastic){
-				rc=m*ntimesteps;
-				offset=lower_row*ntimesteps;
-				ISSM_MPI_Allgather(&rc,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
-				ISSM_MPI_Allgather(&offset,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
-
-				ISSM_MPI_Allgatherv(G_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, G_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
 				ISSM_MPI_Allgatherv(U_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, U_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
 				if(horiz)ISSM_MPI_Allgatherv(H_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, H_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
@@ -460,7 +455,9 @@
 			/*Avoid singularity at 0: */
 			G_gravi[0]=G_gravi[1];
+			for(int t=0;t<ntimesteps;t++){
+				G_viscoelastic[t]=G_viscoelastic[ntimesteps+t];
+			}
 			if(elastic){
 				for(int t=0;t<ntimesteps;t++){
-					G_viscoelastic[t]=G_viscoelastic[ntimesteps+t];
 					U_viscoelastic[t]=U_viscoelastic[ntimesteps+t];
 					if(horiz)H_viscoelastic[t]=H_viscoelastic[ntimesteps+t];
@@ -475,5 +472,5 @@
 				G_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
 				U_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
-				if(horiz)H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
+				if(horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
 				if(rotation){
 					Pmtf_col_interpolated=xNew<IssmDouble>(nt,"t");
@@ -487,5 +484,5 @@
 				G_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
 				U_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
-				if(horiz)H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
+				if(horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
 				if(rotation){
 					Pmtf_col_interpolated=xNew<IssmDouble>(nt);
@@ -497,10 +494,8 @@
 				}
 #endif
-
 				for(int t=0;t<nt;t++){
 					IssmDouble lincoeff;
 					IssmDouble viscoelastic_time=t*timeacc;
 					int        timeindex2=-1;
-
 					/*Find a way to interpolate precomputed Gkernels to our solution time stepping:*/
 					if(t!=0){
@@ -521,5 +516,4 @@
 
 					for(int index=0;index<M;index++){
-
 						int timeindex=index*nt+t;
 						int timepreindex= index*ntimesteps+timeindex2;
@@ -539,53 +533,57 @@
 					}
 				}
-
-			}
-			else if(elastic){
+			}
+			else {
+
 				nt=1; //in elastic, or if we run only selfattraction, we need only one step
 #ifdef _HAVE_AD_
 				G_viscoelastic_interpolated=xNew<IssmDouble>(M,"t");
-				U_viscoelastic_interpolated=xNew<IssmDouble>(M,"t");
-				if(horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M,"t");
 #else
 				G_viscoelastic_interpolated=xNew<IssmDouble>(M);
-				U_viscoelastic_interpolated=xNew<IssmDouble>(M);
-				if(horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M);
 #endif
-
-
 				xMemCpy<IssmDouble>(G_viscoelastic_interpolated,G_viscoelastic,M);
-				xMemCpy<IssmDouble>(U_viscoelastic_interpolated,U_viscoelastic,M);
-				if(horiz) xMemCpy<IssmDouble>(H_viscoelastic_interpolated,H_viscoelastic,M);
-
-				if(rotation){ //if this cpu handles degree 2
+
+				if(elastic){
 #ifdef _HAVE_AD_
-					Pmtf_col_interpolated=xNew<IssmDouble>(1,"t");
-					Pmtf_ortho_interpolated=xNew<IssmDouble>(1,"t");
-					Pmtf_z_interpolated=xNew<IssmDouble>(1,"t");
-					Love_tk2_interpolated=xNew<IssmDouble>(1,"t");
-					Love_th2_interpolated=xNew<IssmDouble>(1,"t");
-					if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(1,"t");
+					U_viscoelastic_interpolated=xNew<IssmDouble>(M,"t");
+					if (horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M,"t");
 #else
-					Pmtf_col_interpolated=xNew<IssmDouble>(1);
-					Pmtf_ortho_interpolated=xNew<IssmDouble>(1);
-					Pmtf_z_interpolated=xNew<IssmDouble>(1);
-					Love_tk2_interpolated=xNew<IssmDouble>(1);
-					Love_th2_interpolated=xNew<IssmDouble>(1);
-					if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(1);
+					U_viscoelastic_interpolated=xNew<IssmDouble>(M);
+					if (horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M);
 #endif
-
-					Pmtf_col_interpolated[0]=love_pmtf_colinear[0];
-					Pmtf_ortho_interpolated[0]=love_pmtf_ortho[0];
-					Pmtf_z_interpolated[0]=1.0+love_k[2];
-					Love_tk2_interpolated[0]=love_tk[2];
-					Love_th2_interpolated[0]=love_th[2];
-					if (horiz) Love_tl2_interpolated[0]=love_tl[2];
+					xMemCpy<IssmDouble>(U_viscoelastic_interpolated,U_viscoelastic,M);
+					if (horiz) xMemCpy<IssmDouble>(H_viscoelastic_interpolated,H_viscoelastic,M);
+
+					if(rotation){ //if this cpu handles degree 2
+#ifdef _HAVE_AD_
+						Pmtf_col_interpolated=xNew<IssmDouble>(1,"t");
+						Pmtf_ortho_interpolated=xNew<IssmDouble>(1,"t");
+						Pmtf_z_interpolated=xNew<IssmDouble>(1,"t");
+						Love_tk2_interpolated=xNew<IssmDouble>(1,"t");
+						Love_th2_interpolated=xNew<IssmDouble>(1,"t");
+						if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(1,"t");
+#else
+						Pmtf_col_interpolated=xNew<IssmDouble>(1);
+						Pmtf_ortho_interpolated=xNew<IssmDouble>(1);
+						Pmtf_z_interpolated=xNew<IssmDouble>(1);
+						Love_tk2_interpolated=xNew<IssmDouble>(1);
+						Love_th2_interpolated=xNew<IssmDouble>(1);
+						if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(1);
+#endif
+
+						Pmtf_col_interpolated[0]=love_pmtf_colinear[0];
+						Pmtf_ortho_interpolated[0]=love_pmtf_ortho[0];
+						Pmtf_z_interpolated[0]=1.0+love_k[2];
+						Love_tk2_interpolated[0]=love_tk[2];
+						Love_th2_interpolated[0]=love_th[2];
+						if (horiz) Love_tl2_interpolated[0]=love_tl[2];
+					}
+
 				}
 			}	
-
 			/*Save our precomputed tables into parameters*/
 			parameters->AddObject(new DoubleVecParam(SealevelchangeGSelfAttractionEnum,G_gravi,M));
+			parameters->AddObject(new DoubleVecParam(SealevelchangeGViscoElasticEnum,G_viscoelastic_interpolated,M*nt));
 			if(viscous || elastic){
-				parameters->AddObject(new DoubleVecParam(SealevelchangeGViscoElasticEnum,G_viscoelastic_interpolated,M*nt));
 				parameters->AddObject(new DoubleVecParam(SealevelchangeUViscoElasticEnum,U_viscoelastic_interpolated,M*nt));
 				if(horiz)parameters->AddObject(new DoubleVecParam(SealevelchangeHViscoElasticEnum,H_viscoelastic_interpolated,M*nt));
@@ -599,9 +597,12 @@
 				}
 			}
-
 			/*free resources: */
 			xDelete<IssmDouble>(G_gravi);
 			xDelete<IssmDouble>(G_gravi_local);
+			xDelete<IssmDouble>(G_viscoelastic);
+			xDelete<IssmDouble>(G_viscoelastic_local);
+			xDelete<IssmDouble>(G_viscoelastic_interpolated);
 			if(elastic){
+				xDelete<IssmDouble>(love_timefreq);
 				xDelete<IssmDouble>(love_h);
 				xDelete<IssmDouble>(love_k);
@@ -610,28 +611,25 @@
 				xDelete<IssmDouble>(love_tk);
 				xDelete<IssmDouble>(love_tl);
-				xDelete<IssmDouble>(G_viscoelastic);
-				xDelete<IssmDouble>(G_viscoelastic_local);
+
 				xDelete<IssmDouble>(G_viscoelastic_interpolated);
 				xDelete<IssmDouble>(U_viscoelastic);
 				xDelete<IssmDouble>(U_viscoelastic_interpolated);
 				xDelete<IssmDouble>(U_viscoelastic_local);
+				xDelete<IssmDouble>(U_viscoelastic_interpolated);
 				if(horiz){
 					xDelete<IssmDouble>(H_viscoelastic);
 					xDelete<IssmDouble>(H_viscoelastic_interpolated);
 					xDelete<IssmDouble>(H_viscoelastic_local);
+					xDelete<IssmDouble>(H_viscoelastic_interpolated);
 				}
 				if(rotation){
 					xDelete<IssmDouble>(love_pmtf_colinear);
 					xDelete<IssmDouble>(love_pmtf_ortho);
-					
-					xDelete<IssmDouble>(Love_tk2_interpolated);
-					xDelete<IssmDouble>(Love_th2_interpolated);
 					xDelete<IssmDouble>(Pmtf_col_interpolated);
 					xDelete<IssmDouble>(Pmtf_ortho_interpolated);
 					xDelete<IssmDouble>(Pmtf_z_interpolated);
-					if(horiz){
-						xDelete<IssmDouble>(Love_tl2_interpolated);
-					}
-
+					xDelete<IssmDouble>(Love_tk2_interpolated);
+					xDelete<IssmDouble>(Love_th2_interpolated);
+					if (horiz) xDelete<IssmDouble>(Love_tl2_interpolated);
 				}
 			}
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 27130)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 27131)
@@ -22,4 +22,5 @@
 #include "../Inputs/ControlInput.h"
 #include "../Inputs/ArrayInput.h"
+#include "../Inputs/IntArrayInput.h"
 /*}}}*/
 #define MAXVERTICES 6 /*Maximum number of vertices per element, currently Penta, to avoid dynamic mem allocation*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 27130)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 27131)
@@ -404,5 +404,5 @@
 		virtual void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom)=0;
 		virtual void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom)=0;
-		virtual void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
+		virtual void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids)=0;
 		virtual void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
 		virtual void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 27130)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 27131)
@@ -227,5 +227,5 @@
 		void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
 		void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
-		void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
+		void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");};
 		void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
 		void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 27130)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 27131)
@@ -180,5 +180,5 @@
 		void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
 		void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
-		void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
+		void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");};
 		void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
 		void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 27130)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 27131)
@@ -187,5 +187,5 @@
 		void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
 		void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
-		void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
+		void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");};
 		void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
 		void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 27130)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 27131)
@@ -6286,5 +6286,5 @@
 }
 /*}}}*/
-void       Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/
+void       Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){ /*{{{*/
 
 	/*Declarations:{{{*/
@@ -6308,21 +6308,7 @@
 
 	#ifdef _HAVE_RESTRICT_
-	IssmDouble* __restrict__ G=NULL;
-	IssmDouble* __restrict__ GU=NULL;
-	IssmDouble* __restrict__ GN=NULL;
-	IssmDouble* __restrict__ GE=NULL;
-	IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
 	IssmDouble* __restrict__ G_gravi_precomputed=NULL;
-	IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
-	IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
 	#else
-	IssmDouble* G=NULL;
-	IssmDouble* GU=NULL;
-	IssmDouble* GN=NULL;
-	IssmDouble* GE=NULL;
-	IssmDouble* G_viscoelastic_precomputed=NULL;
 	IssmDouble* G_gravi_precomputed=NULL;
-	IssmDouble* U_viscoelastic_precomputed=NULL;
-	IssmDouble* H_viscoelastic_precomputed=NULL;
 	#endif
 
@@ -6330,4 +6316,5 @@
 	int index;
 	int M;
+	IssmDouble degacc;
 	IssmDouble doubleindex,lincoef;
 
@@ -6346,4 +6333,8 @@
 	/*Rotational:*/
 	#ifdef _HAVE_RESTRICT_
+	IssmDouble* __restrict__ Grot=NULL;
+	IssmDouble* __restrict__ GUrot=NULL;
+	IssmDouble* __restrict__ GNrot=NULL;
+	IssmDouble* __restrict__ GErot=NULL;
 	IssmDouble* __restrict__ tide_love_h  = NULL;
 	IssmDouble* __restrict__ tide_love_k  = NULL;
@@ -6352,9 +6343,11 @@
 	IssmDouble* __restrict__ LoveRotU     = NULL;
 	IssmDouble* __restrict__ LoveRothoriz = NULL;
-	IssmDouble* __restrict__  Grot        = NULL;
-	IssmDouble* __restrict__ GUrot        = NULL;
-	IssmDouble* __restrict__ GNrot        = NULL;
-	IssmDouble* __restrict__ GErot        = NULL;
+	int* __restrict__ AplhaIndex   = NULL;
+	int* __restrict__ AzimuthIndex = NULL;
 	#else
+	IssmDouble* Grot=NULL;
+	IssmDouble* GUrot=NULL;
+	IssmDouble* GNrot=NULL;
+	IssmDouble* GErot=NULL;
 	IssmDouble* tide_love_h  = NULL;
 	IssmDouble* tide_love_k  = NULL;
@@ -6363,8 +6356,6 @@
 	IssmDouble* LoveRotU     = NULL;
 	IssmDouble* LoveRothoriz = NULL;
-	IssmDouble*  Grot        = NULL;
-	IssmDouble* GUrot        = NULL;
-	IssmDouble* GNrot        = NULL;
-	IssmDouble* GErot        = NULL;
+	int* AlphaIndex   = NULL;
+	int* AzimuthIndex = NULL;
 	#endif
 
@@ -6405,19 +6396,9 @@
 
 	/*Recover precomputed green function kernels:{{{*/
-	DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGSelfAttractionEnum)); _assert_(parameter);
-	parameter->GetParameterValueByPointer((IssmDouble**)&G_gravi_precomputed,&M);
-
+	parameters->FindParam(&degacc,SolidearthSettingsDegreeAccuracyEnum);
+	M=reCast<int,IssmDouble>(180.0/degacc+1.);
+
+	DoubleVecParam* parameter;
 	if(computeelastic){
-		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
-		parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL);
-
-		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
-		parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL);
-
-		if(horiz){
-			parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
-			parameter->GetParameterValueByPointer((IssmDouble**)&H_viscoelastic_precomputed,NULL);
-		}
-
 		if(computerotation){
 			parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeTidalH2Enum)); _assert_(parameter);
@@ -6455,92 +6436,45 @@
 	}
 
-	constant=3/rho_earth/planetarea;
-
-	G=xNewZeroInit<IssmDouble>(3*nel*nt);
-	if(computeelastic){
-		GU=xNewZeroInit<IssmDouble>(3*nel*nt);
-		if(horiz){
-			GN=xNewZeroInit<IssmDouble>(3*nel*nt);
-			GE=xNewZeroInit<IssmDouble>(3*nel*nt);
-		}
-	}
-
-	for(int e=0;e<nel;e++){
-		ratioe=constant*areae[e];
-		for (int i=0;i<3;i++){
-			IssmDouble alpha;
-			IssmDouble delPhi,delLambda;
-			/*recovers info for this element and vertex:*/
-			IssmDouble late= asin(zze[e]/sqrt( pow(xxe[e],2.0)+ pow(yye[e],2.0)+ pow(zze[e],2.0)));
-			IssmDouble longe= atan2(yye[e],xxe[e]);
-
-			lati=latitude[i];
-			longi=longitude[i];
-
-			/*Computes alpha angle between centroid and current vertex, and indexes alpha in precomputed tables: */
-			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)));
-			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){
-				/*Compute azimuths, both north and east components: */
-				x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
-				if(lati==M_PI/2){
-					x=1e-12; y=1e-12;
+	AlphaIndex=xNewZeroInit<int>(3*nel);
+	if (horiz) AzimuthIndex=xNewZeroInit<int>(3*nel);
+	int intmax=pow(2,16)-1;
+
+
+	for (int i=0;i<3;i++){
+		if(lids[this->vertices[i]->lid]==this->lid){
+			for(int e=0;e<nel;e++){
+				IssmDouble alpha;
+				IssmDouble delPhi,delLambda;
+				/*recovers info for this element and vertex:*/
+				IssmDouble late= asin(zze[e]/sqrt( pow(xxe[e],2.0)+ pow(yye[e],2.0)+ pow(zze[e],2.0)));
+				IssmDouble longe= atan2(yye[e],xxe[e]);
+
+				lati=latitude[i];
+				longi=longitude[i];
+
+				/*Computes alpha angle between centroid and current vertex, and indexes alpha in precomputed tables: */
+				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)));
+				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
+				if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour
+				_assert_(index>=0 && index<M);
+
+				if(horiz){
+					/*Compute azimuths*/
+					dx=cos(lati)*sin(late)-sin(lati)*cos(late)*cos(longe-longi);
+					dy=sin(longe-longi)*cos(late);
+					//angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax]
+					AzimuthIndex[i*nel+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
 				}
-				if(lati==-M_PI/2){
-					x=1e-12; y=1e-12;
-				}
-				dx = xxe[e]-x; dy = yye[e]-y; dz = zze[e]-z;
-				N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
-				E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
-			}
-
-			for(int t=0;t<nt;t++){
-				int timeindex=i*nel*nt+e*nt+t; 
-
-				/*Rigid earth gravitational perturbation: */
-				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+0)*nt+t]*(1.-lincoef)
-						     +G_viscoelastic_precomputed[(index+1)*nt+t]*lincoef; //linear interpolation
-				}
-				G[timeindex]=G[timeindex]*ratioe;
-
-				/*Elastic components:*/
-				if(computeelastic){
-					GU[timeindex] =  ratioe *(U_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
-								 +U_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
-					if(horiz){
-						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);
-					}
-				}
-			} //for(int t=0;t<nt;t++)
+				AlphaIndex[i*nel+e]=index;
+			}			
 		} //for (int i=0;i<3;i++)
 	} //for(int e=0;e<nel;e++)
 
 	/*Add in inputs:*/
-	this->inputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,nel*3*nt);
-	if(computeelastic){
-		this->inputs->SetArrayInput(SealevelchangeGUEnum,this->lid,GU,nel*3*nt);
-		if(horiz){
-			this->inputs->SetArrayInput(SealevelchangeGNEnum,this->lid,GN,nel*3*nt);
-			this->inputs->SetArrayInput(SealevelchangeGEEnum,this->lid,GE,nel*3*nt);
-		}
-	}
+	this->inputs->SetIntArrayInput(SealevelchangeAlphaIndexEnum,this->lid,AlphaIndex,nel*3);
+	if(horiz) this->inputs->SetIntArrayInput(SealevelchangeAzimuthIndexEnum,this->lid,AzimuthIndex,nel*3);
+
 	/*}}}*/
 	/*Compute rotation kernel:{{{*/
@@ -6660,38 +6594,183 @@
 	/*Free allocations:{{{*/
 	#ifdef _HAVE_RESTRICT_
-	delete G;
-	if(computeelastic){
-		delete GU;
-		if(horiz){
-			delete GN;
-			delete GE;
-		}
-		if(computerotation){
-			delete Grot;
-			delete GUrot;
-			if (horiz){
-				delete GNrot;
-				delete GErot;
-			}
-		}
-	}
+	delete AlphaIndex;
+	if(horiz) AzimuthIndex;
+
+	if(computerotation){
+		delete Grot;
+		delete GUrot;
+		if (horiz){
+			delete GNrot;
+			delete GErot;
+		}
+	}
+
 	#else
-	xDelete(G);
-	if(computeelastic){
-		xDelete(GU);
-		if(horiz){
-			xDelete(GN);
-			xDelete(GE);
-		}
-		if(computerotation){
-			xDelete(Grot);
-			xDelete(GUrot);
-			if (horiz){
-				xDelete(GNrot);
-				xDelete(GErot);
-			}
+	xDelete<int>(AlphaIndex);
+	if(horiz){
+		xDelete<int>(AzimuthIndex);
+	}
+	if(computerotation){
+		xDelete<IssmDouble>(Grot);
+		xDelete<IssmDouble>(GUrot);
+		if (horiz){
+			xDelete<IssmDouble>(GNrot);
+			xDelete<IssmDouble>(GErot);
 		}
 	}
 	#endif
+	/*}}}*/
+	return;
+
+}
+/*}}}*/
+
+void       Tria::SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){ /*{{{*/
+
+	/*Declarations:{{{*/
+	int nel;
+	IssmDouble planetarea,planetradius;
+	IssmDouble constant,ratioe;
+	IssmDouble rho_earth;
+	IssmDouble lati,longi;
+	IssmDouble latitude[NUMVERTICES];
+	IssmDouble longitude[NUMVERTICES];
+	IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	#ifdef _HAVE_RESTRICT_
+	int** __restrict__ AlphaIndex=NULL;
+	int** __restrict__ AzimIndex=NULL;
+
+	#else
+	int** AlphaIndex=NULL;
+	int** AzimIndex=NULL;
+	#endif
+
+	/*viscoelastic green function:*/
+	int index;
+	int M;
+	IssmDouble doubleindex,lincoef, degacc;
+
+	/*Computational flags:*/
+	bool computeselfattraction = false;
+	bool computeelastic = false;
+	bool computeviscous = false;
+	int  horiz;
+	int grd, grdmodel;
+
+	bool istime=true;
+	IssmDouble timeacc=0;
+	IssmDouble start_time,final_time;
+	int  nt,precomputednt;
+	int intmax=pow(2,16)-1;
+
+	/*}}}*/
+	/*Recover parameters:{{{ */
+	rho_earth=FindParam(MaterialsEarthDensityEnum);
+	this->parameters->FindParam(&computeselfattraction,SolidearthSettingsSelfAttractionEnum);
+	this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
+	this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
+	this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
+	this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
+	this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
+	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;
+
+	/*Recover precomputed green function kernels:{{{*/
+	parameters->FindParam(&degacc,SolidearthSettingsDegreeAccuracyEnum);
+	M=reCast<int,IssmDouble>(180.0/degacc+1.);
+
+	/*}}}*/
+	/*Compute lat long of all vertices in the element:{{{*/
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	for(int i=0;i<NUMVERTICES;i++){
+		latitude[i]= asin(xyz_list[i][2]/planetradius);
+		longitude[i]= atan2(xyz_list[i][1],xyz_list[i][0]);
+	}
+	/*}}}*/
+	/*Compute green functions:{{{ */
+
+	if(computeviscous){
+		this->parameters->FindParam(&istime,LoveIsTimeEnum);
+		if(!istime)_error_("Frequency love numbers not supported yet!");
+		this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
+		this->parameters->FindParam(&start_time,TimesteppingStartTimeEnum);
+		this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
+		nt=reCast<int,IssmDouble>((final_time-start_time)/timeacc)+1;
+	}
+	else{
+		nt=1; //in elastic, or if we run only selfattraction, we need only one step
+	}
+	AlphaIndex=xNew<int*>(SLGEOM_NUMLOADS);
+	if(horiz) AzimIndex=xNew<int*>(SLGEOM_NUMLOADS);
+
+
+	//Allocate: 
+	for(int l=0;l<SLGEOM_NUMLOADS;l++){
+		int nbar=slgeom->nbar[l];
+		AlphaIndex[l]=xNewZeroInit<int>(3*nbar);
+		if(horiz) AzimIndex[l]=xNewZeroInit<int>(3*nbar);
+
+
+		for (int i=0;i<3;i++){
+			if(slgeom->lids[this->vertices[i]->lid]==this->lid){
+				for(int e=0;e<nbar;e++){
+					IssmDouble alpha;
+					IssmDouble delPhi,delLambda;
+					/*recover info for this element and vertex:*/
+					IssmDouble late= slgeom->latbarycentre[l][e]; 
+					IssmDouble longe= slgeom->longbarycentre[l][e]; 
+					late=late/180*M_PI;
+					longe=longe/180*M_PI;
+
+					lati=latitude[i];
+					longi=longitude[i];
+
+					if(horiz){
+						/*Compute azimuths*/
+						dx=cos(lati)*sin(late)-sin(lati)*cos(late)*cos(longe-longi);
+						dy=sin(longe-longi)*cos(late);
+						//angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax]
+						AzimIndex[l][i*nbar+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
+					}
+
+					/*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
+					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)));
+					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
+
+					if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour
+					if (index==M-1){ //avoids out of bound case
+						index-=1;
+						lincoef=1;
+					}
+					AlphaIndex[l][i*nbar+e]=index;
+				}
+			}
+		}
+	}
+
+	/*Save all these arrayins for each element:*/
+	for (int l=0;l<SLGEOM_NUMLOADS;l++){
+		this->inputs->SetIntArrayInput(slgeom->AlphaIndexEnum(l),this->lid,AlphaIndex[l],slgeom->nbar[l]*3);
+		if(horiz) this->inputs->SetIntArrayInput(slgeom->AzimuthIndexEnum(l),this->lid,AzimIndex[l],slgeom->nbar[l]*3);
+	}
+	/*}}}*/
+	/*Free memory:{{{*/
+	for (int l=0;l<SLGEOM_NUMLOADS;l++){
+		xDelete<int>(AlphaIndex[l]);
+		if(horiz) xDelete<int>(AzimIndex[l]);
+	}
+	xDelete<int*>(AlphaIndex);
+	if(horiz) xDelete<int*>(AzimIndex); 
+	
 	/*}}}*/
 	return;
@@ -6906,4 +6985,89 @@
 }
 /*}}}*/
+void       Tria::SealevelchangeBarystaticLoads(GrdLoads* loads,  BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/
+
+	int nel;
+
+	/*Inputs:*/
+	IssmDouble I[NUMVERTICES]; 
+	IssmDouble W[NUMVERTICES];
+	IssmDouble BP[NUMVERTICES];
+	IssmDouble* areae=NULL;
+
+	/*output: */
+	IssmDouble bslcice=0;
+	IssmDouble bslchydro=0;
+	IssmDouble bslcbp=0;
+	IssmDouble BPavg=0;
+	IssmDouble Iavg=0;
+	IssmDouble Wavg=0;
+
+	/*ice properties: */
+	IssmDouble rho_ice,rho_water,rho_freshwater;
+
+	/*recover some parameters:*/
+	this->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
+	this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
+	this->parameters->FindParam(&rho_freshwater,MaterialsRhoFreshwaterEnum);
+	this->parameters->FindParam(&areae,&nel,AreaeEnum);
+
+	/*Retrieve inputs:*/
+	Element::GetInputListOnVertices(&I[0],DeltaIceThicknessEnum);
+	Element::GetInputListOnVertices(&W[0],DeltaTwsEnum);
+	Element::GetInputListOnVertices(&BP[0],DeltaBottomPressureEnum);
+
+	for(int i=0;i<NUMVERTICES;i++){
+		Iavg+=I[i]*slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid]*slgeom->LoadArea[SLGEOM_ICE][this->lid];
+		Wavg+=W[i]*slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]*slgeom->LoadArea[SLGEOM_WATER][this->lid];
+		BPavg+=BP[i]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]*slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
+	}
+
+	/*convert from m^3 to kg:*/
+	Iavg*=rho_ice;
+	Wavg*=rho_freshwater;
+	BPavg*=rho_water;
+
+	#ifdef _ISSM_DEBUG_ 
+	this->AddInput(SealevelBarystaticIceLoadEnum,&Iavg,P0Enum);
+	this->AddInput(SealevelBarystaticHydroLoadEnum,&Wavg,P0Enum);
+	this->AddInput(SealevelBarystaticBpLoadEnum,&BPavg,P0Enum);
+	#endif
+
+	/*Compute barystatic component in kg:*/
+	// Note: Iavg, etc, already include partial area factor phi for subelement loading
+	bslcice =   -Iavg;
+	bslchydro = -Wavg;
+	bslcbp =    -BPavg;
+
+	_assert_(!xIsNan<IssmDouble>(bslcice));
+	_assert_(!xIsNan<IssmDouble>(bslchydro));
+	_assert_(!xIsNan<IssmDouble>(bslcbp));
+
+	/*Plug values into subelement load vector:*/
+	if(slgeom->issubelement[SLGEOM_ICE][this->lid]){
+		int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid];
+		loads->vsubloads[SLGEOM_ICE]->SetValue(intj,Iavg,INS_VAL);
+		Iavg=0; //avoid double counting centroid loads and subelement loads!
+	}
+	if(slgeom->issubelement[SLGEOM_WATER][this->lid]){
+		int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid];
+		loads->vsubloads[SLGEOM_WATER]->SetValue(intj,Wavg,INS_VAL);
+		Wavg=0;
+	}
+	if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
+		int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
+		loads->vsubloads[SLGEOM_OCEAN]->SetValue(intj,BPavg,INS_VAL); 
+		BPavg=0;
+	}
+	/*Plug remaining values into centroid load vector:*/
+	loads->vloads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL);
+
+	/*Keep track of barystatic contributions:*/
+	barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
+
+	/*Free ressources*/
+	xDelete<IssmDouble>(areae);
+
+}/*}}}*/
 void       Tria::SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){ /*{{{*/
 
@@ -7006,242 +7170,4 @@
 }
 /*}}}*/
-void       Tria::SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){ /*{{{*/
-
-	/*Declarations:{{{*/
-	int nel;
-	IssmDouble planetarea,planetradius;
-	IssmDouble constant,ratioe;
-	IssmDouble rho_earth;
-	IssmDouble lati,longi;
-	IssmDouble latitude[NUMVERTICES];
-	IssmDouble longitude[NUMVERTICES];
-	IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
-	IssmDouble xyz_list[NUMVERTICES][3];
-
-	#ifdef _HAVE_RESTRICT_
-	IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
-	IssmDouble* __restrict__ G_gravi_precomputed=NULL;
-	IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
-	IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
-	IssmDouble** __restrict__ Gsubel=NULL;
-	IssmDouble** __restrict__ GUsubel=NULL;
-	IssmDouble** __restrict__ GNsubel=NULL;
-	IssmDouble** __restrict__ GEsubel=NULL;
-
-	#else
-	IssmDouble* G_viscoelastic_precomputed=NULL;
-	IssmDouble* G_gravi_precomputed=NULL;
-	IssmDouble* U_viscoelastic_precomputed=NULL;
-	IssmDouble* H_viscoelastic_precomputed=NULL;
-	IssmDouble** Gsubel=NULL;
-	IssmDouble** GUsubel=NULL;
-	IssmDouble** GNsubel=NULL;
-	IssmDouble** GEsubel=NULL;
-	#endif
-
-	/*viscoelastic green function:*/
-	int index;
-	int M;
-	IssmDouble doubleindex,lincoef;
-
-	/*Computational flags:*/
-	bool computeselfattraction = false;
-	bool computeelastic = false;
-	bool computeviscous = false;
-	int  horiz;
-	int grd, grdmodel;
-
-	bool istime=true;
-	IssmDouble timeacc=0;
-	IssmDouble start_time,final_time;
-	int  nt,precomputednt;
-
-	/*}}}*/
-	/*Recover parameters:{{{ */
-	rho_earth=FindParam(MaterialsEarthDensityEnum);
-	this->parameters->FindParam(&computeselfattraction,SolidearthSettingsSelfAttractionEnum);
-	this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
-	this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
-	this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
-	this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
-	this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
-	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;
-
-	/*Recover precomputed green function kernels:{{{*/
-	DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGSelfAttractionEnum)); _assert_(parameter);
-	parameter->GetParameterValueByPointer((IssmDouble**)&G_gravi_precomputed,&M);
-
-	if(computeelastic){
-		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
-		parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL);
-
-		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
-		parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL);
-
-		if(horiz){
-			parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
-			parameter->GetParameterValueByPointer((IssmDouble**)&H_viscoelastic_precomputed,NULL);
-
-		}
-	}
-	/*}}}*/
-	/*Compute lat long of all vertices in the element:{{{*/
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	for(int i=0;i<NUMVERTICES;i++){
-		latitude[i]= asin(xyz_list[i][2]/planetradius);
-		longitude[i]= atan2(xyz_list[i][1],xyz_list[i][0]);
-	}
-	/*}}}*/
-	/*Compute green functions:{{{ */
-	constant=3/rho_earth/planetarea;
-
-	if(computeviscous){
-		this->parameters->FindParam(&istime,LoveIsTimeEnum);
-		if(!istime)_error_("Frequency love numbers not supported yet!");
-		this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
-		this->parameters->FindParam(&start_time,TimesteppingStartTimeEnum);
-		this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
-		nt=reCast<int,IssmDouble>((final_time-start_time)/timeacc)+1;
-	}
-	else{
-		nt=1; //in elastic, or if we run only selfattraction, we need only one step
-	}
-	Gsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
-	if(computeelastic){
-		GUsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
-		if(horiz){
-			GNsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
-			GEsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
-		}
-	}
-
-	//Allocate: 
-	for(int l=0;l<SLGEOM_NUMLOADS;l++){
-		int nbar=slgeom->nbar[l];
-		Gsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
-		if(computeelastic){
-			GUsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
-			if(horiz){
-				GNsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
-				GEsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
-			}
-		}
-
-		for(int e=0;e<nbar;e++){
-			ratioe=constant*slgeom->area_subel[l][e];
-			for (int i=0;i<3;i++){
-				IssmDouble alpha;
-				IssmDouble delPhi,delLambda;
-				/*recover info for this element and vertex:*/
-				IssmDouble late= slgeom->latbarycentre[l][e]; 
-				IssmDouble longe= slgeom->longbarycentre[l][e]; 
-				late=late/180*M_PI;
-				longe=longe/180*M_PI;
-
-				lati=latitude[i];
-				longi=longitude[i];
-
-				if(horiz){
-					/*Compute azimuths, both north and east components: */
-					x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
-					if(lati==90){
-						x=1e-12; y=1e-12;
-					}
-					if(lati==-90){
-						x=1e-12; y=1e-12;
-					}
-					IssmDouble xbar=planetradius*cos(late)*cos(longe);
-					IssmDouble ybar=planetradius*cos(late)*sin(longe);
-					IssmDouble zbar=planetradius*sin(late);
-
-					dx = xbar-x; dy = ybar-y; dz = zbar-z;
-					N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
-					E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
-				}
-
-				/*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
-				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)));
-				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++){
-					int timeindex=i*nbar*nt+e*nt+t; 
-
-					/*Rigid earth gravitational perturbation: */
-					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+0)*nt+t]*(1.-lincoef)
-								     +G_viscoelastic_precomputed[(index+1)*nt+t]*lincoef; //linear interpolation
-					}
-					Gsubel[l][timeindex]*=ratioe;
-
-					/*Elastic components:*/
-					if(computeelastic){
-						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] = 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);
-						}
-					}
-				}
-			}
-		}
-	}
-
-	/*Save all these arrayins for each element:*/
-	for (int l=0;l<SLGEOM_NUMLOADS;l++){
-		this->inputs->SetArrayInput(slgeom->GEnum(l),this->lid,Gsubel[l],slgeom->nbar[l]*3*nt);
-		if(computeelastic){
-			this->inputs->SetArrayInput(slgeom->GUEnum(l),this->lid,GUsubel[l],slgeom->nbar[l]*3*nt);
-			if(horiz){
-				this->inputs->SetArrayInput(slgeom->GNEnum(l),this->lid,GNsubel[l],slgeom->nbar[l]*3*nt);
-				this->inputs->SetArrayInput(slgeom->GEEnum(l),this->lid,GEsubel[l],slgeom->nbar[l]*3*nt);
-			}
-		}
-	}
-	/*}}}*/
-	/*Free memory:{{{*/
-	for (int l=0;l<SLGEOM_NUMLOADS;l++){
-		xDelete<IssmDouble>(Gsubel[l]);
-		if(computeelastic){
-			xDelete<IssmDouble>(GUsubel[l]);
-			if(horiz){
-				xDelete<IssmDouble>(GNsubel[l]);
-				xDelete<IssmDouble>(GEsubel[l]);
-			}
-		}
-	}
-	xDelete<IssmDouble*>(Gsubel);
-	if(computeelastic){
-		xDelete<IssmDouble*>(GUsubel);
-		if(horiz){
-			xDelete<IssmDouble*>(GNsubel); 
-			xDelete<IssmDouble*>(GEsubel);
-		}
-	}
-	/*}}}*/
-	return;
-
-}
-/*}}}*/
 void       Tria::SealevelchangeUpdateViscousFields(IssmDouble lincoeff, int newindex, int offset){ /*{{{*/
 
@@ -7252,5 +7178,5 @@
 	IssmDouble* viscousE=NULL;
 	int         viscousnumsteps;
-	int         dummy;
+	int         size;
 	bool        viscous=false;
 	int	    horiz=0;
@@ -7262,9 +7188,9 @@
 		this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
 
-		this->inputs->GetArrayPtr(SealevelchangeViscousRSLEnum,this->lid,&viscousRSL,&dummy);
-		this->inputs->GetArrayPtr(SealevelchangeViscousUEnum,this->lid,&viscousU,&dummy);
+		this->inputs->GetArrayPtr(SealevelchangeViscousRSLEnum,this->lid,&viscousRSL,&size);
+		this->inputs->GetArrayPtr(SealevelchangeViscousUEnum,this->lid,&viscousU,&size);
 		if(horiz){
-			this->inputs->GetArrayPtr(SealevelchangeViscousNEnum,this->lid,&viscousN,&dummy);
-			this->inputs->GetArrayPtr(SealevelchangeViscousEEnum,this->lid,&viscousE,&dummy);
+			this->inputs->GetArrayPtr(SealevelchangeViscousNEnum,this->lid,&viscousN,&size);
+			this->inputs->GetArrayPtr(SealevelchangeViscousEEnum,this->lid,&viscousE,&size);
 		}
 
@@ -7277,93 +7203,40 @@
 			}
 		}
-
-	}
-
-}
-/*}}}*/
-void       Tria::SealevelchangeBarystaticLoads(GrdLoads* loads,  BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/
-
-	int nel;
-
-	/*Inputs:*/
-	IssmDouble I[NUMVERTICES]; 
-	IssmDouble W[NUMVERTICES];
-	IssmDouble BP[NUMVERTICES];
-	IssmDouble* areae=NULL;
-
-	/*output: */
-	IssmDouble bslcice=0;
-	IssmDouble bslchydro=0;
-	IssmDouble bslcbp=0;
-	IssmDouble BPavg=0;
-	IssmDouble Iavg=0;
-	IssmDouble Wavg=0;
-
-	/*ice properties: */
-	IssmDouble rho_ice,rho_water,rho_freshwater;
-
-	/*recover some parameters:*/
-	this->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
+	}
+}
+/*}}}*/
+void       Tria::SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){ /*{{{*/
+
+	IssmDouble oceanaverage=0;
+	IssmDouble oceanarea=0;
+	IssmDouble rho_water;
+
 	this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
-	this->parameters->FindParam(&rho_freshwater,MaterialsRhoFreshwaterEnum);
-	this->parameters->FindParam(&areae,&nel,AreaeEnum);
-
-	/*Retrieve inputs:*/
-	Element::GetInputListOnVertices(&I[0],DeltaIceThicknessEnum);
-	Element::GetInputListOnVertices(&W[0],DeltaTwsEnum);
-	Element::GetInputListOnVertices(&BP[0],DeltaBottomPressureEnum);
-
+
+	/*retrieve ocean average and area:*/
 	for(int i=0;i<NUMVERTICES;i++){
-		Iavg+=I[i]*slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid];
-		Wavg+=W[i]*slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid];
-		BPavg+=BP[i]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
-	}
-
-	/*convert from m to kg/m^2:*/
-	Iavg*=rho_ice;
-	Wavg*=rho_freshwater;
-	BPavg*=rho_water;
-
+		oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
+	}
 	#ifdef _ISSM_DEBUG_ 
-	this->AddInput(SealevelBarystaticIceLoadEnum,&Iavg,P0Enum);
-	this->AddInput(SealevelBarystaticHydroLoadEnum,&Wavg,P0Enum);
-	this->AddInput(SealevelBarystaticBpLoadEnum,&BPavg,P0Enum);
+	this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum);
 	#endif
-
-	/*Compute barystatic component in kg:*/
-	// Note: Iavg, etc, already include partial area factor phi for subelement loading
-	bslcice =   -slgeom->LoadArea[SLGEOM_ICE][this->lid]*Iavg;
-	bslchydro = -slgeom->LoadArea[SLGEOM_WATER][this->lid]*Wavg;
-	bslcbp =    -slgeom->LoadArea[SLGEOM_OCEAN][this->lid]*BPavg;
-
-	_assert_(!xIsNan<IssmDouble>(bslcice));
-	_assert_(!xIsNan<IssmDouble>(bslchydro));
-	_assert_(!xIsNan<IssmDouble>(bslcbp));
-
-	/*Plug values into subelement load vector:*/
-	if(slgeom->issubelement[SLGEOM_ICE][this->lid]){
-		int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid];
-		loads->vsubloads[SLGEOM_ICE]->SetValue(intj,Iavg,INS_VAL);
-		Iavg=0; //avoid double counting centroid loads and subelement loads!
-	}
-	if(slgeom->issubelement[SLGEOM_WATER][this->lid]){
-		int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid];
-		loads->vsubloads[SLGEOM_WATER]->SetValue(intj,Wavg,INS_VAL);
-		Wavg=0;
-	}
+	oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
+
+	/*add ocean average in the global sealevelloads vector:*/
 	if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
 		int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
-		loads->vsubloads[SLGEOM_OCEAN]->SetValue(intj,BPavg,INS_VAL); 
-		BPavg=0;
-	}
-	/*Plug remaining values into centroid load vector:*/
-	loads->vloads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL);
-
-	/*Keep track of barystatic contributions:*/
-	barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
-
-	/*Free resources:*/
-	xDelete<IssmDouble>(areae);
-
+		loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water*oceanarea,INS_VAL);
+		loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL);
+	}
+	else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water*oceanarea,INS_VAL);
+
+	/*add ocean area into a global oceanareas vector:*/
+	if(!loads->sealevelloads){
+		oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
+		if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
+			int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
+			subelementoceanareas->SetValue(intj,oceanarea,INS_VAL);
+		}
+	}
 }
 /*}}}*/
@@ -7371,8 +7244,11 @@
 
 	/*sal green function:*/
+	int* AlphaIndex=NULL;
+	int* AlphaIndexsub[SLGEOM_NUMLOADS];
 	IssmDouble* G=NULL;
 	IssmDouble* Grot=NULL;
-	IssmDouble* Gsub[SLGEOM_NUMLOADS];
+	DoubleVecParam* parameter;
 	bool computefuture=false;
+	int spatial_component=0;
 
 	bool sal = false;
@@ -7390,47 +7266,12 @@
 
 	if(sal){
-		this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
-		this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size);
-		this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size);
-		this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
+		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
+		parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL);
+
+		this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
+		for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
 		if (rotation)	this->inputs->GetArrayPtr(SealevelchangeGrotEnum,this->lid,&Grot,&size);
 
-		this->SealevelchangeGxL(sealevelpercpu, G, Gsub, Grot, loads, polarmotionvector, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
-	}
-
-	return;
-} /*}}}*/
-void       Tria::SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){ /*{{{*/
-
-	IssmDouble oceanaverage=0;
-	IssmDouble oceanarea=0;
-	IssmDouble rho_water;
-
-	this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
-
-	/*retrieve ocean average and area:*/
-	for(int i=0;i<NUMVERTICES;i++){
-		oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
-	}
-	#ifdef _ISSM_DEBUG_ 
-	this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum);
-	#endif
-	oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
-
-	/*add ocean average in the global sealevelloads vector:*/
-	if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
-		int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
-		loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water,INS_VAL);
-		loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL);
-	}
-	else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
-
-	/*add ocean area into a global oceanareas vector:*/
-	if(!loads->sealevelloads){
-		oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
-		if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
-			int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
-			subelementoceanareas->SetValue(intj,oceanarea,INS_VAL);
-		}
+		this->SealevelchangeGxL(sealevelpercpu, spatial_component=0,AlphaIndex, AlphaIndexsub, NULL, NULL, G, Grot, loads, polarmotionvector, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
 	}
 
@@ -7446,16 +7287,18 @@
 	int nel,nbar;
 	bool sal = false;
+	int* AlphaIndex=NULL;
+	int* AzimIndex=NULL;
+	int* AlphaIndexsub[SLGEOM_NUMLOADS];
+	int* AzimIndexsub[SLGEOM_NUMLOADS];
+	int spatial_component=0;
 	IssmDouble* G=NULL;
 	IssmDouble* GU=NULL;
-	IssmDouble* GE=NULL;
-	IssmDouble* GN=NULL;
+	IssmDouble* GH=NULL;
 	IssmDouble* Grot=NULL;
 	IssmDouble* GUrot=NULL;
 	IssmDouble* GNrot=NULL;
 	IssmDouble* GErot=NULL;
-	IssmDouble* Gsub[SLGEOM_NUMLOADS];
-	IssmDouble* GUsub[SLGEOM_NUMLOADS];
-	IssmDouble* GNsub[SLGEOM_NUMLOADS];
-	IssmDouble* GEsub[SLGEOM_NUMLOADS];
+
+	DoubleVecParam* parameter;
 	bool computefuture=false;
 
@@ -7476,25 +7319,20 @@
 
 	if(sal){
-
-		this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
-		this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size);
-		this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size);
-		this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
+		this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
+		for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
+
+		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
+		parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL);
 
 		if(elastic){
-			this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&size);
-			this->inputs->GetArrayPtr(SealevelchangeGUsubelIceEnum,this->lid,&GUsub[SLGEOM_ICE],&size);
-			this->inputs->GetArrayPtr(SealevelchangeGUsubelHydroEnum,this->lid,&GUsub[SLGEOM_WATER],&size);
-			this->inputs->GetArrayPtr(SealevelchangeGUsubelOceanEnum,this->lid,&GUsub[SLGEOM_OCEAN],&size);
+			parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
+			parameter->GetParameterValueByPointer((IssmDouble**)&GU,NULL);
+
 			if(horiz){
-				this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size);
-				this->inputs->GetArrayPtr(SealevelchangeGNsubelIceEnum,this->lid,&GNsub[SLGEOM_ICE],&size);
-				this->inputs->GetArrayPtr(SealevelchangeGNsubelHydroEnum,this->lid,&GNsub[SLGEOM_WATER],&size);
-				this->inputs->GetArrayPtr(SealevelchangeGNsubelOceanEnum,this->lid,&GNsub[SLGEOM_OCEAN],&size);
-
-				this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&size);
-				this->inputs->GetArrayPtr(SealevelchangeGEsubelIceEnum,this->lid,&GEsub[SLGEOM_ICE],&size);
-				this->inputs->GetArrayPtr(SealevelchangeGEsubelHydroEnum,this->lid,&GEsub[SLGEOM_WATER],&size);
-				this->inputs->GetArrayPtr(SealevelchangeGEsubelOceanEnum,this->lid,&GEsub[SLGEOM_OCEAN],&size);
+				this->inputs->GetIntArrayPtr(SealevelchangeAzimuthIndexEnum,this->lid,&AzimIndex,&size);
+				for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AzimuthIndexEnum(l),this->lid,&AzimIndexsub[l],&size);
+
+				parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
+				parameter->GetParameterValueByPointer((IssmDouble**)&GH,NULL);
 			}
 			if (rotation) {
@@ -7507,11 +7345,12 @@
 			}
 		}
-		this->SealevelchangeGxL(&RSLGrd[0],G, Gsub, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
+
+		this->SealevelchangeGxL(&RSLGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL,G, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
 
 		if(elastic){
-			this->SealevelchangeGxL(&UGrd[0],GU, GUsub, GUrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
+			this->SealevelchangeGxL(&UGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL, GU, GUrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
 			if(horiz){
-				this->SealevelchangeGxL(&NGrd[0],GN, GNsub, GNrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
-				this->SealevelchangeGxL(&EGrd[0],GE, GEsub, GErot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
+				this->SealevelchangeGxL(&NGrd[0],spatial_component=1,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GNrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
+				this->SealevelchangeGxL(&EGrd[0],spatial_component=2,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GErot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
 			}
 		}
@@ -7538,16 +7377,29 @@
 
 } /*}}}*/
-void       Tria::SealevelchangeGxL(IssmDouble* grdfieldout, IssmDouble* G, IssmDouble** Gsub, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/
+void       Tria::SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/
 
 	//This function performs the actual convolution between Green functions and surface Loads for a particular grd field
 
 	IssmDouble* grdfield=NULL;
-	int i,e,l,t,nbar;
+	int i,e,l,t,a, index, nbar;
+	bool rotation=false;
+	IssmDouble* Centroid_loads=NULL;
+	IssmDouble* Centroid_loads_copy=NULL;
+	IssmDouble* Subelement_loads[SLGEOM_NUMLOADS];
+	IssmDouble* Subelement_loads_copy[SLGEOM_NUMLOADS];
+	IssmDouble* horiz_projection=NULL;
+	IssmDouble* horiz_projectionsub[SLGEOM_NUMLOADS];
+	int nt=1; //important, ensures there is a defined value if computeviscous is false
+
+	//viscous
 	bool computeviscous=false;
-	bool rotation=false;
-	IssmDouble* viscousfield=NULL;
-	int nt=1; //important, ensures there is a defined value if computeviscous is false
 	int viscousindex=0; //important
 	int viscousnumsteps=1; //important
+	IssmDouble* viscousfield=NULL;
+	IssmDouble* grdfieldinterp=NULL;
+	IssmDouble* viscoustimes=NULL;
+	IssmDouble  final_time;
+	IssmDouble  lincoeff;
+	IssmDouble  timeacc;
 
 	this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
@@ -7556,7 +7408,12 @@
 		this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
 		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(viscousenum,this->lid,&viscousfield,NULL);
 		if(computefuture) {
 			nt=viscousnumsteps-viscousindex+2; //number of time steps remaining to reach final_time, +1 is sufficient with no adaptative time stepping, +2 necessary otherwise; we assume the safe choice here for the sake of simplicity
 			if (nt>viscousnumsteps) nt=viscousnumsteps;
+			grdfieldinterp=xNewZeroInit<IssmDouble>(3*viscousnumsteps); 
 		}
 		else nt=1;
@@ -7578,53 +7435,101 @@
 	}
 
-
-	if(loads->sealevelloads){ // general case: loads + sealevel loads
-		for(i=0;i<NUMVERTICES;i++) {
-			if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
-			for (e=0;e<nel;e++){
-				for(t=0;t<nt;t++){
-					int index=i*nel*viscousnumsteps+e*viscousnumsteps+t;
-					grdfield[i*nt+t]+=G[index]*(loads->sealevelloads[e]+loads->loads[e]);
+	//Determine loads /*{{{*/
+	Centroid_loads=xNewZeroInit<IssmDouble>(nel);
+	for (e=0;e<nel;e++){
+		Centroid_loads[e]=loads->loads[e];
+	}
+	for(l=0;l<SLGEOM_NUMLOADS;l++){
+		nbar=slgeom->nbar[l];
+		Subelement_loads[l]=xNewZeroInit<IssmDouble>(nbar);
+		for (e=0;e<nbar;e++){
+			Subelement_loads[l][e]=(loads->subloads[l][e]);
+		}
+	}
+	if(loads->sealevelloads){
+		for (e=0;e<nel;e++){
+			Centroid_loads[e]+=(loads->sealevelloads[e]);
+		}
+		nbar=slgeom->nbar[SLGEOM_OCEAN];
+		for (e=0;e<nbar;e++){
+			Subelement_loads[SLGEOM_OCEAN][e]+=(loads->subsealevelloads[e]);
+		}
+	}
+
+	//Copy loads if dealing with a horizontal component: the result will need to be projected against the North or East axis for each vertex
+	if (spatial_component!=0){
+		horiz_projection=xNewZeroInit<IssmDouble>(3*nel);
+		Centroid_loads_copy=xNewZeroInit<IssmDouble>(nel);
+		for (e=0;e<nel;e++){
+			Centroid_loads_copy[e]=Centroid_loads[e];
+		}
+
+		for(l=0;l<SLGEOM_NUMLOADS;l++){
+			nbar=slgeom->nbar[l];
+			Subelement_loads_copy[l]=xNewZeroInit<IssmDouble>(nbar);
+			horiz_projectionsub[l]=xNewZeroInit<IssmDouble>(3*nbar);
+			for (e=0;e<nbar;e++){
+				Subelement_loads_copy[l][e]=Subelement_loads[l][e];
+			}
+		}
+	}
+	/*}}}*/
+
+	//Convolution
+	for(i=0;i<NUMVERTICES;i++) { /*{{{*/
+		if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
+
+		if (spatial_component!=0){ //horizontals /*{{{*/
+			//GxL needs to be projected on the right axis before summation into the grd field
+			//here we apply the projection scalar to the load prior to the actual convolution loop for more efficiency
+			if (spatial_component==1){ //north
+				for (e=0;e<nel;e++){
+					horiz_projection[i*nel+e]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[i*nel+e])/65535.0); // 65535=2^16-1 = max value of 16 bits unsigned int
 				}
-			}
+				for(l=0;l<SLGEOM_NUMLOADS;l++){
+					nbar=slgeom->nbar[l];
+					for (e=0;e<nbar;e++){
+						horiz_projectionsub[l][i*nbar+e]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][i*nbar+e])/65535.0);;
+					}
+				}
+			}
+			else if (spatial_component==2){ //east
+				for (e=0;e<nel;e++){
+					horiz_projection[i*nel+e]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[i*nel+e])/65535.0);
+				}
+				for(l=0;l<SLGEOM_NUMLOADS;l++){
+					nbar=slgeom->nbar[l];
+					for (e=0;e<nbar;e++){
+						horiz_projectionsub[l][i*nbar+e]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][i*nbar+e])/65535.0);;
+					}
+				}
+			}
+			for (e=0;e<nel;e++) Centroid_loads[e]=Centroid_loads_copy[e]*horiz_projection[i*nel+e];
 			for(l=0;l<SLGEOM_NUMLOADS;l++){
 				nbar=slgeom->nbar[l];
 				for (e=0;e<nbar;e++){
-					for(t=0;t<nt;t++){
-						int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
-						grdfield[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
-					}
+					Subelement_loads[l][e]=Subelement_loads_copy[l][e]*horiz_projectionsub[l][i*nbar+e];
 				}
-				if(l==SLGEOM_OCEAN){
-					for (e=0;e<nbar;e++){
-						for(t=0;t<nt;t++){
-							int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
-							grdfield[i*nt+t]+=Gsub[l][index]*(loads->subsealevelloads[e]);
-						}
-					}
+			}
+		} /*}}}*/
+
+		for (e=0;e<nel;e++){
+			for(t=0;t<nt;t++){
+				a=AlphaIndex[i*nel+e];
+				grdfield[i*nt+t]+=G[a*viscousnumsteps+t]*Centroid_loads[e];
+			}
+		}
+		for(l=0;l<SLGEOM_NUMLOADS;l++){
+			nbar=slgeom->nbar[l];
+			for (e=0;e<nbar;e++){
+				for(t=0;t<nt;t++){
+					a=AlphaIndexsub[l][i*nbar+e];
+					grdfield[i*nt+t]+=G[a*viscousnumsteps+t]*Subelement_loads[l][e];
 				}
 			}
 		}
-	}
-	else{  //this is the initial convolution where only loads are provided
-		for(i=0;i<NUMVERTICES;i++) {
-			if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
-			for (e=0;e<nel;e++){
-				for(t=0;t<nt;t++){
-					int index=i*nel*viscousnumsteps+e*viscousnumsteps+t;
-					grdfield[i*nt+t]+=G[index]*(loads->loads[e]);
-				}
-			}
-			for(l=0;l<SLGEOM_NUMLOADS;l++){
-				nbar=slgeom->nbar[l];
-				for (e=0;e<nbar;e++){
-					for(t=0;t<nt;t++){
-						int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
-						grdfield[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
-					}
-				}
-			}
-		}
-	}
+	} /*}}}*/
+
+
 
 	if(computeviscous){ /*{{{*/
@@ -7634,17 +7539,6 @@
 		// 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);
-		this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
-		this->inputs->GetArrayPtr(viscousenum,this->lid,&viscousfield,NULL);
 		/* Map new grdfield generated by present-day loads onto viscous time vector*/
 		if(computefuture){
-			grdfieldinterp=xNewZeroInit<IssmDouble>(3*viscousnumsteps); 
 			//viscousindex time and first time step of grdfield coincide, so just copy that value
 			for(int i=0;i<NUMVERTICES;i++){
@@ -7658,5 +7552,6 @@
 					for(int i=0;i<NUMVERTICES;i++){
 						if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
-						grdfieldinterp[i*viscousnumsteps+t]=  (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]+lincoeff*grdfield[i*nt+(t-viscousindex)];
+						grdfieldinterp[i*viscousnumsteps+t] = (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]
+											 +lincoeff*grdfield[i*nt+(t-viscousindex)];
 					}
 				}
@@ -7672,9 +7567,11 @@
 		/*update viscous stack with future deformation from present load: */
 		if(computefuture){
-			for(int t=viscousnumsteps-1;t>=viscousindex;t--){
+			for(int t=viscousnumsteps-1;t>=viscousindex;t--){ //we need to go backwards so as not to zero out viscousfield[i*viscousnumsteps+viscousindex] until the end
 				for(int i=0;i<NUMVERTICES;i++){
 					if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
 					//offset viscousfield to remove all deformation that has already been added
-					viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*viscousnumsteps+t]-grdfieldinterp[i*viscousnumsteps+viscousindex]-viscousfield[i*viscousnumsteps+viscousindex];
+					viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*viscousnumsteps+t]
+									  -grdfieldinterp[i*viscousnumsteps+viscousindex]
+									  -viscousfield[i*viscousnumsteps+viscousindex];
 				}
 			}
@@ -7685,7 +7582,4 @@
 			xDelete<IssmDouble>(grdfieldinterp);
 		}
-
-		/*Free allocatoins:*/
-		xDelete<IssmDouble>(viscoustimes);
 	} 
 	/*}}}*/
@@ -7702,6 +7596,22 @@
 		for(i=0;i<NUMVERTICES;i++) grdfieldout[i]=grdfield[i*nt+0];
 	}
-	/*Free resources:*/
+	//free resources
 	xDelete<IssmDouble>(grdfield);
+	xDelete<IssmDouble>(Centroid_loads);
+	for(l=0;l<SLGEOM_NUMLOADS;l++) xDelete<IssmDouble>(Subelement_loads[l]);
+	if (spatial_component!=0){
+		xDelete<IssmDouble>(horiz_projection);
+		xDelete<IssmDouble>(Centroid_loads_copy);
+		for(l=0;l<SLGEOM_NUMLOADS;l++) {
+			xDelete<IssmDouble>(Subelement_loads_copy[l]);
+			xDelete<IssmDouble>(horiz_projectionsub[l]);
+		}
+	}
+	if (computeviscous){
+		xDelete<IssmDouble>(viscoustimes);
+		if (computefuture){
+			xDelete<IssmDouble>(grdfieldinterp);
+		}
+	}
 
 } /*}}}*/
@@ -7709,4 +7619,5 @@
 void       Tria::SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
 
+	offset*=slgeom->LoadArea[SLGEOM_OCEAN][this->lid]; //kg.m^-2 to kg 
 	if(slgeom->isoceanin[this->lid]){
 		if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 27130)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 27131)
@@ -171,5 +171,5 @@
 		#ifdef _HAVE_SEALEVELCHANGE_
 		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       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids);
 		void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom);
 		void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae);
@@ -244,5 +244,5 @@
 		void           UpdateConstraintsExtrudeFromBase(void);
 		void           UpdateConstraintsExtrudeFromTop(void);
-		void           SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector,SealevelGeometry* slgeom, int nel, bool percpu,int stackenum,bool computefuture);
+		void           SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture);
 		/*}}}*/
 
Index: /issm/trunk-jpl/src/c/classes/GrdLoads.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/GrdLoads.cpp	(revision 27130)
+++ /issm/trunk-jpl/src/c/classes/GrdLoads.cpp	(revision 27131)
@@ -90,7 +90,8 @@
 void GrdLoads::SHDegree2Coefficients(IssmDouble* deg2coeff, FemModel* femmodel, SealevelGeometry* slgeom){ /*{{{*/
 
-	IssmDouble area,re, S;
+	IssmDouble re, S;
 	int ylmindex, intj;
 	IssmDouble deg2coeff_local[5];
+	//IssmDouble area;
 
 	femmodel->parameters->FindParam(&re,SolidearthPlanetRadiusEnum);
@@ -106,9 +107,8 @@
 		if(sealevelloads) S+=sealevelloads[element->Sid()];
 		if(S!=0){
-			element->Element::GetInputValue(&area,AreaEnum);
 
 			for (int c=0;c<5;c++){ //degree l=2 has 2*l+1=5 coefficients: 2,0; 2,1cos; 2,1sin; 2,2cos; 2,2sin
 				ylmindex=(4+c)*slgeom->localnel+element->lid; // starting at index=l^2
-				deg2coeff_local[c] += S*area/re/re*slgeom->Ylm[ylmindex];
+				deg2coeff_local[c] += S/re/re*slgeom->Ylm[ylmindex];
 			}
 		}
@@ -121,8 +121,8 @@
 				if(i==SLGEOM_OCEAN && sealevelloads) S+=subsealevelloads[intj];
 				if(S!=0){
-					area=slgeom->area_subel[i][intj];
+					//area=slgeom->area_subel[i][intj];
 					for (int c=0;c<5;c++){ //degree l=2 has 2*l+1=5 coefficients
 						ylmindex=(4+c)*slgeom->localnel+element->lid; // starting at index=l^2
-						deg2coeff_local[c] += S*area/re/re*slgeom->Ylm_subel[i][ylmindex];
+						deg2coeff_local[c] += S/re/re*slgeom->Ylm_subel[i][ylmindex];
 					}
 				}
Index: /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp	(revision 27130)
+++ /issm/trunk-jpl/src/c/classes/SealevelGeometry.cpp	(revision 27131)
@@ -115,11 +115,11 @@
 
 } /*}}}*/
-int SealevelGeometry::GEnum(int l){ /*{{{*/
+int SealevelGeometry::AlphaIndexEnum(int l){ /*{{{*/
 
 	int output = -1;
 	switch(l){
-		case SLGEOM_OCEAN: output=SealevelchangeGsubelOceanEnum; break;
-		case SLGEOM_ICE:   output=SealevelchangeGsubelIceEnum;   break;
-		case SLGEOM_WATER: output=SealevelchangeGsubelHydroEnum; break;
+		case SLGEOM_OCEAN: output=SealevelchangeAlphaIndexOceanEnum; break;
+		case SLGEOM_ICE:   output=SealevelchangeAlphaIndexIceEnum;   break;
+		case SLGEOM_WATER: output=SealevelchangeAlphaIndexHydroEnum; break;
 		default: _error_("not supported");
 	}
@@ -127,35 +127,11 @@
 
 } /*}}}*/
-int SealevelGeometry::GUEnum(int l){ /*{{{*/
+int SealevelGeometry::AzimuthIndexEnum(int l){ /*{{{*/
 
 	int output = -1;
 	switch(l){
-		case SLGEOM_OCEAN: output=SealevelchangeGUsubelOceanEnum; break;
-		case SLGEOM_ICE:   output=SealevelchangeGUsubelIceEnum;   break;
-		case SLGEOM_WATER: output=SealevelchangeGUsubelHydroEnum; break;
-		default: _error_("not supported");
-	}
-	return output;
-
-} /*}}}*/
-int SealevelGeometry::GNEnum(int l){ /*{{{*/
-
-	int output = -1;
-	switch(l){
-		case SLGEOM_OCEAN: output=SealevelchangeGNsubelOceanEnum; break;
-		case SLGEOM_ICE:   output=SealevelchangeGNsubelIceEnum;   break;
-		case SLGEOM_WATER: output=SealevelchangeGNsubelHydroEnum; break;
-		default: _error_("not supported");
-	}
-
-	return output;
-} /*}}}*/
-int SealevelGeometry::GEEnum(int l){ /*{{{*/
-
-	int output = -1;
-	switch(l){
-		case SLGEOM_OCEAN: output=SealevelchangeGEsubelOceanEnum; break;
-		case SLGEOM_ICE:   output=SealevelchangeGEsubelIceEnum;   break;
-		case SLGEOM_WATER: output=SealevelchangeGEsubelHydroEnum; break;
+		case SLGEOM_OCEAN: output=SealevelchangeAzimuthIndexOceanEnum; break;
+		case SLGEOM_ICE:   output=SealevelchangeAzimuthIndexIceEnum;   break;
+		case SLGEOM_WATER: output=SealevelchangeAzimuthIndexHydroEnum; break;
 		default: _error_("not supported");
 	}
Index: /issm/trunk-jpl/src/c/classes/SealevelGeometry.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/SealevelGeometry.h	(revision 27130)
+++ /issm/trunk-jpl/src/c/classes/SealevelGeometry.h	(revision 27131)
@@ -45,8 +45,6 @@
 		void InitializeMappingsAndBarycentres(void);
 		void Assemble(void);
-		int GEnum(int l);
-		int GUEnum(int l);
-		int GNEnum(int l);
-		int GEEnum(int l);
+		int AlphaIndexEnum(int l);
+		int AzimuthIndexEnum(int l);
 		void BuildSphericalHarmonics(void);
 };
Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 27130)
+++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 27131)
@@ -29,4 +29,5 @@
 void ivins_deformation_core(FemModel* femmodel);
 IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel);
+void slc_geometry_cleanup(SealevelGeometry* slgeom, FemModel* femmodel);
 /*}}}*/
 
@@ -83,5 +84,5 @@
 
 	/*Free resources:*/
-	delete slgeom;
+	slc_geometry_cleanup(slgeom, femmodel);
 }
 /*}}}*/
@@ -361,5 +362,5 @@
 		//Conserve ocean mass: 
 		oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, totaloceanarea);
-		ConserveOceanMass(femmodel,loads,barycontrib->Total()/totaloceanarea - oceanaverage,slgeom);
+		ConserveOceanMass(femmodel,loads,barycontrib->Total()/totaloceanarea -oceanaverage,slgeom);
 
 		//broadcast sea level loads 
@@ -443,5 +444,4 @@
 	delete oceanareas;
 	xDelete<IssmDouble>(sealevelpercpu); 
-
 }
 /*}}}*/
@@ -588,4 +588,5 @@
 	IssmDouble* areae  = NULL;
 	int  nel;
+	int* lids;
 	int  grdmodel=0;
 
@@ -603,8 +604,21 @@
 	ElementCoordinatesx(&xxe,&yye,&zze,&areae,femmodel->elements);
 
+
+	/*Compute element ids, used to speed up computations in convolution phase:{{{*/
+	lids=xNew<int>(femmodel->vertices->Size());
+
+	for(Object* & object : femmodel->elements->objects){
+		Element*   element=xDynamicCast<Element*>(object);
+		for(int i=0;i<3;i++){
+			lids[element->vertices[i]->lid]=element->lid;
+		}
+	}
+
+	/*}}}*/
+
 	/*Run sealevel geometry routine in elements:*/
 	for(Object* & object : femmodel->elements->objects){
 		Element*   element=xDynamicCast<Element*>(object);
-		element->SealevelchangeGeometryInitial(xxe,yye,zze,areae);
+		element->SealevelchangeGeometryInitial(xxe,yye,zze,areae,lids);
 	}
 
@@ -625,4 +639,5 @@
 	xDelete<IssmDouble>(zze);
 	xDelete<IssmDouble>(areae);
+	xDelete(lids);
 
 	return;
@@ -642,9 +657,11 @@
 	int nel;
 	int  grdmodel=0;
+	int isgrd=0;
 	SealevelGeometry* slgeom=NULL;
 
 	/*early return?:*/
 	femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
-	if(grdmodel==IvinsEnum) return NULL;
+	femmodel->parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum);
+	if(grdmodel!=ElasticEnum || !isgrd) return NULL;
 
 	/*retrieve parameters:*/
@@ -701,4 +718,22 @@
 
 }/*}}}*/
+void slc_geometry_cleanup(SealevelGeometry* slgeom, FemModel* femmodel){  /*{{{*/
+	int  grdmodel=0;
+	int isgrd=0;
+	int horiz=0;
+
+	/*early return?:*/
+	femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
+	femmodel->parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum);
+	femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
+	if(grdmodel!=ElasticEnum || !isgrd) return;
+
+	for (int l=0;l<SLGEOM_NUMLOADS;l++){
+		femmodel->inputs->DeleteInput(slgeom->AlphaIndexEnum(l));
+		if (horiz) femmodel->inputs->DeleteInput(slgeom->AzimuthIndexEnum(l));
+	}
+
+	delete slgeom;
+} /*}}}*/
 
 /*subroutines:*/
@@ -756,7 +791,4 @@
 	Vector<IssmDouble>* vsealevelloadsvolume=loads->vsealevelloads->Duplicate();
 	Vector<IssmDouble>* vsubsealevelloadsvolume=loads->vsubsealevelloads->Duplicate();
-
-	vsealevelloadsvolume->PointwiseMult(loads->vsealevelloads,oceanareas);
-	vsubsealevelloadsvolume->PointwiseMult(loads->vsubsealevelloads,suboceanareas); 
 
 	vsealevelloadsvolume->Sum(&sealevelloadsaverage);
@@ -965,4 +997,5 @@
 		/*free allocations:*/
 		xDelete<IssmDouble>(viscoustimes);
+		if (rotation) 	xDelete<IssmDouble>(viscouspolarmotion);
 	}
 
Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 27130)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 27131)
@@ -289,5 +289,7 @@
 syn keyword cConstant LoveInnerCoreBoundaryEnum
 syn keyword cConstant LoveComplexComputationEnum
-syn keyword cConstant LoveIntStepsPerLayerEnum
+syn keyword cConstant LoveQuadPrecisionEnum
+syn keyword cConstant LoveMinIntegrationStepsEnum
+syn keyword cConstant LoveMaxIntegrationdrEnum
 syn keyword cConstant LoveKernelsEnum
 syn keyword cConstant LoveMu0Enum
@@ -301,4 +303,7 @@
 syn keyword cConstant LoveUnderflowTolEnum
 syn keyword cConstant LovePostWidderThresholdEnum
+syn keyword cConstant LoveDebugEnum
+syn keyword cConstant LoveHypergeomNZEnum
+syn keyword cConstant LoveHypergeomNAlphaEnum
 syn keyword cConstant MassFluxSegmentsEnum
 syn keyword cConstant MassFluxSegmentsPresentEnum
@@ -387,4 +392,6 @@
 syn keyword cConstant SolidearthSettingsElasticEnum
 syn keyword cConstant SolidearthSettingsViscousEnum
+syn keyword cConstant SolidearthSettingsSatelliteGraviEnum
+syn keyword cConstant SolidearthSettingsDegreeAccuracyEnum
 syn keyword cConstant SealevelchangeGeometryDoneEnum
 syn keyword cConstant SealevelchangeViscousNumStepsEnum
@@ -409,4 +416,7 @@
 syn keyword cConstant LoveTimeFreqEnum
 syn keyword cConstant LoveIsTimeEnum
+syn keyword cConstant LoveHypergeomZEnum
+syn keyword cConstant LoveHypergeomTable1Enum
+syn keyword cConstant LoveHypergeomTable2Enum
 syn keyword cConstant SealevelchangeGSelfAttractionEnum
 syn keyword cConstant SealevelchangeGViscoElasticEnum
@@ -851,4 +861,5 @@
 syn keyword cConstant SealevelEnum
 syn keyword cConstant SealevelGRDEnum
+syn keyword cConstant SatGraviGRDEnum
 syn keyword cConstant SealevelBarystaticMaskEnum
 syn keyword cConstant SealevelBarystaticIceMaskEnum
@@ -890,25 +901,19 @@
 syn keyword cConstant SealevelUNorthEsaEnum
 syn keyword cConstant SealevelchangeIndicesEnum
-syn keyword cConstant SealevelchangeGEnum
-syn keyword cConstant SealevelchangeGUEnum
-syn keyword cConstant SealevelchangeGEEnum
-syn keyword cConstant SealevelchangeGNEnum
+syn keyword cConstant SealevelchangeAlphaIndexEnum
+syn keyword cConstant SealevelchangeAzimuthIndexEnum
 syn keyword cConstant SealevelchangeGrotEnum
+syn keyword cConstant SealevelchangeGSatGravirotEnum
 syn keyword cConstant SealevelchangeGUrotEnum
 syn keyword cConstant SealevelchangeGNrotEnum
 syn keyword cConstant SealevelchangeGErotEnum
-syn keyword cConstant SealevelchangeGsubelOceanEnum
-syn keyword cConstant SealevelchangeGUsubelOceanEnum
-syn keyword cConstant SealevelchangeGEsubelOceanEnum
-syn keyword cConstant SealevelchangeGNsubelOceanEnum
-syn keyword cConstant SealevelchangeGsubelIceEnum
-syn keyword cConstant SealevelchangeGUsubelIceEnum
-syn keyword cConstant SealevelchangeGEsubelIceEnum
-syn keyword cConstant SealevelchangeGNsubelIceEnum
-syn keyword cConstant SealevelchangeGsubelHydroEnum
-syn keyword cConstant SealevelchangeGUsubelHydroEnum
-syn keyword cConstant SealevelchangeGEsubelHydroEnum
-syn keyword cConstant SealevelchangeGNsubelHydroEnum
+syn keyword cConstant SealevelchangeAlphaIndexOceanEnum
+syn keyword cConstant SealevelchangeAlphaIndexIceEnum
+syn keyword cConstant SealevelchangeAlphaIndexHydroEnum
+syn keyword cConstant SealevelchangeAzimuthIndexOceanEnum
+syn keyword cConstant SealevelchangeAzimuthIndexIceEnum
+syn keyword cConstant SealevelchangeAzimuthIndexHydroEnum
 syn keyword cConstant SealevelchangeViscousRSLEnum
+syn keyword cConstant SealevelchangeViscousSGEnum
 syn keyword cConstant SealevelchangeViscousUEnum
 syn keyword cConstant SealevelchangeViscousNEnum
@@ -1295,4 +1300,5 @@
 syn keyword cConstant DoubleArrayInputEnum
 syn keyword cConstant ArrayInputEnum
+syn keyword cConstant IntArrayInputEnum
 syn keyword cConstant DoubleExternalResultEnum
 syn keyword cConstant DoubleMatArrayParamEnum
@@ -1410,4 +1416,6 @@
 syn keyword cConstant LovePMTF1tEnum
 syn keyword cConstant LovePMTF2tEnum
+syn keyword cConstant LoveYiEnum
+syn keyword cConstant LoveRhsEnum
 syn keyword cConstant LoveSolutionEnum
 syn keyword cConstant MINIEnum
@@ -1624,4 +1632,5 @@
 syn keyword cType Cfsurfacesquare
 syn keyword cType Channel
+syn keyword cType classes
 syn keyword cType Constraint
 syn keyword cType Constraints
@@ -1630,6 +1639,6 @@
 syn keyword cType ControlInput
 syn keyword cType Covertree
+syn keyword cType DatasetInput
 syn keyword cType DataSetParam
-syn keyword cType DatasetInput
 syn keyword cType Definition
 syn keyword cType DependentObject
@@ -1644,6 +1653,6 @@
 syn keyword cType ElementInput
 syn keyword cType ElementMatrix
+syn keyword cType Elements
 syn keyword cType ElementVector
-syn keyword cType Elements
 syn keyword cType ExponentialVariogram
 syn keyword cType ExternalResult
@@ -1652,9 +1661,10 @@
 syn keyword cType Friction
 syn keyword cType Gauss
+syn keyword cType GaussianVariogram
+syn keyword cType gaussobjects
 syn keyword cType GaussPenta
 syn keyword cType GaussSeg
 syn keyword cType GaussTetra
 syn keyword cType GaussTria
-syn keyword cType GaussianVariogram
 syn keyword cType GenericExternalResult
 syn keyword cType GenericOption
@@ -1665,4 +1675,5 @@
 syn keyword cType Input
 syn keyword cType Inputs
+syn keyword cType IntArrayInput
 syn keyword cType IntInput
 syn keyword cType IntMatParam
@@ -1672,4 +1683,5 @@
 syn keyword cType IssmDirectApplicInterface
 syn keyword cType IssmParallelDirectApplicInterface
+syn keyword cType krigingobjects
 syn keyword cType Load
 syn keyword cType Loads
@@ -1682,4 +1694,5 @@
 syn keyword cType Matice
 syn keyword cType Matlitho
+syn keyword cType matrixobjects
 syn keyword cType MatrixParam
 syn keyword cType Misfit
@@ -1694,6 +1707,6 @@
 syn keyword cType Observations
 syn keyword cType Option
+syn keyword cType Options
 syn keyword cType OptionUtilities
-syn keyword cType Options
 syn keyword cType Param
 syn keyword cType Parameters
@@ -1709,11 +1722,11 @@
 syn keyword cType Regionaloutput
 syn keyword cType Results
+syn keyword cType Riftfront
 syn keyword cType RiftStruct
-syn keyword cType Riftfront
 syn keyword cType SealevelGeometry
 syn keyword cType Seg
 syn keyword cType SegInput
+syn keyword cType Segment
 syn keyword cType SegRef
-syn keyword cType Segment
 syn keyword cType SpcDynamic
 syn keyword cType SpcStatic
@@ -1734,8 +1747,4 @@
 syn keyword cType Vertex
 syn keyword cType Vertices
-syn keyword cType classes
-syn keyword cType gaussobjects
-syn keyword cType krigingobjects
-syn keyword cType matrixobjects
 syn keyword cType AdjointBalancethickness2Analysis
 syn keyword cType AdjointBalancethicknessAnalysis
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27130)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27131)
@@ -283,5 +283,7 @@
 	LoveInnerCoreBoundaryEnum,
 	LoveComplexComputationEnum,
-	LoveIntStepsPerLayerEnum,
+	LoveQuadPrecisionEnum,
+	LoveMinIntegrationStepsEnum,
+	LoveMaxIntegrationdrEnum,
 	LoveKernelsEnum,
 	LoveMu0Enum,
@@ -295,4 +297,7 @@
 	LoveUnderflowTolEnum,
 	LovePostWidderThresholdEnum,
+	LoveDebugEnum,
+	LoveHypergeomNZEnum,
+	LoveHypergeomNAlphaEnum,
 	MassFluxSegmentsEnum,
 	MassFluxSegmentsPresentEnum,
@@ -381,4 +386,6 @@
 	SolidearthSettingsElasticEnum,
 	SolidearthSettingsViscousEnum,
+	SolidearthSettingsSatelliteGraviEnum,
+	SolidearthSettingsDegreeAccuracyEnum,
 	SealevelchangeGeometryDoneEnum,
 	SealevelchangeViscousNumStepsEnum,
@@ -403,4 +410,7 @@
 	LoveTimeFreqEnum,
 	LoveIsTimeEnum,
+	LoveHypergeomZEnum,
+	LoveHypergeomTable1Enum,
+	LoveHypergeomTable2Enum,
 	SealevelchangeGSelfAttractionEnum,
 	SealevelchangeGViscoElasticEnum,
@@ -847,4 +857,5 @@
 	SealevelEnum,
 	SealevelGRDEnum,
+	SatGraviGRDEnum,
 	SealevelBarystaticMaskEnum,
 	SealevelBarystaticIceMaskEnum,
@@ -886,25 +897,19 @@
 	SealevelUNorthEsaEnum,
 	SealevelchangeIndicesEnum,
-	SealevelchangeGEnum,
-	SealevelchangeGUEnum,
-	SealevelchangeGEEnum,
-	SealevelchangeGNEnum,
+	SealevelchangeAlphaIndexEnum,
+	SealevelchangeAzimuthIndexEnum,
 	SealevelchangeGrotEnum,
+	SealevelchangeGSatGravirotEnum,
 	SealevelchangeGUrotEnum,
 	SealevelchangeGNrotEnum,
 	SealevelchangeGErotEnum,
-	SealevelchangeGsubelOceanEnum,
-	SealevelchangeGUsubelOceanEnum,
-	SealevelchangeGEsubelOceanEnum,
-	SealevelchangeGNsubelOceanEnum,
-	SealevelchangeGsubelIceEnum,
-	SealevelchangeGUsubelIceEnum,
-	SealevelchangeGEsubelIceEnum,
-	SealevelchangeGNsubelIceEnum,
-	SealevelchangeGsubelHydroEnum,
-	SealevelchangeGUsubelHydroEnum,
-	SealevelchangeGEsubelHydroEnum,
-	SealevelchangeGNsubelHydroEnum,
+	SealevelchangeAlphaIndexOceanEnum,
+	SealevelchangeAlphaIndexIceEnum,
+	SealevelchangeAlphaIndexHydroEnum,
+	SealevelchangeAzimuthIndexOceanEnum,
+	SealevelchangeAzimuthIndexIceEnum,
+	SealevelchangeAzimuthIndexHydroEnum,
 	SealevelchangeViscousRSLEnum,
+	SealevelchangeViscousSGEnum,
 	SealevelchangeViscousUEnum,
 	SealevelchangeViscousNEnum,
@@ -1410,4 +1415,6 @@
 	LovePMTF1tEnum,
 	LovePMTF2tEnum,
+	LoveYiEnum,
+	LoveRhsEnum,
 	LoveSolutionEnum,
 	MINIEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 27130)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 27131)
@@ -291,5 +291,7 @@
 		case LoveInnerCoreBoundaryEnum : return "LoveInnerCoreBoundary";
 		case LoveComplexComputationEnum : return "LoveComplexComputation";
-		case LoveIntStepsPerLayerEnum : return "LoveIntStepsPerLayer";
+		case LoveQuadPrecisionEnum : return "LoveQuadPrecision";
+		case LoveMinIntegrationStepsEnum : return "LoveMinIntegrationSteps";
+		case LoveMaxIntegrationdrEnum : return "LoveMaxIntegrationdr";
 		case LoveKernelsEnum : return "LoveKernels";
 		case LoveMu0Enum : return "LoveMu0";
@@ -303,4 +305,7 @@
 		case LoveUnderflowTolEnum : return "LoveUnderflowTol";
 		case LovePostWidderThresholdEnum : return "LovePostWidderThreshold";
+		case LoveDebugEnum : return "LoveDebug";
+		case LoveHypergeomNZEnum : return "LoveHypergeomNZ";
+		case LoveHypergeomNAlphaEnum : return "LoveHypergeomNAlpha";
 		case MassFluxSegmentsEnum : return "MassFluxSegments";
 		case MassFluxSegmentsPresentEnum : return "MassFluxSegmentsPresent";
@@ -389,4 +394,6 @@
 		case SolidearthSettingsElasticEnum : return "SolidearthSettingsElastic";
 		case SolidearthSettingsViscousEnum : return "SolidearthSettingsViscous";
+		case SolidearthSettingsSatelliteGraviEnum : return "SolidearthSettingsSatelliteGravi";
+		case SolidearthSettingsDegreeAccuracyEnum : return "SolidearthSettingsDegreeAccuracy";
 		case SealevelchangeGeometryDoneEnum : return "SealevelchangeGeometryDone";
 		case SealevelchangeViscousNumStepsEnum : return "SealevelchangeViscousNumSteps";
@@ -411,4 +418,7 @@
 		case LoveTimeFreqEnum : return "LoveTimeFreq";
 		case LoveIsTimeEnum : return "LoveIsTime";
+		case LoveHypergeomZEnum : return "LoveHypergeomZ";
+		case LoveHypergeomTable1Enum : return "LoveHypergeomTable1";
+		case LoveHypergeomTable2Enum : return "LoveHypergeomTable2";
 		case SealevelchangeGSelfAttractionEnum : return "SealevelchangeGSelfAttraction";
 		case SealevelchangeGViscoElasticEnum : return "SealevelchangeGViscoElastic";
@@ -853,4 +863,5 @@
 		case SealevelEnum : return "Sealevel";
 		case SealevelGRDEnum : return "SealevelGRD";
+		case SatGraviGRDEnum : return "SatGraviGRD";
 		case SealevelBarystaticMaskEnum : return "SealevelBarystaticMask";
 		case SealevelBarystaticIceMaskEnum : return "SealevelBarystaticIceMask";
@@ -892,25 +903,19 @@
 		case SealevelUNorthEsaEnum : return "SealevelUNorthEsa";
 		case SealevelchangeIndicesEnum : return "SealevelchangeIndices";
-		case SealevelchangeGEnum : return "SealevelchangeG";
-		case SealevelchangeGUEnum : return "SealevelchangeGU";
-		case SealevelchangeGEEnum : return "SealevelchangeGE";
-		case SealevelchangeGNEnum : return "SealevelchangeGN";
+		case SealevelchangeAlphaIndexEnum : return "SealevelchangeAlphaIndex";
+		case SealevelchangeAzimuthIndexEnum : return "SealevelchangeAzimuthIndex";
 		case SealevelchangeGrotEnum : return "SealevelchangeGrot";
+		case SealevelchangeGSatGravirotEnum : return "SealevelchangeGSatGravirot";
 		case SealevelchangeGUrotEnum : return "SealevelchangeGUrot";
 		case SealevelchangeGNrotEnum : return "SealevelchangeGNrot";
 		case SealevelchangeGErotEnum : return "SealevelchangeGErot";
-		case SealevelchangeGsubelOceanEnum : return "SealevelchangeGsubelOcean";
-		case SealevelchangeGUsubelOceanEnum : return "SealevelchangeGUsubelOcean";
-		case SealevelchangeGEsubelOceanEnum : return "SealevelchangeGEsubelOcean";
-		case SealevelchangeGNsubelOceanEnum : return "SealevelchangeGNsubelOcean";
-		case SealevelchangeGsubelIceEnum : return "SealevelchangeGsubelIce";
-		case SealevelchangeGUsubelIceEnum : return "SealevelchangeGUsubelIce";
-		case SealevelchangeGEsubelIceEnum : return "SealevelchangeGEsubelIce";
-		case SealevelchangeGNsubelIceEnum : return "SealevelchangeGNsubelIce";
-		case SealevelchangeGsubelHydroEnum : return "SealevelchangeGsubelHydro";
-		case SealevelchangeGUsubelHydroEnum : return "SealevelchangeGUsubelHydro";
-		case SealevelchangeGEsubelHydroEnum : return "SealevelchangeGEsubelHydro";
-		case SealevelchangeGNsubelHydroEnum : return "SealevelchangeGNsubelHydro";
+		case SealevelchangeAlphaIndexOceanEnum : return "SealevelchangeAlphaIndexOcean";
+		case SealevelchangeAlphaIndexIceEnum : return "SealevelchangeAlphaIndexIce";
+		case SealevelchangeAlphaIndexHydroEnum : return "SealevelchangeAlphaIndexHydro";
+		case SealevelchangeAzimuthIndexOceanEnum : return "SealevelchangeAzimuthIndexOcean";
+		case SealevelchangeAzimuthIndexIceEnum : return "SealevelchangeAzimuthIndexIce";
+		case SealevelchangeAzimuthIndexHydroEnum : return "SealevelchangeAzimuthIndexHydro";
 		case SealevelchangeViscousRSLEnum : return "SealevelchangeViscousRSL";
+		case SealevelchangeViscousSGEnum : return "SealevelchangeViscousSG";
 		case SealevelchangeViscousUEnum : return "SealevelchangeViscousU";
 		case SealevelchangeViscousNEnum : return "SealevelchangeViscousN";
@@ -1413,4 +1418,6 @@
 		case LovePMTF1tEnum : return "LovePMTF1t";
 		case LovePMTF2tEnum : return "LovePMTF2t";
+		case LoveYiEnum : return "LoveYi";
+		case LoveRhsEnum : return "LoveRhs";
 		case LoveSolutionEnum : return "LoveSolution";
 		case MINIEnum : return "MINI";
Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 27130)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 27131)
@@ -282,5 +282,7 @@
 syn keyword juliaConstC LoveInnerCoreBoundaryEnum
 syn keyword juliaConstC LoveComplexComputationEnum
-syn keyword juliaConstC LoveIntStepsPerLayerEnum
+syn keyword juliaConstC LoveQuadPrecisionEnum
+syn keyword juliaConstC LoveMinIntegrationStepsEnum
+syn keyword juliaConstC LoveMaxIntegrationdrEnum
 syn keyword juliaConstC LoveKernelsEnum
 syn keyword juliaConstC LoveMu0Enum
@@ -294,4 +296,7 @@
 syn keyword juliaConstC LoveUnderflowTolEnum
 syn keyword juliaConstC LovePostWidderThresholdEnum
+syn keyword juliaConstC LoveDebugEnum
+syn keyword juliaConstC LoveHypergeomNZEnum
+syn keyword juliaConstC LoveHypergeomNAlphaEnum
 syn keyword juliaConstC MassFluxSegmentsEnum
 syn keyword juliaConstC MassFluxSegmentsPresentEnum
@@ -380,4 +385,6 @@
 syn keyword juliaConstC SolidearthSettingsElasticEnum
 syn keyword juliaConstC SolidearthSettingsViscousEnum
+syn keyword juliaConstC SolidearthSettingsSatelliteGraviEnum
+syn keyword juliaConstC SolidearthSettingsDegreeAccuracyEnum
 syn keyword juliaConstC SealevelchangeGeometryDoneEnum
 syn keyword juliaConstC SealevelchangeViscousNumStepsEnum
@@ -402,4 +409,7 @@
 syn keyword juliaConstC LoveTimeFreqEnum
 syn keyword juliaConstC LoveIsTimeEnum
+syn keyword juliaConstC LoveHypergeomZEnum
+syn keyword juliaConstC LoveHypergeomTable1Enum
+syn keyword juliaConstC LoveHypergeomTable2Enum
 syn keyword juliaConstC SealevelchangeGSelfAttractionEnum
 syn keyword juliaConstC SealevelchangeGViscoElasticEnum
@@ -844,4 +854,5 @@
 syn keyword juliaConstC SealevelEnum
 syn keyword juliaConstC SealevelGRDEnum
+syn keyword juliaConstC SatGraviGRDEnum
 syn keyword juliaConstC SealevelBarystaticMaskEnum
 syn keyword juliaConstC SealevelBarystaticIceMaskEnum
@@ -883,25 +894,19 @@
 syn keyword juliaConstC SealevelUNorthEsaEnum
 syn keyword juliaConstC SealevelchangeIndicesEnum
-syn keyword juliaConstC SealevelchangeGEnum
-syn keyword juliaConstC SealevelchangeGUEnum
-syn keyword juliaConstC SealevelchangeGEEnum
-syn keyword juliaConstC SealevelchangeGNEnum
+syn keyword juliaConstC SealevelchangeAlphaIndexEnum
+syn keyword juliaConstC SealevelchangeAzimuthIndexEnum
 syn keyword juliaConstC SealevelchangeGrotEnum
+syn keyword juliaConstC SealevelchangeGSatGravirotEnum
 syn keyword juliaConstC SealevelchangeGUrotEnum
 syn keyword juliaConstC SealevelchangeGNrotEnum
 syn keyword juliaConstC SealevelchangeGErotEnum
-syn keyword juliaConstC SealevelchangeGsubelOceanEnum
-syn keyword juliaConstC SealevelchangeGUsubelOceanEnum
-syn keyword juliaConstC SealevelchangeGEsubelOceanEnum
-syn keyword juliaConstC SealevelchangeGNsubelOceanEnum
-syn keyword juliaConstC SealevelchangeGsubelIceEnum
-syn keyword juliaConstC SealevelchangeGUsubelIceEnum
-syn keyword juliaConstC SealevelchangeGEsubelIceEnum
-syn keyword juliaConstC SealevelchangeGNsubelIceEnum
-syn keyword juliaConstC SealevelchangeGsubelHydroEnum
-syn keyword juliaConstC SealevelchangeGUsubelHydroEnum
-syn keyword juliaConstC SealevelchangeGEsubelHydroEnum
-syn keyword juliaConstC SealevelchangeGNsubelHydroEnum
+syn keyword juliaConstC SealevelchangeAlphaIndexOceanEnum
+syn keyword juliaConstC SealevelchangeAlphaIndexIceEnum
+syn keyword juliaConstC SealevelchangeAlphaIndexHydroEnum
+syn keyword juliaConstC SealevelchangeAzimuthIndexOceanEnum
+syn keyword juliaConstC SealevelchangeAzimuthIndexIceEnum
+syn keyword juliaConstC SealevelchangeAzimuthIndexHydroEnum
 syn keyword juliaConstC SealevelchangeViscousRSLEnum
+syn keyword juliaConstC SealevelchangeViscousSGEnum
 syn keyword juliaConstC SealevelchangeViscousUEnum
 syn keyword juliaConstC SealevelchangeViscousNEnum
@@ -1288,4 +1293,5 @@
 syn keyword juliaConstC DoubleArrayInputEnum
 syn keyword juliaConstC ArrayInputEnum
+syn keyword juliaConstC IntArrayInputEnum
 syn keyword juliaConstC DoubleExternalResultEnum
 syn keyword juliaConstC DoubleMatArrayParamEnum
@@ -1403,4 +1409,6 @@
 syn keyword juliaConstC LovePMTF1tEnum
 syn keyword juliaConstC LovePMTF2tEnum
+syn keyword juliaConstC LoveYiEnum
+syn keyword juliaConstC LoveRhsEnum
 syn keyword juliaConstC LoveSolutionEnum
 syn keyword juliaConstC MINIEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 27130)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 27131)
@@ -297,5 +297,7 @@
 	      else if (strcmp(name,"LoveInnerCoreBoundary")==0) return LoveInnerCoreBoundaryEnum;
 	      else if (strcmp(name,"LoveComplexComputation")==0) return LoveComplexComputationEnum;
-	      else if (strcmp(name,"LoveIntStepsPerLayer")==0) return LoveIntStepsPerLayerEnum;
+	      else if (strcmp(name,"LoveQuadPrecision")==0) return LoveQuadPrecisionEnum;
+	      else if (strcmp(name,"LoveMinIntegrationSteps")==0) return LoveMinIntegrationStepsEnum;
+	      else if (strcmp(name,"LoveMaxIntegrationdr")==0) return LoveMaxIntegrationdrEnum;
 	      else if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum;
 	      else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum;
@@ -309,4 +311,7 @@
 	      else if (strcmp(name,"LoveUnderflowTol")==0) return LoveUnderflowTolEnum;
 	      else if (strcmp(name,"LovePostWidderThreshold")==0) return LovePostWidderThresholdEnum;
+	      else if (strcmp(name,"LoveDebug")==0) return LoveDebugEnum;
+	      else if (strcmp(name,"LoveHypergeomNZ")==0) return LoveHypergeomNZEnum;
+	      else if (strcmp(name,"LoveHypergeomNAlpha")==0) return LoveHypergeomNAlphaEnum;
 	      else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
 	      else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
@@ -378,13 +383,13 @@
 	      else if (strcmp(name,"Modelname")==0) return ModelnameEnum;
 	      else if (strcmp(name,"SamplingAlpha")==0) return SamplingAlphaEnum;
-	      else if (strcmp(name,"SamplingNumRequestedOutputs")==0) return SamplingNumRequestedOutputsEnum;
+         else stage=4;
+   }
+   if(stage==4){
+	      if (strcmp(name,"SamplingNumRequestedOutputs")==0) return SamplingNumRequestedOutputsEnum;
 	      else if (strcmp(name,"SamplingRequestedOutputs")==0) return SamplingRequestedOutputsEnum;
 	      else if (strcmp(name,"SamplingRobin")==0) return SamplingRobinEnum;
 	      else if (strcmp(name,"SamplingSeed")==0) return SamplingSeedEnum;
 	      else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
-         else stage=4;
-   }
-   if(stage==4){
-	      if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum;
+	      else if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum;
 	      else if (strcmp(name,"SolidearthPartitionHydro")==0) return SolidearthPartitionHydroEnum;
 	      else if (strcmp(name,"SolidearthPartitionOcean")==0) return SolidearthPartitionOceanEnum;
@@ -398,4 +403,6 @@
 	      else if (strcmp(name,"SolidearthSettingsElastic")==0) return SolidearthSettingsElasticEnum;
 	      else if (strcmp(name,"SolidearthSettingsViscous")==0) return SolidearthSettingsViscousEnum;
+	      else if (strcmp(name,"SolidearthSettingsSatelliteGravi")==0) return SolidearthSettingsSatelliteGraviEnum;
+	      else if (strcmp(name,"SolidearthSettingsDegreeAccuracy")==0) return SolidearthSettingsDegreeAccuracyEnum;
 	      else if (strcmp(name,"SealevelchangeGeometryDone")==0) return SealevelchangeGeometryDoneEnum;
 	      else if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum;
@@ -420,4 +427,7 @@
 	      else if (strcmp(name,"LoveTimeFreq")==0) return LoveTimeFreqEnum;
 	      else if (strcmp(name,"LoveIsTime")==0) return LoveIsTimeEnum;
+	      else if (strcmp(name,"LoveHypergeomZ")==0) return LoveHypergeomZEnum;
+	      else if (strcmp(name,"LoveHypergeomTable1")==0) return LoveHypergeomTable1Enum;
+	      else if (strcmp(name,"LoveHypergeomTable2")==0) return LoveHypergeomTable2Enum;
 	      else if (strcmp(name,"SealevelchangeGSelfAttraction")==0) return SealevelchangeGSelfAttractionEnum;
 	      else if (strcmp(name,"SealevelchangeGViscoElastic")==0) return SealevelchangeGViscoElasticEnum;
@@ -496,5 +506,8 @@
 	      else if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum;
 	      else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum;
-	      else if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum;
+         else stage=5;
+   }
+   if(stage==5){
+	      if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum;
 	      else if (strcmp(name,"SmbIsmungsm")==0) return SmbIsmungsmEnum;
 	      else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum;
@@ -506,8 +519,5 @@
 	      else if (strcmp(name,"SmbK")==0) return SmbKEnum;
 	      else if (strcmp(name,"SmbLapseRates")==0) return SmbLapseRatesEnum;
-         else stage=5;
-   }
-   if(stage==5){
-	      if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum;
+	      else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum;
 	      else if (strcmp(name,"SmbNumElevationBins")==0) return SmbNumElevationBinsEnum;
 	      else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
@@ -619,5 +629,8 @@
 	      else if (strcmp(name,"Approximation")==0) return ApproximationEnum;
 	      else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
-	      else if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum;
+         else stage=6;
+   }
+   if(stage==6){
+	      if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum;
 	      else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum;
 	      else if (strcmp(name,"BalancethicknessSpcthickness")==0) return BalancethicknessSpcthicknessEnum;
@@ -629,8 +642,5 @@
 	      else if (strcmp(name,"BasalforcingsFloatingiceMeltingRate")==0) return BasalforcingsFloatingiceMeltingRateEnum;
 	      else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum;
-         else stage=6;
-   }
-   if(stage==6){
-	      if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum;
+	      else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum;
 	      else if (strcmp(name,"BasalforcingsLinearBasinId")==0) return BasalforcingsLinearBasinIdEnum;
 	      else if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum;
@@ -742,5 +752,8 @@
 	      else if (strcmp(name,"EsaYmotion")==0) return EsaYmotionEnum;
 	      else if (strcmp(name,"EtaDiff")==0) return EtaDiffEnum;
-	      else if (strcmp(name,"FlowequationBorderFS")==0) return FlowequationBorderFSEnum;
+         else stage=7;
+   }
+   if(stage==7){
+	      if (strcmp(name,"FlowequationBorderFS")==0) return FlowequationBorderFSEnum;
 	      else if (strcmp(name,"FrictionAs")==0) return FrictionAsEnum;
 	      else if (strcmp(name,"FrictionC")==0) return FrictionCEnum;
@@ -752,8 +765,5 @@
 	      else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
 	      else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;
-         else stage=7;
-   }
-   if(stage==7){
-	      if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
+	      else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
 	      else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
 	      else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
@@ -865,5 +875,8 @@
 	      else if (strcmp(name,"SampleOld")==0) return SampleOldEnum;
 	      else if (strcmp(name,"SampleNoise")==0) return SampleNoiseEnum;
-	      else if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum;
+         else stage=8;
+   }
+   if(stage==8){
+	      if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum;
 	      else if (strcmp(name,"SamplingKappa")==0) return SamplingKappaEnum;
 	      else if (strcmp(name,"SamplingPhi")==0) return SamplingPhiEnum;
@@ -871,12 +884,10 @@
 	      else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
 	      else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum;
+	      else if (strcmp(name,"SatGraviGRD")==0) return SatGraviGRDEnum;
 	      else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum;
 	      else if (strcmp(name,"SealevelBarystaticIceMask")==0) return SealevelBarystaticIceMaskEnum;
 	      else if (strcmp(name,"SealevelBarystaticIceWeights")==0) return SealevelBarystaticIceWeightsEnum;
 	      else if (strcmp(name,"SealevelBarystaticIceArea")==0) return SealevelBarystaticIceAreaEnum;
-         else stage=8;
-   }
-   if(stage==8){
-	      if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum;
+	      else if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum;
 	      else if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum;
 	      else if (strcmp(name,"SealevelBarystaticIceLoad")==0) return SealevelBarystaticIceLoadEnum;
@@ -913,25 +924,19 @@
 	      else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
 	      else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
-	      else if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum;
-	      else if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum;
-	      else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum;
-	      else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum;
+	      else if (strcmp(name,"SealevelchangeAlphaIndex")==0) return SealevelchangeAlphaIndexEnum;
+	      else if (strcmp(name,"SealevelchangeAzimuthIndex")==0) return SealevelchangeAzimuthIndexEnum;
 	      else if (strcmp(name,"SealevelchangeGrot")==0) return SealevelchangeGrotEnum;
+	      else if (strcmp(name,"SealevelchangeGSatGravirot")==0) return SealevelchangeGSatGravirotEnum;
 	      else if (strcmp(name,"SealevelchangeGUrot")==0) return SealevelchangeGUrotEnum;
 	      else if (strcmp(name,"SealevelchangeGNrot")==0) return SealevelchangeGNrotEnum;
 	      else if (strcmp(name,"SealevelchangeGErot")==0) return SealevelchangeGErotEnum;
-	      else if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum;
-	      else if (strcmp(name,"SealevelchangeGUsubelOcean")==0) return SealevelchangeGUsubelOceanEnum;
-	      else if (strcmp(name,"SealevelchangeGEsubelOcean")==0) return SealevelchangeGEsubelOceanEnum;
-	      else if (strcmp(name,"SealevelchangeGNsubelOcean")==0) return SealevelchangeGNsubelOceanEnum;
-	      else if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum;
-	      else if (strcmp(name,"SealevelchangeGUsubelIce")==0) return SealevelchangeGUsubelIceEnum;
-	      else if (strcmp(name,"SealevelchangeGEsubelIce")==0) return SealevelchangeGEsubelIceEnum;
-	      else if (strcmp(name,"SealevelchangeGNsubelIce")==0) return SealevelchangeGNsubelIceEnum;
-	      else if (strcmp(name,"SealevelchangeGsubelHydro")==0) return SealevelchangeGsubelHydroEnum;
-	      else if (strcmp(name,"SealevelchangeGUsubelHydro")==0) return SealevelchangeGUsubelHydroEnum;
-	      else if (strcmp(name,"SealevelchangeGEsubelHydro")==0) return SealevelchangeGEsubelHydroEnum;
-	      else if (strcmp(name,"SealevelchangeGNsubelHydro")==0) return SealevelchangeGNsubelHydroEnum;
+	      else if (strcmp(name,"SealevelchangeAlphaIndexOcean")==0) return SealevelchangeAlphaIndexOceanEnum;
+	      else if (strcmp(name,"SealevelchangeAlphaIndexIce")==0) return SealevelchangeAlphaIndexIceEnum;
+	      else if (strcmp(name,"SealevelchangeAlphaIndexHydro")==0) return SealevelchangeAlphaIndexHydroEnum;
+	      else if (strcmp(name,"SealevelchangeAzimuthIndexOcean")==0) return SealevelchangeAzimuthIndexOceanEnum;
+	      else if (strcmp(name,"SealevelchangeAzimuthIndexIce")==0) return SealevelchangeAzimuthIndexIceEnum;
+	      else if (strcmp(name,"SealevelchangeAzimuthIndexHydro")==0) return SealevelchangeAzimuthIndexHydroEnum;
 	      else if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum;
+	      else if (strcmp(name,"SealevelchangeViscousSG")==0) return SealevelchangeViscousSGEnum;
 	      else if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum;
 	      else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum;
@@ -993,13 +998,13 @@
 	      else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
 	      else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
-	      else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
+         else stage=9;
+   }
+   if(stage==9){
+	      if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
 	      else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum;
 	      else if (strcmp(name,"SmbGdn")==0) return SmbGdnEnum;
 	      else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum;
 	      else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum;
-         else stage=9;
-   }
-   if(stage==9){
-	      if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
+	      else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
 	      else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
 	      else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
@@ -1116,13 +1121,13 @@
 	      else if (strcmp(name,"VxBase")==0) return VxBaseEnum;
 	      else if (strcmp(name,"Vx")==0) return VxEnum;
-	      else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
+         else stage=10;
+   }
+   if(stage==10){
+	      if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
 	      else if (strcmp(name,"VxObs")==0) return VxObsEnum;
 	      else if (strcmp(name,"VxShear")==0) return VxShearEnum;
 	      else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum;
 	      else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
-         else stage=10;
-   }
-   if(stage==10){
-	      if (strcmp(name,"VyBase")==0) return VyBaseEnum;
+	      else if (strcmp(name,"VyBase")==0) return VyBaseEnum;
 	      else if (strcmp(name,"Vy")==0) return VyEnum;
 	      else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
@@ -1239,13 +1244,13 @@
 	      else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
 	      else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
-	      else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
+         else stage=11;
+   }
+   if(stage==11){
+	      if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
 	      else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum;
 	      else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum;
 	      else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum;
 	      else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
-         else stage=11;
-   }
-   if(stage==11){
-	      if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
+	      else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
 	      else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
 	      else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum;
@@ -1362,13 +1367,13 @@
 	      else if (strcmp(name,"FrontalForcingsRignotAutoregression")==0) return FrontalForcingsRignotAutoregressionEnum;
 	      else if (strcmp(name,"Fset")==0) return FsetEnum;
-	      else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum;
+         else stage=12;
+   }
+   if(stage==12){
+	      if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum;
 	      else if (strcmp(name,"GLheightadvectionAnalysis")==0) return GLheightadvectionAnalysisEnum;
 	      else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
 	      else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
 	      else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
-         else stage=12;
-   }
-   if(stage==12){
-	      if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
+	      else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
 	      else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
 	      else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
@@ -1446,4 +1451,6 @@
 	      else if (strcmp(name,"LovePMTF1t")==0) return LovePMTF1tEnum;
 	      else if (strcmp(name,"LovePMTF2t")==0) return LovePMTF2tEnum;
+	      else if (strcmp(name,"LoveYi")==0) return LoveYiEnum;
+	      else if (strcmp(name,"LoveRhs")==0) return LoveRhsEnum;
 	      else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum;
 	      else if (strcmp(name,"MINI")==0) return MINIEnum;
@@ -1483,5 +1490,8 @@
 	      else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum;
 	      else if (strcmp(name,"Moulin")==0) return MoulinEnum;
-	      else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
+         else stage=13;
+   }
+   if(stage==13){
+	      if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
 	      else if (strcmp(name,"Mpi")==0) return MpiEnum;
 	      else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
@@ -1490,8 +1500,5 @@
 	      else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum;
 	      else if (strcmp(name,"Nodal")==0) return NodalEnum;
-         else stage=13;
-   }
-   if(stage==13){
-	      if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
+	      else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
 	      else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
 	      else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
@@ -1606,5 +1613,8 @@
 	      else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
 	      else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;
-	      else if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
+         else stage=14;
+   }
+   if(stage==14){
+	      if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
 	      else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum;
 	      else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum;
@@ -1613,8 +1623,5 @@
 	      else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
 	      else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum;
-         else stage=14;
-   }
-   if(stage==14){
-	      if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
+	      else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
 	      else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
 	      else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/solidearthsettings.m	(revision 27130)
+++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m	(revision 27131)
@@ -20,7 +20,7 @@
 		compute_bp_grd         = 0; %will GRD patterns for bottom pressures be computed? 
 		degacc                 = 0; %degree increment for resolution of Green tables.
-		timeacc                = 0; %time step accuracy required to compute Green tables
+		timeacc                = 1; %time step accuracy required to compute Green tables
 		horiz                  = 0; %compute horizontal deformation
-		grdmodel               = 0; %grd model (0 by default, 1 for (visco-)elastic, 2 for Ivins)
+		grdmodel               = 1; %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
@@ -89,27 +89,31 @@
 
 			%no grd model by default:
-			self.grdmodel=0;
+			self.grdmodel=1;
 
 		end % }}}
 		function disp(self) % {{{
 			disp(sprintf('   solidearth settings:'));
-
-			fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)');
-			fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, 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,'sealevelloading','enables surface loading from sea-level change (default: 1)');
+			disp(sprintf('      core:'));
 			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,'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,'runfrequency','How many time steps we let masstransport core accumulate changes before each run of the sealevelchange core (default: 1, i.e run slc every time step)');
+			disp(sprintf('      computational flags:'));
 			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,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)');
+			fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore. Used only for grdmodel=2 only');
+			disp(sprintf('      resolution:'));
+			fielddisplay(self,'degacc','spatial 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, 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');
+			disp(sprintf('      sea-level equation:'));
+			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,'sealevelloading','enables surface loading from sea-level change (default: 1)');
+			fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
+			fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)');
+			fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)');
+			fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'); 
+
 		end % }}}
 		function md = checkconsistency(self,md,solution,analyses) % {{{
@@ -133,4 +137,7 @@
 			if self.viscous==1 & self.elastic==0,
 				error('solidearthsettings checkconsistency error message: need elastic on if viscous flag is set');
+			end
+			if self.rotation==1 & self.elastic==0,
+				error('solidearthsettings checkconsistency error message: need elastic on if rotation flag is set');
 			end
 
@@ -164,4 +171,5 @@
 			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','rotation','name','md.solidearth.settings.satellitegravity','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');
@@ -187,5 +195,5 @@
 			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.grdocean'],self.grdocean);
 			writejsdouble(fid,[modelname '.solidearth.settings.ocean_area_scaling'],self.ocean_area_scaling);
 			writejsdouble(fid,[modelname '.solidearth.settings.run_frequency'],self.run_frequency);