Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp	(revision 25654)
+++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp	(revision 25655)
@@ -175,5 +175,6 @@
 	IssmDouble  planetradius=0;
 	IssmDouble  planetarea=0;
-	bool         elastic=false;
+	bool		rigid=false;
+	bool		elastic=false;
 
 	int     numoutputs;
@@ -231,154 +232,216 @@
 	} /*}}}*/
 	/*Deal with elasticity {{{*/
+	iomodel->FetchData(&rigid,"md.solidearth.settings.rigid");
 	iomodel->FetchData(&elastic,"md.solidearth.settings.elastic");
-// if(elastic){
-		/*love numbers: */
-		iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
-		iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
-		iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
-		iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th");
-		iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
-		iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
-		parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
-
-		parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
-		parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1));
-		parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1));
-		parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1));
-		parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1));
-		parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
-
-		/*compute elastic green function for a range of angles*/
-		iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
-		M=reCast<int,IssmDouble>(180./degacc+1.);
-
-		// AD performance is sensitive to calls to ensurecontiguous.
-		// // Providing "t" will cause ensurecontiguous to be called.
-		#ifdef _HAVE_AD_
-		G_rigid=xNew<IssmDouble>(M,"t");
-		G_elastic=xNew<IssmDouble>(M,"t");
-		U_elastic=xNew<IssmDouble>(M,"t");
-		H_elastic=xNew<IssmDouble>(M,"t");
-		#else
-		G_rigid=xNew<IssmDouble>(M);
-		G_elastic=xNew<IssmDouble>(M);
-		U_elastic=xNew<IssmDouble>(M);
-		H_elastic=xNew<IssmDouble>(M);
-		#endif
-
-		/*compute combined legendre + love number (elastic green function:*/
-		m=DetermineLocalSize(M,IssmComm::GetComm());
-		GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
-		// AD performance is sensitive to calls to ensurecontiguous.
-		// // Providing "t" will cause ensurecontiguous to be called.
-		#ifdef _HAVE_AD_
-		G_elastic_local=xNew<IssmDouble>(m,"t");
-		G_rigid_local=xNew<IssmDouble>(m,"t");
-		U_elastic_local=xNew<IssmDouble>(m,"t");
-		H_elastic_local=xNew<IssmDouble>(m,"t");
-		#else
-		G_elastic_local=xNew<IssmDouble>(m);
-		G_rigid_local=xNew<IssmDouble>(m);
-		U_elastic_local=xNew<IssmDouble>(m);
-		H_elastic_local=xNew<IssmDouble>(m);
-		#endif
-
-		for(int i=lower_row;i<upper_row;i++){
-			IssmDouble alpha,x;
-			alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
-
-			G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
-			G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
-			U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
-			H_elastic_local[i-lower_row]= 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<nl;n++) {
-				IssmDouble deltalove_G;
-				IssmDouble deltalove_U;
-
-				deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
-				deltalove_U = (love_h[n]-love_h[nl-1]);
-
-				/*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
-				if(n==0){
-					Pn=1; 
-					Pn_p=0; 
+	if (rigid) {
+		if (elastic) {
+			/*love numbers: */
+			iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
+			iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
+			iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
+			iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th");
+			iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
+			iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
+			parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
+
+			parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
+			parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1));
+			parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1));
+			parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1));
+			parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1));
+			parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
+
+			/*compute elastic green function for a range of angles*/
+			iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
+			M=reCast<int,IssmDouble>(180./degacc+1.);
+
+			// AD performance is sensitive to calls to ensurecontiguous.
+			// // Providing "t" will cause ensurecontiguous to be called.
+			#ifdef _HAVE_AD_
+			G_rigid=xNew<IssmDouble>(M,"t");
+			G_elastic=xNew<IssmDouble>(M,"t");
+			U_elastic=xNew<IssmDouble>(M,"t");
+			H_elastic=xNew<IssmDouble>(M,"t");
+			#else
+			G_rigid=xNew<IssmDouble>(M);
+			G_elastic=xNew<IssmDouble>(M);
+			U_elastic=xNew<IssmDouble>(M);
+			H_elastic=xNew<IssmDouble>(M);
+			#endif
+
+			/*compute combined legendre + love number (elastic green function:*/
+			m=DetermineLocalSize(M,IssmComm::GetComm());
+			GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
+			// AD performance is sensitive to calls to ensurecontiguous.
+			// // Providing "t" will cause ensurecontiguous to be called.
+			#ifdef _HAVE_AD_
+			G_elastic_local=xNew<IssmDouble>(m,"t");
+			G_rigid_local=xNew<IssmDouble>(m,"t");
+			U_elastic_local=xNew<IssmDouble>(m,"t");
+			H_elastic_local=xNew<IssmDouble>(m,"t");
+			#else
+			G_elastic_local=xNew<IssmDouble>(m);
+			G_rigid_local=xNew<IssmDouble>(m);
+			U_elastic_local=xNew<IssmDouble>(m);
+			H_elastic_local=xNew<IssmDouble>(m);
+			#endif
+
+			for(int i=lower_row;i<upper_row;i++){
+				IssmDouble alpha,x;
+				alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
+
+				G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
+				G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
+				U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
+				H_elastic_local[i-lower_row]= 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<nl;n++) {
+					IssmDouble deltalove_G;
+					IssmDouble deltalove_U;
+
+					deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
+					deltalove_U = (love_h[n]-love_h[nl-1]);
+
+					/*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;
+
+					G_elastic_local[i-lower_row] += deltalove_G*Pn;		// gravitational potential 
+					U_elastic_local[i-lower_row] += deltalove_U*Pn;		// vertical (up) displacement 
+					H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p;		// horizontal displacements 
 				}
-				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;
-
-				G_elastic_local[i-lower_row] += deltalove_G*Pn;		// gravitational potential 
-				U_elastic_local[i-lower_row] += deltalove_U*Pn;		// vertical (up) displacement 
-				H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p;		// horizontal displacements 
 			}
-		}
-
-		/*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
-		int* recvcounts=xNew<int>(IssmComm::GetSize());
-		int* displs=xNew<int>(IssmComm::GetSize());
-
-		//recvcounts:
-		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_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-		ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-		ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-		ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-		
-		/*free resources: */
-		xDelete<int>(recvcounts);
-		xDelete<int>(displs);
-
-		/*}}}*/
-
-		/*Avoid singularity at 0: */
-		G_rigid[0]=G_rigid[1];
-		parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
-		G_elastic[0]=G_elastic[1];
-		parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
-		U_elastic[0]=U_elastic[1];
-		parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
-		H_elastic[0]=H_elastic[1];
-		parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
-
-		/*free resources: */
-		xDelete<IssmDouble>(love_h);
-		xDelete<IssmDouble>(love_k);
-		xDelete<IssmDouble>(love_l);
-		xDelete<IssmDouble>(love_th);
-		xDelete<IssmDouble>(love_tk);
-		xDelete<IssmDouble>(love_tl);
-		xDelete<IssmDouble>(G_rigid);
-		xDelete<IssmDouble>(G_rigid_local);
-		xDelete<IssmDouble>(G_elastic);
-		xDelete<IssmDouble>(G_elastic_local);
-		xDelete<IssmDouble>(U_elastic);
-		xDelete<IssmDouble>(U_elastic_local);
-		xDelete<IssmDouble>(H_elastic);
-		xDelete<IssmDouble>(H_elastic_local);
-	// }
+
+			/*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
+			int* recvcounts=xNew<int>(IssmComm::GetSize());
+			int* displs=xNew<int>(IssmComm::GetSize());
+
+			//recvcounts:
+			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_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+			ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+			ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+			ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+			
+			/*free resources: */
+			xDelete<int>(recvcounts);
+			xDelete<int>(displs);
+
+			/*}}}*/
+
+			/*Avoid singularity at 0: */
+			G_rigid[0]=G_rigid[1];
+			parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
+			G_elastic[0]=G_elastic[1];
+			parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
+			U_elastic[0]=U_elastic[1];
+			parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
+			H_elastic[0]=H_elastic[1];
+			parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
+
+			/*free resources: */
+			xDelete<IssmDouble>(love_h);
+			xDelete<IssmDouble>(love_k);
+			xDelete<IssmDouble>(love_l);
+			xDelete<IssmDouble>(love_th);
+			xDelete<IssmDouble>(love_tk);
+			xDelete<IssmDouble>(love_tl);
+			xDelete<IssmDouble>(G_rigid);
+			xDelete<IssmDouble>(G_rigid_local);
+			xDelete<IssmDouble>(G_elastic);
+			xDelete<IssmDouble>(G_elastic_local);
+			xDelete<IssmDouble>(U_elastic);
+			xDelete<IssmDouble>(U_elastic_local);
+			xDelete<IssmDouble>(H_elastic);
+			xDelete<IssmDouble>(H_elastic_local);
+		} else { /*elastic==false*/
+			/*compute elastic green function for a range of angles*/
+			iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
+			M=reCast<int,IssmDouble>(180./degacc+1.);
+			
+			// AD performance is sensitive to calls to ensurecontiguous.
+			// // Providing "t" will cause ensurecontiguous to be called.
+			#ifdef _HAVE_AD_
+			G_rigid=xNew<IssmDouble>(M,"t");
+			#else
+			G_rigid=xNew<IssmDouble>(M);
+			#endif
+
+			/*compute combined legendre + love number (elastic green function:*/
+			m=DetermineLocalSize(M,IssmComm::GetComm());
+			GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
+			// AD performance is sensitive to calls to ensurecontiguous.
+			// // Providing "t" will cause ensurecontiguous to be called.
+			#ifdef _HAVE_AD_
+			G_rigid_local=xNew<IssmDouble>(m,"t");
+			#else
+			G_rigid_local=xNew<IssmDouble>(m);
+			#endif
+
+			for(int i=lower_row;i<upper_row;i++){
+				IssmDouble alpha,x;
+				alpha=reCast<IssmDouble>(i)*degacc*PI/180.0;
+				G_rigid_local[i-lower_row]=.5/sin(alpha/2.0);
+			}
+
+			/*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
+			int* recvcounts=xNew<int>(IssmComm::GetSize());
+			int* displs=xNew<int>(IssmComm::GetSize());
+
+			//recvcounts:
+			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_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+			
+			/*free resources: */
+			xDelete<int>(recvcounts);
+			xDelete<int>(displs);
+
+			/*}}}*/
+
+			/*Avoid singularity at 0: */
+			G_rigid[0]=G_rigid[1];
+			parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
+
+			/*free resources: */
+			xDelete<IssmDouble>(G_rigid);
+			xDelete<IssmDouble>(G_rigid_local);
+		}
+	} else { /*rigid==false*/
+		if (elastic) {
+			_error_("Must set md.solidearth.settings.rigid=1 when setting md.solidearth.settings.elastic=1");
+		}
+	}
 	
-	/*Indicate we have not yet ran the Geometry Core module: */
+	/*Indicate we have not yet run the Geometry Core module: */
 	parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false));
 
-/*}}}*/
 	/*Transitions:{{{ */
 	iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.solidearth.transitions");
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25654)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25655)
@@ -5475,113 +5475,160 @@
 
 	/*recover precomputed green function kernels:*/
-	DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
-	parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
-
-	parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
-	parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
-
-	parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
-	parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
-
-	parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
-	parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
-
-	/*allocate indices:*/
-	indices=xNew<IssmDouble>(gsize);
-	G=xNewZeroInit<IssmDouble>(gsize);
-	GU=xNewZeroInit<IssmDouble>(gsize);
-	if(horiz){
-		GN=xNewZeroInit<IssmDouble>(gsize);
-		GE=xNewZeroInit<IssmDouble>(gsize);
-	}
-
-	
-	/* Where is the centroid of this element:*/
-
-	/*retrieve coordinates: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*figure out gravity center of our element (Cartesian): */
-	x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
-	y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
-	z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
-
-	/*compute gravity center in lat,long: */
-	late= asin(z_element/planetradius);
-	longe = atan2(y_element,x_element);
-
-	constant=3/rho_earth*area/planetarea;
-
-	for(int i=0;i<gsize;i++){
-		IssmDouble alpha;
-		IssmDouble delPhi,delLambda;
-
-		/*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
-		lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
-		delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)delLambda=2*PI-delLambda;
-		alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
-		indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
-		index=reCast<int,IssmDouble>(indices[i]);
-
-		/*Rigid earth gravitational perturbation: */
-		if(computerigid){
-			G[i]+=G_rigid_precomputed[index];
-		}
+	if(computerigid){
+		DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
+		parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
+
+		/*allocate indices:*/
+		indices=xNew<IssmDouble>(gsize);
+		G=xNewZeroInit<IssmDouble>(gsize);
+
 		if(computeelastic){
-			G[i]+=G_elastic_precomputed[index];
-		}
-		G[i]=G[i]*constant;
-
-		/*Elastic components:*/
-		GU[i] =  constant * U_elastic_precomputed[index];
-		if(horiz){
-			/*Compute azimuths, both north and east components: */
-			x = xx[i]; y = yy[i]; z = zz[i];
-			if(latitude[i]==90){
-				x=1e-12; y=1e-12;
-			}
-			if(latitude[i]==-90){
-				x=1e-12; y=1e-12;
-			}
-			dx = x_element-x; dy = y_element-y; dz = z_element-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);
-
-			GN[i] = constant*H_elastic_precomputed[index]*N_azim;
-			GE[i] = constant*H_elastic_precomputed[index]*E_azim;
-		}
-	}
-
-	/*Add in inputs:*/
-    this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
-    this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
-    this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
-    if(horiz){
-		this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize);
-		this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
-	}
-
-	/* NOTE: May need to remove these two lines (from merge) */
-	this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
-	this->AddInput(SealevelAreaEnum,&area,P0Enum);
-
-	/*Free allocations:*/
-	#ifdef _HAVE_RESTRICT_
-	delete indices;
-	delete G;
-	delete GU;
-	if(horiz){
-		delete GN;
-		delete GE;
-	}
-	#else
-	xDelete(indices);
-	xDelete(G);
-	xDelete(GU);
-	if(horiz){
-		xDelete(GN);
-		xDelete(GE);
-	}
-	#endif
+			parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
+			parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
+
+			parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
+			parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
+
+			parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
+			parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
+
+			/*allocate indices:*/
+			GU=xNewZeroInit<IssmDouble>(gsize);
+			if(horiz){
+				GN=xNewZeroInit<IssmDouble>(gsize);
+				GE=xNewZeroInit<IssmDouble>(gsize);
+			}
+
+			/* Where is the centroid of this element:*/
+
+			/*retrieve coordinates: */
+			::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+			/*figure out gravity center of our element (Cartesian): */
+			x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
+			y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
+			z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
+
+			/*compute gravity center in lat,long: */
+			late= asin(z_element/planetradius);
+			longe = atan2(y_element,x_element);
+
+			constant=3/rho_earth*area/planetarea;
+
+			for(int i=0;i<gsize;i++){
+				IssmDouble alpha;
+				IssmDouble delPhi,delLambda;
+
+				/*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
+				lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
+				delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)delLambda=2*PI-delLambda;
+				alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
+				indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
+				index=reCast<int,IssmDouble>(indices[i]);
+
+				/*Rigid earth gravitational perturbation: */
+				if(computerigid){
+					G[i]+=G_rigid_precomputed[index];
+				}
+				if(computeelastic){
+					G[i]+=G_elastic_precomputed[index];
+				}
+				G[i]=G[i]*constant;
+
+				/*Elastic components:*/
+				GU[i] =  constant * U_elastic_precomputed[index];
+				if(horiz){
+					/*Compute azimuths, both north and east components: */
+					x = xx[i]; y = yy[i]; z = zz[i];
+					if(latitude[i]==90){
+						x=1e-12; y=1e-12;
+					}
+					if(latitude[i]==-90){
+						x=1e-12; y=1e-12;
+					}
+					dx = x_element-x; dy = y_element-y; dz = z_element-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);
+
+					GN[i] = constant*H_elastic_precomputed[index]*N_azim;
+					GE[i] = constant*H_elastic_precomputed[index]*E_azim;
+				}
+			}
+
+			/*Add in inputs:*/
+		    this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
+		    this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
+		    this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
+		    if(horiz){
+				this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize);
+				this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
+			}
+			this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
+			this->AddInput(SealevelAreaEnum,&area,P0Enum);
+
+			/*Free allocations:*/
+			#ifdef _HAVE_RESTRICT_
+			delete GU;
+			if(horiz){
+				delete GN;
+				delete GE;
+			}
+			#else
+			xDelete(GU);
+			if(horiz){
+				xDelete(GN);
+				xDelete(GE);
+			}
+			#endif
+		} else { /*computeelastic==false*/
+			/* Where is the centroid of this element:*/
+
+			/*retrieve coordinates: */
+			::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+			/*figure out gravity center of our element (Cartesian): */
+			x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
+			y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
+			z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
+
+			/*compute gravity center in lat,long: */
+			late= asin(z_element/planetradius);
+			longe = atan2(y_element,x_element);
+
+			constant=3/rho_earth*area/planetarea;
+			for(int i=0;i<gsize;i++){
+				IssmDouble alpha;
+				IssmDouble delPhi,delLambda;
+
+				/*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
+				lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
+				delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)delLambda=2*PI-delLambda;
+				alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
+				indices[i]=alpha/PI*reCast<IssmDouble,int>(M-1);
+				index=reCast<int,IssmDouble>(indices[i]);
+
+				/*Rigid earth gravitational perturbation: */
+				G[i]+=G_rigid_precomputed[index];
+				G[i]=G[i]*constant;
+			}
+
+			/*Add in inputs:*/
+		    this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
+		    this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
+			this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
+			this->AddInput(SealevelAreaEnum,&area,P0Enum);
+
+			/*Free allocations:*/
+			#ifdef _HAVE_RESTRICT_
+			delete indices;
+			delete G;
+			#else
+			xDelete(indices);
+			xDelete(G);
+			#endif
+		}
+	} else { /*computerigid==false*/
+		// do nothing
+	}
 
 	return;
@@ -5646,6 +5693,4 @@
 
 	/*retrieve precomputed G:*/
-
-	/* NOTE: May need to remove condition here (from merge) */
 	if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
 
@@ -5752,6 +5797,4 @@
 
 	/*retrieve precomputed G:*/
-
-	/* NOTE: May need to remove condition here (from merge) */
 	if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);