source: issm/oecreview/Archive/22819-23185/ISSM-23155-23156.diff

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

CHG: added Archive/22819-23185
File size: 12.5 KB
RevLine 
[23186]1Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
2===================================================================
3--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23155)
4+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23156)
5@@ -450,7 +450,7 @@
6 FrictionPressureAdjustedTemperatureEnum,
7 FrictionQEnum,
8 FrictionWaterLayerEnum,
9- FrictionWatercolumnMaxEnum,
10+ HydrologyWatercolumnMaxEnum,
11 FrictionTillFrictionAngleEnum,
12 FrictionSedimentCompressibilityCoefficientEnum,
13 GeometryHydrostaticRatioEnum,
14Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
15===================================================================
16--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23155)
17+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23156)
18@@ -456,7 +456,7 @@
19 case FrictionPressureAdjustedTemperatureEnum : return "FrictionPressureAdjustedTemperature";
20 case FrictionQEnum : return "FrictionQ";
21 case FrictionWaterLayerEnum : return "FrictionWaterLayer";
22- case FrictionWatercolumnMaxEnum : return "FrictionWatercolumnMax";
23+ case HydrologyWatercolumnMaxEnum : return "HydrologyWatercolumnMax";
24 case FrictionTillFrictionAngleEnum : return "FrictionTillFrictionAngle";
25 case FrictionSedimentCompressibilityCoefficientEnum : return "FrictionSedimentCompressibilityCoefficient";
26 case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio";
27Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
28===================================================================
29--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23155)
30+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23156)
31@@ -465,7 +465,7 @@
32 else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;
33 else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
34 else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
35- else if (strcmp(name,"FrictionWatercolumnMax")==0) return FrictionWatercolumnMaxEnum;
36+ else if (strcmp(name,"HydrologyWatercolumnMax")==0) return HydrologyWatercolumnMaxEnum;
37 else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
38 else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
39 else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
40Index: ../trunk-jpl/src/c/classes/Loads/Friction.cpp
41===================================================================
42--- ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23155)
43+++ ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23156)
44@@ -709,21 +709,11 @@
45 element->parameters->FindParam(&e0,FrictionVoidRatioEnum);
46 element->GetInputValue(&Cc,gauss,FrictionSedimentCompressibilityCoefficientEnum);
47 element->GetInputValue(&W,gauss,WatercolumnEnum);
48- element->GetInputValue(&Wmax,gauss,FrictionWatercolumnMaxEnum);
49+ element->GetInputValue(&Wmax,gauss,HydrologyWatercolumnMaxEnum);
50
51-// /*Check that water column height is within 0 and upper bound, correct if needed*/
52-// // if watercolumn height is higher than the maximum allowed height, set height to upper bound
53-// if(W>Wmax){
54-// W=Wmax;
55-// }
56-// // if watercolumn height is negative (shouldn't happen), set it to 0
57-// if(W<0){
58-// W=0;
59-// }
60-// // if watercolumn height is within 0 and upper bound, nothing to be done
61-// else{
62-// //do nothing
63-// }
64+ /*Check that water column height is within 0 and upper bound, correct if needed*/
65+ if(W>Wmax) W=Wmax;
66+ if(W<0) W=0.;
67
68 N = delta*P0*pow(10.,(e0/Cc)*(1.-W/Wmax));
69
70@@ -745,7 +735,7 @@
71 IssmDouble u0,q;
72 element->parameters->FindParam(&u0,FrictionThresholdSpeedEnum);
73 element->parameters->FindParam(&q,FrictionPseudoplasticityExponentEnum);
74- IssmDouble alpha2 = tau_c/(pow(ub,1.-q)*pow(u0,q));
75+ IssmDouble alpha2 = tau_c/(pow(ub+1.e-10,1.-q)*pow(u0,q));
76
77 /*Final checks in debuging mode*/
78 _assert_(!xIsNan<IssmDouble>(alpha2));
79Index: ../trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp
80===================================================================
81--- ../trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp (revision 23155)
82+++ ../trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp (revision 23156)
83@@ -35,6 +35,7 @@
84 iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
85 iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
86 iomodel->FetchDataToInput(elements,"md.hydrology.drainage_rate",HydrologyDrainageRateEnum);
87+ iomodel->FetchDataToInput(elements,"md.hydrology.watercolumn_max",HydrologyWatercolumnMaxEnum);
88 iomodel->FetchDataToInput(elements,"md.initialization.watercolumn",WatercolumnEnum,0.);
89 }/*}}}*/
90 void HydrologyPismAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
91@@ -114,28 +115,18 @@
92 IssmDouble* watercolumn = xNew<IssmDouble>(numvertices);
93 IssmDouble* drainagerate = xNew<IssmDouble>(numvertices);
94 IssmDouble* meltingrate = xNew<IssmDouble>(numvertices);
95-// IssmDouble* watercolumn_max = xNew<IssmDouble>(numvertices);
96+ IssmDouble* watercolumn_max = xNew<IssmDouble>(numvertices);
97 element->GetInputListOnVertices(&watercolumn[0],WatercolumnEnum);
98 element->GetInputListOnVertices(&drainagerate[0],HydrologyDrainageRateEnum);
99 element->GetInputListOnVertices(&meltingrate[0],BasalforcingsGroundediceMeltingRateEnum);
100-// element->GetInputListOnVertices(&watercolumn_max[0],FrictionWatercolumnMaxEnum);
101+ element->GetInputListOnVertices(&watercolumn_max[0],HydrologyWatercolumnMaxEnum);
102
103 /*Add water*/
104-// /*Check that water column height is within 0 and upper bound, correct if needed*/
105 for(int i=0;i<numvertices;i++){
106 watercolumn[i] += (meltingrate[i]/rho_ice*rho_water-drainagerate[i])*dt;
107-// // if watercolumn height is higher than the maximum allowed height, set height to upper bound
108-// if(watercolumn[i]>watercolumn_max[i]){
109-// watercolumn[i]=watercolumn_max[i];
110-// }
111-// // if watercolumn height is negative (shouldn't happen), set it to 0
112-// if(watercolumn[i]<0){
113-// watercolumn[i]=0;
114-// }
115-// // if watercolumn height is within 0 and upper bound, nothing to be done
116-// else{
117-// //do nothing
118-// }
119+ /*Check that water column height is within 0 and upper bound, correct if needed*/
120+ if(watercolumn[i]>watercolumn_max[i]) watercolumn[i]=watercolumn_max[i];
121+ if(watercolumn[i]<0) watercolumn[i]=0.;
122 }
123
124 /* Divide by connectivity, add degree of channelization as an input */
125@@ -143,5 +134,7 @@
126
127 /*Clean up and return*/
128 xDelete<IssmDouble>(watercolumn);
129+ xDelete<IssmDouble>(meltingrate);
130+ xDelete<IssmDouble>(watercolumn_max);
131 xDelete<IssmDouble>(drainagerate);
132 }/*}}}*/
133Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
134===================================================================
135--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23155)
136+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23156)
137@@ -847,7 +847,6 @@
138 case 10:
139 iomodel->FetchDataToInput(elements,"md.friction.till_friction_angle",FrictionTillFrictionAngleEnum);
140 iomodel->FetchDataToInput(elements,"md.friction.sediment_compressibility_coefficient",FrictionSedimentCompressibilityCoefficientEnum);
141- iomodel->FetchDataToInput(elements,"md.friction.watercolumn_max",FrictionWatercolumnMaxEnum);
142 break;
143 default:
144 _error_("friction law "<< frictionlaw <<" not supported");
145Index: ../trunk-jpl/src/m/classes/hydrologypism.m
146===================================================================
147--- ../trunk-jpl/src/m/classes/hydrologypism.m (revision 23155)
148+++ ../trunk-jpl/src/m/classes/hydrologypism.m (revision 23156)
149@@ -5,7 +5,8 @@
150
151 classdef hydrologypism
152 properties (SetAccess=public)
153- drainage_rate = NaN;
154+ drainage_rate = NaN;
155+ watercolumn_max = NaN;
156 end
157 methods
158 function self = extrude(self,md) % {{{
159@@ -35,10 +36,12 @@
160 end
161
162 md = checkfield(md,'fieldname','hydrology.drainage_rate','Inf',1,'NaN',1,'>=',0,'size',[md.mesh.numberofvertices 1]);
163+ md = checkfield(md,'fieldname','friction.watercolumn_max','NaN',1,'Inf',1,'>',0.,'size',[md.mesh.numberofvertices 1]);
164 end % }}}
165 function disp(self) % {{{
166 disp(sprintf(' hydrologypism solution parameters:'));
167 fielddisplay(self,'drainage_rate','fixed drainage rate [mm/yr]');
168+ fielddisplay(self,'watercolumn_max','maximum water column height [m], recommended default: 2 m');
169 end % }}}
170 function marshall(self,prefix,md,fid) % {{{
171
172@@ -46,6 +49,7 @@
173
174 WriteData(fid,prefix,'name','md.hydrology.model','data',4,'format','Integer');
175 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','drainage_rate','format','DoubleMat','mattype',1,'scale',1./(1000.*yts)); %from mm/yr to m/s
176+ WriteData(fid,prefix,'class','hydrology','object',self,'fieldname','watercolumn_max','format','DoubleMat','mattype',1);
177 WriteData(fid,prefix,'data',{'Watercolumn'},'name','md.hydrology.requested_outputs','format','StringArray');
178 end % }}}
179 end
180Index: ../trunk-jpl/src/m/classes/frictionpism.m
181===================================================================
182--- ../trunk-jpl/src/m/classes/frictionpism.m (revision 23155)
183+++ ../trunk-jpl/src/m/classes/frictionpism.m (revision 23156)
184@@ -11,7 +11,6 @@
185 void_ratio = 0.;
186 till_friction_angle = NaN;
187 sediment_compressibility_coefficient = NaN;
188- watercolumn_max = NaN;
189 end
190 methods
191 function self = extrude(self,md) % {{{
192@@ -48,7 +47,6 @@
193 md = checkfield(md,'fieldname','friction.void_ratio','numel',[1],'>',0,'<',1,'NaN',1,'Inf',1);
194 md = checkfield(md,'fieldname','friction.till_friction_angle','NaN',1,'Inf',1,'<',360.,'>',0.,'size',[md.mesh.numberofvertices 1]); %User should give angle in degrees, Matlab calculates in rad
195 md = checkfield(md,'fieldname','friction.sediment_compressibility_coefficient','NaN',1,'Inf',1,'<',1.,'>',0.,'size',[md.mesh.numberofvertices 1]);
196- md = checkfield(md,'fieldname','friction.watercolumn_max','NaN',1,'Inf',1,'>',0.,'size',[md.mesh.numberofvertices 1]);
197 end % }}}
198 function disp(self) % {{{
199 disp(sprintf('Basal shear stress parameters for the PISM friction law (See Aschwanden et al. 2016 for more details)'));
200@@ -58,7 +56,6 @@
201 fielddisplay(self,'void_ratio','void ratio at a reference effective pressure [dimensionless]');
202 fielddisplay(self,'till_friction_angle','till friction angle [deg], recommended default: 30 deg');
203 fielddisplay(self,'sediment_compressibility_coefficient','coefficient of compressibility of the sediment [dimensionless], recommended default: 0.12');
204- fielddisplay(self,'watercolumn_max','maximum water column height [m], recommended default: 2 m');
205 end % }}}
206 function marshall(self,prefix,md,fid) % {{{
207 yts=md.constants.yts;
208@@ -70,7 +67,6 @@
209 WriteData(fid,prefix,'class','friction','object',self,'fieldname','void_ratio','format','Double');
210 WriteData(fid,prefix,'class','friction','object',self,'fieldname','till_friction_angle','format','DoubleMat','mattype',1);
211 WriteData(fid,prefix,'class','friction','object',self,'fieldname','sediment_compressibility_coefficient','format','DoubleMat','mattype',1);
212- WriteData(fid,prefix,'class','friction','object',self,'fieldname','watercolumn_max','format','DoubleMat','mattype',1);
213 end % }}}
214 function savemodeljs(self,fid,modelname) % {{{
215 error('not implemented yet!');
216Index: ../trunk-jpl/test/NightlyRun/test612.m
217===================================================================
218--- ../trunk-jpl/test/NightlyRun/test612.m (revision 23155)
219+++ ../trunk-jpl/test/NightlyRun/test612.m (revision 23156)
220@@ -6,7 +6,8 @@
221
222 %Hydrology
223 md.hydrology = hydrologypism();
224-md.hydrology.drainage_rate = 0.001*ones(md.mesh.numberofvertices,1);
225+md.hydrology.drainage_rate = 0.001*ones(md.mesh.numberofvertices,1);
226+md.hydrology.watercolumn_max=2*ones(md.mesh.numberofvertices,1);
227 md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
228 md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1);
229 md.basalforcings.groundedice_melting_rate = [1:md.mesh.numberofvertices]';
230@@ -15,7 +16,6 @@
231 md.friction=frictionpism();
232 md.friction.till_friction_angle=30*ones(md.mesh.numberofvertices,1);
233 md.friction.sediment_compressibility_coefficient=0.12*ones(md.mesh.numberofvertices,1);
234-md.friction.watercolumn_max=2*ones(md.mesh.numberofvertices,1);
235
236 md.transient.ishydrology = 1;
237 md.transient.issmb = 0;
Note: See TracBrowser for help on using the repository browser.