source: issm/oecreview/Archive/19101-20495/ISSM-20043-20044.diff@ 20498

Index: ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp        (revision 20043)
+++ ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp        (revision 20044)
@@ -70,8 +70,29 @@
                iomodel->FetchData(&degacc,SealevelriseDegaccEnum);
                M=reCast<int,IssmDouble>(180./degacc+1.);
                G_elastic=xNew<IssmDouble>(M);
+               
+               /*compute combined legendre + love number (elastic green function:*/
+               #ifdef _HAVE_ADOLC_
+               for(int i=0;i<M;i++){ 
+                       IssmDouble alpha,x;
+                       alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
 
-               /*compute combined legendre + love number (elastic green function:*/
+                       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;
+                       }
+               }
+               #else
                m=DetermineLocalSize(M,IssmComm::GetComm());
                GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
                G_elastic_local=xNew<IssmDouble>(m);
@@ -113,29 +134,8 @@
                xDelete<int>(displs);
 
                /*}}}*/
+               #endif
 
-               /*Old code: {{{*/
-               /*for(int i=0;i<M;i++){ 
-                       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));