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

Last change on this file was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 3.6 KB
RevLine 
[25834]1Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24971)
4+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24972)
5@@ -5930,14 +5930,14 @@
6
7 /*diverse:*/
8 int gsize,dummy;
9- IssmDouble S; //change in water water level(Farrel and Clarke, Equ. 4)
10+ IssmDouble S,BP, Stotal; //change in water water level(Farrel and Clarke, Equ. 4)
11 IssmDouble constant=0;
12 IssmDouble rho_water;
13 IssmDouble* G=NULL;
14+ int bp_compute_fingerprints= 0;
15
16- /*optimization:*/
17- bool store_green_functions=false;
18-
19+ /*retrieve parameters:*/
20+ this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
21 rho_water=FindParam(MaterialsRhoSeawaterEnum);
22
23 /*early return if we are not on the ocean:*/
24@@ -5962,12 +5962,22 @@
25 /*From Sg_old, recover water sea level rise:*/
26 S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
27
28+ /*If we have bottom pressure fingerprinting requested, retrieve BP: */
29+ if(bp_compute_fingerprints){
30+ Input2* bottompressure_change_input=this->GetInput2(DslSeaWaterPressureChangeAtSeaFloor);
31+ if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level rise fingerprint!");
32+ bottompressure_change_input->GetInputAverage(&BP);
33+
34+ } else BP=0;
35+
36 /*convert to kg/m^2: */
37 S=S*rho_water;
38+ BP=BP*rho_water;
39
40- for(int i=0;i<gsize;i++) Sgo[i]+=G[i]*S;
41- //cblas_daxpy(gsize,S,&G[0],1,&Sgo[0],1); //for later.
42+ Stotal=S+BP;
43
44+ for(int i=0;i<gsize;i++) Sgo[i]+=G[i]*Stotal;
45+
46 return;
47 }
48 /*}}}*/
49@@ -5975,9 +5985,10 @@
50
51 /*diverse:*/
52 int gsize;
53- IssmDouble I, S; //change in relative ice thickness and sea level
54+ IssmDouble I, S, BP, Stotal; //change in relative ice thickness and sea level
55 IssmDouble rho_ice,rho_water;
56 int horiz;
57+ int bp_compute_fingerprints= 0;
58
59 /*precomputed elastic green functions:*/
60 IssmDouble* GU=NULL;
61@@ -5996,6 +6007,7 @@
62 /*recover parameters:*/
63 this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
64 this->parameters->FindParam(&horiz,SealevelriseHorizEnum);
65+ this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
66
67 /*early return if elastic not requested:*/
68 if(!computeelastic) return;
69@@ -6018,11 +6030,14 @@
70 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
71 deltathickness_input->GetInputAverage(&I);
72
73+ /*convert to kg/m^2*/
74+ I=I*rho_ice;
75+
76 for(int i=0;i<gsize;i++){
77- Up[i]+=rho_ice*I*GU[i];
78+ Up[i]+=I*GU[i];
79 if(horiz){
80- North[i]+=rho_ice*I*GN[i];
81- East[i]+=rho_ice*I*GE[i];
82+ North[i]+=I*GN[i];
83+ East[i]+=I*GE[i];
84 }
85 }
86 }
87@@ -6030,11 +6045,23 @@
88 /*From Sg, recover water sea level rise:*/
89 S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;
90
91+ /*If we have bottom pressure fingerprinting requested, retrieve BP: */
92+ if(bp_compute_fingerprints){
93+ Input2* bottompressure_change_input=this->GetInput2(DslSeaWaterPressureChangeAtSeaFloor);
94+ if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level rise fingerprint!");
95+ bottompressure_change_input->GetInputAverage(&BP);
96+ } else BP=0;
97+
98+ /*convert to kg/m^2:*/
99+ S=rho_water*S;
100+ BP=rho_water*BP;
101+ Stotal=S+BP;
102+
103 for(int i=0;i<gsize;i++){
104- Up[i]+=rho_water*S*GU[i];
105+ Up[i]+=Stotal*GU[i];
106 if(horiz){
107- North[i]+=rho_water*S*GN[i];
108- East[i]+=rho_water*S*GE[i];
109+ North[i]+=Stotal*GN[i];
110+ East[i]+=Stotal*GE[i];
111 }
112 }
113 }
Note: See TracBrowser for help on using the repository browser.