Changeset 20029

Timestamp:

01/30/16 14:51:22 (9 years ago)

Author:

Eric.Larour

Message:

CHG: do not compute G_elastic if we are not in elastic mode.
Also, avoid singularity due to discretization of low incidence angles.

File:

• issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp

@@ -43 +43 @@
         IssmDouble* love_k=NULL;
         
+        bool elastic=false;
         IssmDouble* G_elastic = NULL;
         int         M;
@@ -54 +55 @@
         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);
 
 }/*}}}*/