source: issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp@ 27510

Last change on this file since 27510 was 27510, checked in by inwoo, 2 years ago

CHG: fix SMB with SEMIC.

1) Update display of SMBsemic.m
2) Enable to use different albedo scheme in SEMIC.
3) Surface mass balance is from SMB_ice in SEMIC.
4) Enable to return surface mass balance of snow in SEMIC.
File size: 29.8 KB
Line 
1#include "./SmbAnalysis.h"
2#include "../toolkits/toolkits.h"
3#include "../classes/classes.h"
4#include "../shared/shared.h"
5#include "../modules/modules.h"
6
7// FIX
8#include "./shared/io/Print/Print.h"
9
10/*Model processing*/
11void SmbAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
12 /*No constraints*/
13}/*}}}*/
14void SmbAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
15 /*No loads*/
16}/*}}}*/
17void SmbAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
18 ::CreateNodes(nodes,iomodel,SmbAnalysisEnum,P1Enum);
19}/*}}}*/
20int SmbAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
21 return 1;
22}/*}}}*/
23void SmbAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
24
25 int smb_model;
26 bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled,isfirnwarming,isstochastic;
27
28 /*Update elements: */
29 int counter=0;
30 for(int i=0;i<iomodel->numberofelements;i++){
31 if(iomodel->my_elements[i]){
32 Element* element=(Element*)elements->GetObjectByOffset(counter);
33 element->Update(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum);
34 counter++;
35 }
36 }
37
38 /*Figure out smb model: */
39 iomodel->FindConstant(&smb_model,"md.smb.model");
40 iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing");
41 switch(smb_model){
42 case SMBforcingEnum:
43 iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
44 if(isstochastic){
45 iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
46 iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",BaselineSmbMassBalanceEnum,0.);
47 }
48 break;
49 case SMBgembEnum:
50 iomodel->FetchDataToInput(inputs,elements,"md.smb.Ta",SmbTaEnum);
51 iomodel->FetchDataToInput(inputs,elements,"md.smb.V",SmbVEnum);
52 iomodel->FetchDataToInput(inputs,elements,"md.smb.dswrf",SmbDswrfEnum);
53 iomodel->FetchDataToInput(inputs,elements,"md.smb.dswdiffrf",SmbDswdiffrfEnum);
54 iomodel->FetchDataToInput(inputs,elements,"md.smb.dlwrf",SmbDlwrfEnum);
55 iomodel->FetchDataToInput(inputs,elements,"md.smb.P",SmbPEnum);
56 iomodel->FetchDataToInput(inputs,elements,"md.smb.eAir",SmbEAirEnum);
57 iomodel->FetchDataToInput(inputs,elements,"md.smb.pAir",SmbPAirEnum);
58 iomodel->FetchDataToInput(inputs,elements,"md.smb.zTop",SmbZTopEnum);
59 iomodel->FetchDataToInput(inputs,elements,"md.smb.dzTop",SmbDzTopEnum);
60 iomodel->FetchDataToInput(inputs,elements,"md.smb.dzMin",SmbDzMinEnum);
61 iomodel->FetchDataToInput(inputs,elements,"md.smb.zY",SmbZYEnum);
62 iomodel->FetchDataToInput(inputs,elements,"md.smb.zMax",SmbZMaxEnum);
63 iomodel->FetchDataToInput(inputs,elements,"md.smb.zMin",SmbZMinEnum);
64 iomodel->FetchDataToInput(inputs,elements,"md.smb.Tmean",SmbTmeanEnum);
65 iomodel->FetchDataToInput(inputs,elements,"md.smb.Vmean",SmbVmeanEnum);
66 iomodel->FetchDataToInput(inputs,elements,"md.smb.C",SmbCEnum);
67 iomodel->FetchDataToInput(inputs,elements,"md.smb.Tz",SmbTzEnum);
68 iomodel->FetchDataToInput(inputs,elements,"md.smb.Vz",SmbVzEnum);
69 InputUpdateFromConstantx(inputs,elements,false,SmbIsInitializedEnum);
70 iomodel->FetchDataToInput(inputs,elements,"md.smb.Dzini",SmbDziniEnum);
71 iomodel->FetchDataToInput(inputs,elements,"md.smb.Dini",SmbDiniEnum);
72 iomodel->FetchDataToInput(inputs,elements,"md.smb.Reini",SmbReiniEnum);
73 iomodel->FetchDataToInput(inputs,elements,"md.smb.Gdnini",SmbGdniniEnum);
74 iomodel->FetchDataToInput(inputs,elements,"md.smb.Gspini",SmbGspiniEnum);
75 iomodel->FetchDataToInput(inputs,elements,"md.smb.ECini",SmbECiniEnum);
76 iomodel->FetchDataToInput(inputs,elements,"md.smb.Wini",SmbWiniEnum);
77 iomodel->FetchDataToInput(inputs,elements,"md.smb.Aini",SmbAiniEnum);
78 iomodel->FetchDataToInput(inputs,elements,"md.smb.Adiffini",SmbAdiffiniEnum);
79 iomodel->FetchDataToInput(inputs,elements,"md.smb.Tini",SmbTiniEnum);
80 iomodel->FetchDataToInput(inputs,elements,"md.smb.Sizeini",SmbSizeiniEnum);
81 iomodel->FetchDataToInput(inputs,elements,"md.smb.aValue",SmbAValueEnum);
82 iomodel->FetchDataToInput(inputs,elements,"md.smb.dulwrfValue",SmbDulwrfValueEnum);
83 iomodel->FetchDataToInput(inputs,elements,"md.smb.teValue",SmbTeValueEnum);
84 iomodel->FetchDataToInput(inputs,elements,"md.smb.szaValue",SmbSzaValueEnum);
85 iomodel->FetchDataToInput(inputs,elements,"md.smb.cotValue",SmbCotValueEnum);
86 iomodel->FetchDataToInput(inputs,elements,"md.smb.ccsnowValue",SmbCcsnowValueEnum);
87 iomodel->FetchDataToInput(inputs,elements,"md.smb.cciceValue",SmbCciceValueEnum);
88 break;
89 case SMBpddEnum:
90 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
91 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
92 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
93 iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
94 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
95 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
96 if(isdelta18o || ismungsm){
97 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_lgm",SmbTemperaturesLgmEnum);
98 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
99 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
100 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum);
101 }else{
102 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitation",SmbPrecipitationEnum);
103 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
104 }
105 if(issetpddfac){
106 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
107 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
108 }
109 break;
110 case SMBpddSicopolisEnum:
111 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
112 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
113 iomodel->FindConstant(&isfirnwarming,"md.smb.isfirnwarming");
114 iomodel->FetchDataToInput(inputs,elements,"md.smb.smb_corr",SmbSmbCorrEnum);
115 iomodel->FetchDataToInput(inputs,elements,"md.smb.precipitation_anomaly",SmbPrecipitationsAnomalyEnum);
116 iomodel->FetchDataToInput(inputs,elements,"md.smb.temperature_anomaly",SmbTemperaturesAnomalyEnum);
117 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
118 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitation",SmbPrecipitationEnum);
119 break;
120 case SMBd18opddEnum:
121 iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
122 iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
123 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
124 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
125 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
126 iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
127 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
128 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
129 if(isd18opd){
130 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
131 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
132 if(!istemperaturescaled){
133 /*Fetch array*/
134 IssmDouble* doublearray = NULL;
135 int M,N;
136 iomodel->FetchData(&doublearray,&M,&N,"md.smb.temperatures_reconstructed");
137 if(M!=iomodel->numberofvertices+1) _error_("md.smb.temperatures_reconstructed should have nbv+1 rows");
138 if(N%12!=0) _error_("md.smb.temperatures_reconstructed should have a multiple of 12 columns (since it is monthly)");
139
140 /*Build indices*/
141 int* ids = xNew<int>(N); for(int i=0;i<N;i++) ids[i] = i;
142
143 for(Object* & object : elements->objects){
144 Element* element=xDynamicCast<Element*>(object);
145 element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs,iomodel,SmbTemperaturesReconstructedEnum);
146 }
147 xDelete<int>(ids);
148 iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
149 }
150 if(!isprecipscaled){
151 /*Fetch array*/
152 IssmDouble* doublearray = NULL;
153 int M,N;
154 iomodel->FetchData(&doublearray,&M,&N,"md.smb.precipitations_reconstructed");
155 if(M!=iomodel->numberofvertices+1) _error_("md.smb.precipitations_reconstructed should have nbv+1 rows");
156 if(N%12!=0) _error_("md.smb.precipitations_reconstructed should have a multiple of 12 columns (since it is monthly)");
157
158 /*Build indices*/
159 int* ids = xNew<int>(N); for(int i=0;i<N;i++) ids[i] = i;
160
161 for(Object* & object : elements->objects){
162 Element* element=xDynamicCast<Element*>(object);
163 element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs,iomodel,SmbPrecipitationsReconstructedEnum);
164 }
165 xDelete<int>(ids);
166 iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
167 }
168 }
169 if(issetpddfac){
170 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
171 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
172 }
173 break;
174 case SMBgradientsEnum:
175 iomodel->FetchDataToInput(inputs,elements,"md.smb.href",SmbHrefEnum);
176 iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum);
177 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_pos",SmbBPosEnum);
178 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_neg",SmbBNegEnum);
179 break;
180 case SMBgradientselaEnum:
181 iomodel->FetchDataToInput(inputs,elements,"md.smb.ela",SmbElaEnum);
182 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_pos",SmbBPosEnum);
183 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_neg",SmbBNegEnum);
184 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_max",SmbBMaxEnum);
185 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_min",SmbBMinEnum);
186 break;
187 case SMBarmaEnum:
188 iomodel->FetchDataToInput(inputs,elements,"md.smb.basin_id",SmbBasinsIdEnum);
189 break;
190 case SMBhenningEnum:
191 iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum,0.);
192 break;
193 case SMBcomponentsEnum:
194 iomodel->FetchDataToInput(inputs,elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
195 iomodel->FetchDataToInput(inputs,elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
196 iomodel->FetchDataToInput(inputs,elements,"md.smb.runoff",SmbRunoffEnum,0.);
197 break;
198 case SMBmeltcomponentsEnum:
199 iomodel->FetchDataToInput(inputs,elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
200 iomodel->FetchDataToInput(inputs,elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
201 iomodel->FetchDataToInput(inputs,elements,"md.smb.melt",SmbMeltEnum,0.);
202 iomodel->FetchDataToInput(inputs,elements,"md.smb.refreeze",SmbRefreezeEnum,0.);
203 break;
204 case SMBgradientscomponentsEnum:
205 /* Nothing to add to input */
206 break;
207 case SMBsemicEnum:
208 int ismethod;
209 //if(VerboseSolution()) _printf0_(" smb semic: UpdateElements.\n");
210 //iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
211 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0gcm",SmbS0gcmEnum);
212 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailysnowfall",SmbDailysnowfallEnum);
213 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyrainfall",SmbDailyrainfallEnum);
214 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailydsradiation",SmbDailydsradiationEnum);
215 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailydlradiation",SmbDailydlradiationEnum);
216 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailywindspeed",SmbDailywindspeedEnum);
217 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailypressure",SmbDailypressureEnum);
218 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyairdensity",SmbDailyairdensityEnum);
219 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyairhumidity",SmbDailyairhumidityEnum);
220 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailytemperature",SmbDailytemperatureEnum);
221 // assign initial SEMIC temperature from initialization class.
222 if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - temperature.\n");
223 iomodel->FetchDataToInput(inputs,elements,"md.initialization.temperature",TemperatureEnum);
224
225 iomodel->FindConstant(&ismethod,"md.smb.ismethod");
226 if (ismethod == 1){
227 // initial albedo: TODO: this is not needed.
228 //if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - albedo.\n");
229 iomodel->FetchDataToInput(inputs,elements,"md.smb.albedo",SmbAlbedoInitEnum);
230 iomodel->FetchDataToInput(inputs,elements,"md.smb.albedo_snow",SmbAlbedoSnowInitEnum);
231 if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - Hice/Hsnow.\n");
232 iomodel->FetchDataToInput(inputs,elements,"md.smb.hice",SmbHIceInitEnum);
233 iomodel->FetchDataToInput(inputs,elements,"md.smb.hsnow",SmbHSnowInitEnum);
234
235 // initial Temperature amplitude.
236 if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - Tamp.\n");
237 iomodel->FetchDataToInput(inputs,elements,"md.smb.Tamp",SmbTampEnum);
238
239 // assign masking
240 iomodel->FetchDataToInput(inputs,elements,"md.smb.mask",SmbMaskEnum);
241 if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - done.\n");
242 }
243 break;
244 case SMBdebrisMLEnum:
245 iomodel->FetchDataToInput(inputs,elements,"md.initialization.debris",DebrisThicknessEnum);
246 break;
247 default:
248 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
249 }
250
251}/*}}}*/
252void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
253
254 int numoutputs;
255 char** requestedoutputs = NULL;
256 bool isdelta18o,ismungsm,isd18opd,issetpddfac,interp,cycle,isfirnwarming;
257 int smb_model, smbslices, averaging;
258 IssmDouble *temp = NULL;
259 int N,M;
260
261 parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
262
263 iomodel->FindConstant(&smb_model,"md.smb.model");
264 iomodel->FindConstant(&interp,"md.timestepping.interp_forcing");
265 iomodel->FindConstant(&cycle,"md.timestepping.cycle_forcing");
266
267 iomodel->FindConstant(&smbslices,"md.smb.steps_per_step");
268 parameters->AddObject(new IntParam(SmbStepsPerStepEnum,smbslices));
269
270 parameters->AddObject(iomodel->CopyConstantObject("md.smb.averaging",SmbAveragingEnum));
271
272 switch(smb_model){
273 case SMBforcingEnum:
274 break;
275 case SMBgembEnum:
276 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIdx",SmbAIdxEnum));
277 parameters->AddObject(iomodel->CopyConstantObject("md.smb.eIdx",SmbEIdxEnum));
278 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tcIdx",SmbTcIdxEnum));
279 parameters->AddObject(iomodel->CopyConstantObject("md.smb.swIdx",SmbSwIdxEnum));
280 parameters->AddObject(iomodel->CopyConstantObject("md.smb.denIdx",SmbDenIdxEnum));
281 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dsnowIdx",SmbDsnowIdxEnum));
282 parameters->AddObject(iomodel->CopyConstantObject("md.smb.cldFrac",SmbCldFracEnum));
283 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0wet",SmbT0wetEnum));
284 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0dry",SmbT0dryEnum));
285 parameters->AddObject(iomodel->CopyConstantObject("md.smb.K",SmbKEnum));
286 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aSnow",SmbASnowEnum));
287 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));
288 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dt",SmbDtEnum));
289 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isgraingrowth",SmbIsgraingrowthEnum));
290 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isalbedo",SmbIsalbedoEnum));
291 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isshortwave",SmbIsshortwaveEnum));
292 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isthermal",SmbIsthermalEnum));
293 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isaccumulation",SmbIsaccumulationEnum));
294 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismelt",SmbIsmeltEnum));
295 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdensification",SmbIsdensificationEnum));
296 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum));
297 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isconstrainsurfaceT",SmbIsconstrainsurfaceTEnum));
298 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdeltaLWup",SmbIsdeltaLWupEnum));
299 parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
300 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
301 parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));
302 parameters->AddObject(iomodel->CopyConstantObject("md.smb.teThresh",SmbTeThreshEnum));
303 break;
304 case SMBpddEnum:
305 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdelta18o",SmbIsdelta18oEnum));
306 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
307 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
308 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
309 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
310
311 if(ismungsm){
312 iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2);
313 parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,cycle,M));
314 iomodel->DeleteData(temp,"md.smb.Pfac");
315
316 iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
317 parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,cycle,M));
318 iomodel->DeleteData(temp,"md.smb.Tdiff");
319
320 iomodel->FetchData(&temp,&N,&M,"md.smb.sealev"); _assert_(N==2);
321 parameters->AddObject(new TransientParam(SmbSealevEnum,&temp[0],&temp[M],interp,cycle,M));
322 iomodel->DeleteData(temp,"md.smb.sealev");
323 }
324 if(isdelta18o){
325 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
326 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,cycle,M));
327 iomodel->DeleteData(temp,"md.smb.delta18o");
328
329 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o_surface"); _assert_(N==2);
330 parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,cycle,M));
331 iomodel->DeleteData(temp,"md.smb.delta18o_surface");
332 }
333
334 break;
335 case SMBpddSicopolisEnum:
336 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isfirnwarming",SmbIsfirnwarmingEnum));
337 break;
338 case SMBd18opddEnum:
339 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
340 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isd18opd",SmbIsd18opdEnum));
341 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
342 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
343 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
344 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
345 if(isd18opd){
346 parameters->AddObject(iomodel->CopyConstantObject("md.smb.f",SmbFEnum));
347 parameters->AddObject(iomodel->CopyConstantObject("md.smb.istemperaturescaled",SmbIstemperaturescaledEnum));
348 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipscaled",SmbIsprecipscaledEnum));
349 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
350 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,cycle,M));
351 iomodel->DeleteData(temp,"md.smb.delta18o");
352
353 IssmDouble yts;
354 bool istemperaturescaled,isprecipscaled;
355 iomodel->FindConstant(&yts,"md.constants.yts");
356 iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
357 iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
358 if(!istemperaturescaled){
359 /*Fetch array*/
360 IssmDouble* doublearray = NULL;
361 int M,N;
362 iomodel->FetchData(&doublearray,&M,&N,"md.smb.temperatures_reconstructed");
363 if(M!=iomodel->numberofvertices+1) _error_("md.smb.temperatures_reconstructed should have nbv+1 rows");
364 if(N%12!=0) _error_("md.smb.temperatures_reconstructed should have a multiple of 12 columns (since it is monthly)");
365 int numyears = N/12; _assert_(numyears*12==N);
366
367 /*Check times*/
368 #ifdef _ISSM_DEBUG_
369 for(int i=0;i<numyears;i++){
370 for(int j=1;j<12;j++){
371 //_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
372 _assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
373 }
374 }
375 #endif
376
377 /*Build time*/
378 IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
379 parameters->AddObject(new DoubleVecParam(SmbTemperaturesReconstructedYearsEnum,times,numyears));
380 xDelete<IssmDouble>(times);
381 iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
382 }
383 if(!isprecipscaled){
384 /*Fetch array*/
385 IssmDouble* doublearray = NULL;
386 int M,N;
387 iomodel->FetchData(&doublearray,&M,&N,"md.smb.precipitations_reconstructed");
388 if(M!=iomodel->numberofvertices+1) _error_("md.smb.precipitations_reconstructed should have nbv+1 rows");
389 if(N%12!=0) _error_("md.smb.precipitations_reconstructed should have a multiple of 12 columns (since it is monthly)");
390 int numyears = N/12; _assert_(numyears*12==N);
391
392 /*Check times*/
393 #ifdef _ISSM_DEBUG_
394 for(int i=0;i<numyears;i++){
395 for(int j=1;j<12;j++){
396 //_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
397 _assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
398 }
399 }
400 #endif
401
402 /*Build time*/
403 IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
404 parameters->AddObject(new DoubleVecParam(SmbPrecipitationsReconstructedYearsEnum,times,numyears));
405 xDelete<IssmDouble>(times);
406 iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
407 }
408 }
409
410 break;
411 case SMBgradientsEnum:
412 /*Nothing to add to parameters*/
413 break;
414 case SMBgradientselaEnum:
415 /*Nothing to add to parameters*/
416 break;
417 case SMBarmaEnum:
418 /*Nothing to add to parameters*/
419 break;
420 case SMBhenningEnum:
421 /*Nothing to add to parameters*/
422 break;
423 case SMBcomponentsEnum:
424 break;
425 case SMBmeltcomponentsEnum:
426 break;
427 case SMBgradientscomponentsEnum:
428 parameters->AddObject(iomodel->CopyConstantObject("md.smb.accualti",SmbAccualtiEnum));
429 parameters->AddObject(iomodel->CopyConstantObject("md.smb.runoffalti",SmbRunoffaltiEnum));
430
431 iomodel->FetchData(&temp,&N,&M,"md.smb.accugrad"); _assert_(N==2);
432 parameters->AddObject(new TransientParam(SmbAccugradEnum,&temp[0],&temp[M],interp,cycle,M));
433 iomodel->DeleteData(temp,"md.smb.accugrad");
434 iomodel->FetchData(&temp,&N,&M,"md.smb.runoffgrad"); _assert_(N==2);
435 parameters->AddObject(new TransientParam(SmbRunoffgradEnum,&temp[0],&temp[M],interp,cycle,M));
436 iomodel->DeleteData(temp,"md.smb.runoffgrad");
437
438 iomodel->FetchData(&temp,&N,&M,"md.smb.accuref"); _assert_(N==2);
439 parameters->AddObject(new TransientParam(SmbAccurefEnum,&temp[0],&temp[M],interp,cycle,M));
440 iomodel->DeleteData(temp,"md.smb.accuref");
441 iomodel->FetchData(&temp,&N,&M,"md.smb.runoffref"); _assert_(N==2);
442 parameters->AddObject(new TransientParam(SmbRunoffrefEnum,&temp[0],&temp[M],interp,cycle,M));
443 iomodel->DeleteData(temp,"md.smb.runoffref");
444 break;
445 case SMBsemicEnum:
446 int ismethod;
447 parameters->FindParam(&ismethod,SmbSemicMethodEnum);
448 if (ismethod == 1){
449 parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfacElevation",SmbDesfacElevEnum));
450 parameters->AddObject(iomodel->CopyConstantObject("md.smb.hcrit",SmbSemicHcritEnum));
451 parameters->AddObject(iomodel->CopyConstantObject("md.smb.rcrit",SmbSemicRcritEnum));
452 /*Define albedo parameters.*/
453 parameters->AddObject(iomodel->CopyConstantObject("md.smb.albedo_scheme",SmbAlbedoSchemeEnum));
454 parameters->AddObject(iomodel->CopyConstantObject("md.smb.alb_smax",SmbAlbedoSnowMaxEnum));
455 parameters->AddObject(iomodel->CopyConstantObject("md.smb.alb_smin",SmbAlbedoSnowMinEnum));
456 parameters->AddObject(iomodel->CopyConstantObject("md.smb.albi",SmbAlbedoIceEnum));
457 parameters->AddObject(iomodel->CopyConstantObject("md.smb.albl",SmbAlbedoLandEnum));
458
459 //albedo parameter - slatter
460 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tmin",SmbSemicTmaxEnum));
461 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tmax",SmbSemicTminEnum));
462
463 //albedo parameter - isba & denby
464 parameters->AddObject(iomodel->CopyConstantObject("md.smb.mcrit",SmbSemicMcritEnum));
465 //albedo parameter - isba
466 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tau_a",SmbSemicTauAEnum));
467 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tau_f",SmbSemicTauFEnum));
468 parameters->AddObject(iomodel->CopyConstantObject("md.smb.wcrit",SmbSemicWcritEnum));
469 //albedo parameter - alex
470 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tmid",SmbSemicTmidEnum));
471 parameters->AddObject(iomodel->CopyConstantObject("md.smb.afac",SmbSemicAfacEnum));
472 }
473 /*Nothing to add to parameters*/
474 break;
475 case SMBdebrisMLEnum:
476 /*Nothing to add to parameters*/
477 break;
478 default:
479 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
480 }
481
482 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.smb.requested_outputs");
483 parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
484 if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
485 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.smb.requested_outputs");
486
487}/*}}}*/
488
489/*Finite Element Analysis*/
490void SmbAnalysis::Core(FemModel* femmodel){/*{{{*/
491
492 int smb_model;
493
494 /*Figure out smb model: */
495 femmodel->parameters->FindParam(&smb_model,SmbEnum);
496
497 /*branch to correct module*/
498 switch(smb_model){
499 case SMBforcingEnum:
500 SmbForcingx(femmodel);
501 break;
502 case SMBgembEnum:
503 Gembx(femmodel);
504 break;
505 case SMBpddEnum:
506 bool isdelta18o,ismungsm;
507 femmodel->parameters->FindParam(&isdelta18o,SmbIsdelta18oEnum);
508 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
509 if(isdelta18o){
510 if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n");
511 Delta18oParameterizationx(femmodel);
512 }
513 if(ismungsm){
514 if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n");
515 MungsmtpParameterizationx(femmodel);
516 }
517 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
518 PositiveDegreeDayx(femmodel);
519 break;
520 case SMBpddSicopolisEnum:
521 if(VerboseSolution()) _printf0_(" call SICOPOLIS positive degree day module\n");
522 PositiveDegreeDaySicopolisx(femmodel);
523 break;
524 case SMBd18opddEnum:
525 bool isd18opd;
526 femmodel->parameters->FindParam(&isd18opd,SmbIsd18opdEnum);
527 if(isd18opd){
528 if(VerboseSolution()) _printf0_(" call Delta18opdParameterization module\n");
529 Delta18opdParameterizationx(femmodel);
530 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
531 PositiveDegreeDayx(femmodel);
532 }
533 break;
534 case SMBgradientsEnum:
535 if(VerboseSolution())_printf0_(" call smb gradients module\n");
536 SmbGradientsx(femmodel);
537 break;
538 case SMBgradientselaEnum:
539 if(VerboseSolution())_printf0_(" call smb gradients ela module\n");
540 SmbGradientsElax(femmodel);
541 break;
542 case SMBarmaEnum:
543 if(VerboseSolution())_printf0_(" call smb arma module\n");
544 Smbarmax(femmodel);
545 break;
546 case SMBhenningEnum:
547 if(VerboseSolution())_printf0_(" call smb Henning module\n");
548 SmbHenningx(femmodel);
549 break;
550 case SMBcomponentsEnum:
551 if(VerboseSolution())_printf0_(" call smb Components module\n");
552 SmbComponentsx(femmodel);
553 break;
554 case SMBmeltcomponentsEnum:
555 if(VerboseSolution())_printf0_(" call smb Melt Components module\n");
556 SmbMeltComponentsx(femmodel);
557 break;
558 case SMBgcmEnum:
559 /*Nothing to be done*/
560 break;
561 case SMBgradientscomponentsEnum:
562 if(VerboseSolution())_printf0_(" call smb gradients components module\n");
563 SmbGradientsComponentsx(femmodel);
564 break;
565 case SMBsemicEnum:
566 #ifdef _HAVE_SEMIC_
567 if(VerboseSolution())_printf0_(" call smb SEMIC module\n");
568 int ismethod;
569 femmodel->parameters->FindParam(&ismethod,SmbSemicMethodEnum);
570 SmbSemicx(femmodel,ismethod);
571 #else
572 _error_("SEMIC not installed");
573 #endif //_HAVE_SEMIC_
574 break;
575 case SMBdebrisMLEnum:
576 if(VerboseSolution())_printf0_(" call smb debris Mayer & Liculli module\n");
577 SmbDebrisMLx(femmodel);
578 break;
579 default:
580 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
581 }
582
583}/*}}}*/
584void SmbAnalysis::PreCore(FemModel* femmodel){/*{{{*/
585 _error_("not implemented");
586}/*}}}*/
587ElementVector* SmbAnalysis::CreateDVector(Element* element){/*{{{*/
588 _error_("not implemented");
589}/*}}}*/
590ElementMatrix* SmbAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
591_error_("Not implemented");
592}/*}}}*/
593ElementMatrix* SmbAnalysis::CreateKMatrix(Element* element){/*{{{*/
594 _error_("not implemented yet");
595}/*}}}*/
596ElementVector* SmbAnalysis::CreatePVector(Element* element){/*{{{*/
597_error_("not implemented yet");
598}/*}}}*/
599void SmbAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
600 _error_("not implemented yet");
601}/*}}}*/
602void SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
603 _error_("Not implemented yet");
604}/*}}}*/
605void SmbAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
606 _error_("not implemented yet");
607}/*}}}*/
608void SmbAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
609 /*Default, do nothing*/
610 return;
611}/*}}}*/
Note: See TracBrowser for help on using the repository browser.