Changeset 24972

Timestamp:

06/05/20 15:34:39 (5 years ago)

Author:

Eric.Larour

Message:

CHG: corrected wrong treatment of bottom pressure fingerprints.

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (6 diffs)

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

@@ -5931 +5931 @@
         /*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;
-
-        /*optimization:*/
-        bool store_green_functions=false;
-        
+        int  bp_compute_fingerprints= 0;
+
+        /*retrieve parameters:*/
+        this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
         rho_water=FindParam(MaterialsRhoSeawaterEnum);
 
@@ -5963 +5963 @@
         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;
-
-        for(int i=0;i<gsize;i++) Sgo[i]+=G[i]*S;
-        //cblas_daxpy(gsize,S,&G[0],1,&Sgo[0],1); //for later.
+        BP=BP*rho_water;
+
+        Stotal=S+BP;
+
+        for(int i=0;i<gsize;i++) Sgo[i]+=G[i]*Stotal;
 
         return;
@@ -5976 +5986 @@
         /*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:*/
@@ -5997 +6008 @@
         this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
         this->parameters->FindParam(&horiz,SealevelriseHorizEnum);
+        this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
 
         /*early return if elastic not requested:*/
@@ -6019 +6031 @@
                 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];
                         }
                 }
@@ -6031 +6046 @@
                 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];
                         }
                 }