Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp	(revision 20028)
+++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp	(revision 20029)
@@ -43,4 +43,5 @@
 	IssmDouble* love_k=NULL;
 	
+	bool elastic=false;
 	IssmDouble* G_elastic = NULL;
 	int         M;
@@ -54,41 +55,48 @@
 	parameters->AddObject(iomodel->CopyConstantObject(SealevelriseEustaticEnum));
 
-	/*love numbers: */
-	iomodel->FetchData(&love_h,&nl,NULL,SealevelriseLoveHEnum);
-	iomodel->FetchData(&love_k,&nl,NULL,SealevelriseLoveKEnum);
+	iomodel->FetchData(&elastic,SealevelriseElasticEnum);
+	if(elastic){
 
-	/*compute elastic green function for a range of angles*/
-	const IssmDouble degacc=.01; M=reCast<int>(180/degacc+1);
-	G_elastic=xNew<IssmDouble>(M);
-	
-	/*compute combined legendre + love number (elastic green function:*/
-	for(int i=0;i<M;i++){ //watch out the <=
-		IssmDouble alpha,x;
-		alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
-		
-		G_elastic[i]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
-		for (int n=0;n<nl;n++) {
-			IssmDouble Pn,Pn1,Pn2;
-			IssmDouble deltalove;
-			
-			deltalove = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
+		/*love numbers: */
+		iomodel->FetchData(&love_h,&nl,NULL,SealevelriseLoveHEnum);
+		iomodel->FetchData(&love_k,&nl,NULL,SealevelriseLoveKEnum);
 
-			if(n==0)Pn=1;
-			else if(n==1)Pn=cos(alpha);
-			else Pn= ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
-			Pn2=Pn1; Pn1=Pn;
+		/*compute elastic green function for a range of angles*/
+		const IssmDouble degacc=.01; M=reCast<int>(180/degacc+1);
+		G_elastic=xNew<IssmDouble>(M);
 
-			G_elastic[i] += deltalove*Pn;
+		/*compute combined legendre + love number (elastic green function:*/
+		for(int i=0;i<M;i++){ //watch out the <=
+			IssmDouble alpha,x;
+			alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
+
+			G_elastic[i]= (love_k[nl-1]-love_h[nl-1])/2.0/sin(alpha/2.0);
+			for (int n=0;n<nl;n++) {
+				IssmDouble Pn,Pn1,Pn2;
+				IssmDouble deltalove;
+
+				deltalove = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
+
+				if(n==0)Pn=1;
+				else if(n==1)Pn=cos(alpha);
+				else Pn= ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
+				Pn2=Pn1; Pn1=Pn;
+
+				G_elastic[i] += deltalove*Pn;
+			}
 		}
+
+		/*Avoid singularity at 0: */
+		G_elastic[0]=G_elastic[1];
+		parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
+
+		/*free ressources:*/
+		xDelete<IssmDouble>(G_elastic);
+
+		/*free ressources: */
+		xDelete<IssmDouble>(love_h);
+		xDelete<IssmDouble>(love_k);
+		xDelete<IssmDouble>(G_elastic);
 	}
-	parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
-
-	/*free ressources:*/
-	xDelete<IssmDouble>(G_elastic);
-
-	/*free ressources: */
-	xDelete<IssmDouble>(love_h);
-	xDelete<IssmDouble>(love_k);
-	xDelete<IssmDouble>(G_elastic);
 
 }/*}}}*/