source: issm/oecreview/Archive/24684-25833/ISSM-24971-24972.diff

../trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -5930 +5930 @@
 
         /*diverse:*/
         int gsize,dummy;
-        IssmDouble S;  //change in water water level(Farrel and Clarke, Equ. 4)
+        IssmDouble S,BP, Stotal;  //change in water water level(Farrel and Clarke, Equ. 4)
         IssmDouble constant=0;
         IssmDouble rho_water;
         IssmDouble* G=NULL;
+        int  bp_compute_fingerprints= 0;
 
-        /*optimization:*/
-        bool store_green_functions=false;
-        
+        /*retrieve parameters:*/
+        this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
         rho_water=FindParam(MaterialsRhoSeawaterEnum);
 
         /*early return if we are not on the ocean:*/
@@ -5962 +5962 @@
         /*From Sg_old, recover water sea level rise:*/
         S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
 
+        /*If we have bottom pressure fingerprinting requested, retrieve BP: */
+        if(bp_compute_fingerprints){
+                Input2* bottompressure_change_input=this->GetInput2(DslSeaWaterPressureChangeAtSeaFloor);
+                if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level rise fingerprint!");
+                bottompressure_change_input->GetInputAverage(&BP);
+
+        } else BP=0;
+
         /*convert to kg/m^2: */
         S=S*rho_water;
+        BP=BP*rho_water;
 
-        for(int i=0;i<gsize;i++) Sgo[i]+=G[i]*S;
-        //cblas_daxpy(gsize,S,&G[0],1,&Sgo[0],1); //for later.
+        Stotal=S+BP;
 
+        for(int i=0;i<gsize;i++) Sgo[i]+=G[i]*Stotal;
+
         return;
 }
 /*}}}*/
@@ -5975 +5985 @@
 
         /*diverse:*/
         int gsize;
-        IssmDouble I, S;                //change in relative ice thickness and sea level
+        IssmDouble I, S, BP, Stotal;            //change in relative ice thickness and sea level
         IssmDouble rho_ice,rho_water;
         int horiz;
+        int  bp_compute_fingerprints= 0;
 
         /*precomputed elastic green functions:*/
         IssmDouble* GU=NULL;
@@ -5996 +6007 @@
         /*recover parameters:*/
         this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
         this->parameters->FindParam(&horiz,SealevelriseHorizEnum);
+        this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
 
         /*early return if elastic not requested:*/
         if(!computeelastic) return;
@@ -6018 +6030 @@
                 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
                 deltathickness_input->GetInputAverage(&I);
 
+                /*convert to kg/m^2*/
+                I=I*rho_ice;
+
                 for(int i=0;i<gsize;i++){
-                        Up[i]+=rho_ice*I*GU[i];
+                        Up[i]+=I*GU[i];
                         if(horiz){
-                                North[i]+=rho_ice*I*GN[i];
-                                East[i]+=rho_ice*I*GE[i];
+                                North[i]+=I*GN[i];
+                                East[i]+=I*GE[i];
                         }
                 }
         }
@@ -6030 +6045 @@
                 /*From Sg, recover water sea level rise:*/
                 S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;
 
+                /*If we have bottom pressure fingerprinting requested, retrieve BP: */
+                if(bp_compute_fingerprints){
+                        Input2* bottompressure_change_input=this->GetInput2(DslSeaWaterPressureChangeAtSeaFloor);
+                        if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level rise fingerprint!");
+                        bottompressure_change_input->GetInputAverage(&BP);
+                } else BP=0;
+
+                /*convert to kg/m^2:*/
+                S=rho_water*S;
+                BP=rho_water*BP;
+                Stotal=S+BP;
+
                 for(int i=0;i<gsize;i++){
-                        Up[i]+=rho_water*S*GU[i];
+                        Up[i]+=Stotal*GU[i];
                         if(horiz){
-                                North[i]+=rho_water*S*GN[i];
-                                East[i]+=rho_water*S*GE[i];
+                                North[i]+=Stotal*GN[i];
+                                East[i]+=Stotal*GE[i];
                         }
                 }
         }