source: issm/oecreview/Archive/21724-22754/ISSM-21793-21794.diff@ 22755

Last change on this file since 22755 was 22755, checked in by Mathieu Morlighem, 7 years ago

CHG: added 21724-22754
File size: 4.3 KB
RevLine 
[22755]1Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 21793)
4+++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 21794)
5@@ -469,7 +469,8 @@
6
7 /*Intermediaries */
8 IssmDouble Jdet,dt;
9- IssmDouble ms,mb,gmb,fmb,thickness,phi;
10+ IssmDouble ms,mb,gmb,fmb,thickness;
11+ IssmDouble gllevelset,phi=1.;
12 IssmDouble* xyz_list = NULL;
13
14 /*Fetch number of nodes and dof for this finite element*/
15@@ -484,7 +485,7 @@
16 element->FindParam(&dt,TimesteppingTimeStepEnum);
17 Input* gmb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
18 Input* fmb_input = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
19- Input* groundedice_input = element->GetInput(MaskGroundediceLevelsetEnum); _assert_(groundedice_input);
20+ Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
21 Input* ms_input = element->GetInput(SmbMassBalanceEnum); _assert_(ms_input);
22 Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
23
24@@ -499,7 +500,7 @@
25 ms_input->GetInputValue(&ms,gauss);
26 gmb_input->GetInputValue(&gmb,gauss);
27 fmb_input->GetInputValue(&fmb,gauss);
28- groundedice_input->GetInputValue(&phi,gauss);
29+ gllevelset_input->GetInputValue(&phi,gauss);
30 thickness_input->GetInputValue(&thickness,gauss);
31 if(phi>0.) mb=gmb;
32 else mb=fmb;
33@@ -519,9 +520,13 @@
34 if(!element->IsIceInElement()) return NULL;
35
36 /*Intermediaries */
37- IssmDouble Jdet,dt;
38- IssmDouble ms,mb,gmb,fmb,thickness,phi;
39+ int migration_style,point1;
40+ bool mainlyfloating;
41+ IssmDouble Jdet,fraction1,fraction2;
42+ IssmDouble dt,gllevelset;
43+ IssmDouble ms,mb,gmb,fmb,thickness,phi=1.;
44 IssmDouble* xyz_list = NULL;
45+ Gauss* gauss = NULL;
46
47 /*Fetch number of nodes and dof for this finite element*/
48 int numnodes = element->GetNumberOfNodes();
49@@ -536,11 +541,21 @@
50 Input* gmb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
51 Input* fmb_input = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
52 Input* ms_input = element->GetInput(SmbMassBalanceEnum); _assert_(ms_input);
53- Input* groundedice_input = element->GetInput(MaskGroundediceLevelsetEnum); _assert_(groundedice_input);
54+ Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
55 Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
56
57+ /*Recover portion of element that is grounded*/
58+ phi=element->GetGroundedPortion(xyz_list);
59+ if(migration_style==SubelementMigration2Enum){
60+ gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
61+ element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
62+ gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,3);
63+ }
64+ else{
65+ gauss = element->NewGauss(3);
66+ }
67+
68 /* Start looping on the number of gaussian points: */
69- Gauss* gauss=element->NewGauss(2);
70 for(int ig=gauss->begin();ig<gauss->end();ig++){
71 gauss->GaussPoint(ig);
72
73@@ -550,10 +565,24 @@
74 ms_input->GetInputValue(&ms,gauss);
75 gmb_input->GetInputValue(&gmb,gauss);
76 fmb_input->GetInputValue(&fmb,gauss);
77- groundedice_input->GetInputValue(&phi,gauss);
78+ gllevelset_input->GetInputValue(&gllevelset,gauss);
79 thickness_input->GetInputValue(&thickness,gauss);
80- if(phi>0) mb=gmb;
81- else mb=fmb;
82+ if(migration_style==SubelementMigrationEnum){
83+ if (phi<0.00000001) mb=gmb;
84+ else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
85+ }
86+ else if(migration_style==SubelementMigration2Enum){
87+ if(gllevelset>0.) mb=gmb;
88+ else mb=fmb;
89+ }
90+ else if(migration_style==SubelementMigration3Enum){
91+ if (phi<0.00000001) mb=gmb;
92+ else mb=gmb;
93+ }
94+ else{
95+ if(gllevelset>0) mb=gmb;
96+ else mb=fmb;
97+ }
98
99 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
100 }
Note: See TracBrowser for help on using the repository browser.