Index: /issm/branches/trunk-larour-SLPS2020/src/c/analyses/SealevelriseAnalysis.cpp
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/c/analyses/SealevelriseAnalysis.cpp	(revision 25565)
+++ /issm/branches/trunk-larour-SLPS2020/src/c/analyses/SealevelriseAnalysis.cpp	(revision 25566)
@@ -175,4 +175,5 @@
 	IssmDouble  planetradius=0;
 	IssmDouble  planetarea=0;
+	bool         elastic=false;
 
 	int     numoutputs;
@@ -230,145 +231,148 @@
 	} /*}}}*/
 	/*Deal with elasticity {{{*/
-	/*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; 
+	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; 
+				}
+				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 ressources: */
+		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 ressources: */
+		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 ressources: */
-	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 ressources: */
-	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);
 	
 	/*Indicate we have not yet ran the Geometry Core module: */
Index: /issm/branches/trunk-larour-SLPS2020/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/c/classes/Elements/Tria.cpp	(revision 25565)
+++ /issm/branches/trunk-larour-SLPS2020/src/c/classes/Elements/Tria.cpp	(revision 25566)
@@ -5539,4 +5539,10 @@
 	this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
 
+	/*compute area and add to inputs:*/
+	area=GetAreaSpherical();
+	this->inputs2->SetDoubleInput(AreaEnum,this->lid,area);
+	this->AddInput2(SealevelAreaEnum,&area,P0Enum);
+	if(!computerigid & !computeelastic)return;
+
 	/*recover precomputed green function kernels:*/
 	DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
@@ -5561,8 +5567,5 @@
 	}
 
-	/*compute area:*/
-	area=GetAreaSpherical();
-
-
+	
 	/* Where is the centroid of this element:*/
 
@@ -5629,7 +5632,5 @@
 		this->inputs2->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
 	}
-	this->inputs2->SetDoubleInput(AreaEnum,this->lid,area);
-	this->AddInput2(SealevelAreaEnum,&area,P0Enum);
-
+	
 	/*Free allocations:*/
 	#ifdef _HAVE_RESTRICT_
@@ -5663,4 +5664,5 @@
 	bool notfullygrounded=false;
 	bool scaleoceanarea= false;
+	bool computeelastic= true;
 
 	/*output: */
@@ -5704,13 +5706,12 @@
 	#endif
 
-	/*recover material parameters: */
+	/*recover some parameters:*/
 	rho_ice=FindParam(MaterialsRhoIceEnum);
 	rho_water=FindParam(MaterialsRhoSeawaterEnum);
-
-	/*recover ocean area scaling: */
+	this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
 	this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
 
 	/*retrieve precomputed G:*/
-	this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
+	if(computeelastic)this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
 
 	/*Get area of element: precomputed in the sealevelrise_core_geometry:*/
@@ -5763,9 +5764,11 @@
 	_assert_(!xIsNan<IssmDouble>(bslrice));
 
-	/*convert from m to kg/m^2:*/
-	I=I*rho_ice*phi;
-
-	/*convolve:*/
-	for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*I;
+	if(computeelastic){
+		/*convert from m to kg/m^2:*/
+		I=I*rho_ice*phi;
+
+		/*convolve:*/
+		for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*I;
+	}
 
 	/*return :*/
@@ -5782,4 +5785,5 @@
 	bool notfullygrounded=false;
 	bool scaleoceanarea= false;
+	bool computeelastic= true;
 
 	/*elastic green function:*/
@@ -5807,13 +5811,12 @@
 	/*If we are here, we are on earth, not on ice: */
 
-	/*recover material parameters: */
+	/*recover parameters: */
 	rho_water=FindParam(MaterialsRhoSeawaterEnum);
 	rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum);
-
-	/*recover ocean area scaling: */
+	this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
 	this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
 
 	/*retrieve precomputed G:*/
-	this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
+	if(computeelastic)this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
 
 	/*Get area of element: precomputed in the sealevelrise_core_geometry:*/
@@ -5831,9 +5834,11 @@
 	_assert_(!xIsNan<IssmDouble>(bslrhydro));
 
-	/*convert from m to kg/m^2:*/
-	W=W*rho_freshwater*phi;
-
-	/*convolve:*/
-	for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*W;
+	if(computeelastic){
+		/*convert from m to kg/m^2:*/
+		W=W*rho_freshwater*phi;
+
+		/*convolve:*/
+		for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*W;
+	}
 
 	/*output:*/
@@ -5848,4 +5853,5 @@
 	IssmDouble BP;  //change in bottom pressure (Farrel and Clarke, Equ. 4)
 	IssmDouble constant;
+	bool computeelastic= true;
 
 	/*elastic green function:*/
@@ -5854,4 +5860,8 @@
 	/*water properties: */
 	IssmDouble rho_water;
+
+	/*exit now?:*/
+	this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
+	if(!computeelastic)return;
 
 	/*we are here to compute fingerprints originating fromn bottom pressure changes:*/
Index: /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25565)
+++ /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25566)
@@ -2,5 +2,4 @@
  * \brief: ISSM POST processing main program
  */ 
-
 /*includes and prototypes:*/
 #include "./issm.h"
@@ -693,5 +692,5 @@
 
 		/*open file: */
-		_printf0_("    opening file: " << file < "\n");
+		_printf0_("    opening file: " << file << "\n");
 		FILE* fid=fopen(file,"rb");
 		if(fid==NULL)_error_("cound not open file: " << file << "\n");