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

Last change on this file since 22448 was 22448, checked in by schlegel, 7 years ago

CHG: add pdd factor input field
File size: 15.6 KB
RevLine 
[19528]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/*Model processing*/
8void SmbAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
9 /*No constraints*/
10}/*}}}*/
11void SmbAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
12 /*No loads*/
13}/*}}}*/
14void SmbAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
15 ::CreateNodes(nodes,iomodel,SmbAnalysisEnum,P1Enum);
16}/*}}}*/
17int SmbAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
18 return 1;
19}/*}}}*/
20void SmbAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
21
22 int smb_model;
[22448]23 bool isdelta18o,ismungsm,isd18opd,issetpddfac;
[19528]24
25 /*Update elements: */
26 int counter=0;
27 for(int i=0;i<iomodel->numberofelements;i++){
28 if(iomodel->my_elements[i]){
29 Element* element=(Element*)elements->GetObjectByOffset(counter);
30 element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
31 counter++;
32 }
33 }
34
[19554]35 /*Figure out smb model: */
[20690]36 iomodel->FindConstant(&smb_model,"md.smb.model");
[19554]37
[19528]38 switch(smb_model){
39 case SMBforcingEnum:
[20690]40 iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
[19528]41 break;
[19554]42 case SMBgembEnum:
[20690]43 iomodel->FetchDataToInput(elements,"md.smb.Ta",SmbTaEnum);
44 iomodel->FetchDataToInput(elements,"md.smb.V",SmbVEnum);
45 iomodel->FetchDataToInput(elements,"md.smb.dswrf",SmbDswrfEnum);
46 iomodel->FetchDataToInput(elements,"md.smb.dlwrf",SmbDlwrfEnum);
47 iomodel->FetchDataToInput(elements,"md.smb.P",SmbPEnum);
48 iomodel->FetchDataToInput(elements,"md.smb.eAir",SmbEAirEnum);
49 iomodel->FetchDataToInput(elements,"md.smb.pAir",SmbPAirEnum);
50 iomodel->FetchDataToInput(elements,"md.smb.zTop",SmbZTopEnum);
51 iomodel->FetchDataToInput(elements,"md.smb.dzTop",SmbDzTopEnum);
52 iomodel->FetchDataToInput(elements,"md.smb.dzMin",SmbDzMinEnum);
53 iomodel->FetchDataToInput(elements,"md.smb.zY",SmbZYEnum);
54 iomodel->FetchDataToInput(elements,"md.smb.zMax",SmbZMaxEnum);
55 iomodel->FetchDataToInput(elements,"md.smb.zMin",SmbZMinEnum);
56 iomodel->FetchDataToInput(elements,"md.smb.Tmean",SmbTmeanEnum);
57 iomodel->FetchDataToInput(elements,"md.smb.C",SmbCEnum);
58 iomodel->FetchDataToInput(elements,"md.smb.Tz",SmbTzEnum);
59 iomodel->FetchDataToInput(elements,"md.smb.Vz",SmbVzEnum);
[21232]60 InputUpdateFromConstantx(elements,0.,SmbIsInitializedEnum);
61 iomodel->FetchDataToInput(elements,"md.smb.Dzini",SmbDziniEnum);
62 iomodel->FetchDataToInput(elements,"md.smb.Dini",SmbDiniEnum);
63 iomodel->FetchDataToInput(elements,"md.smb.Reini",SmbReiniEnum);
64 iomodel->FetchDataToInput(elements,"md.smb.Gdnini",SmbGdniniEnum);
65 iomodel->FetchDataToInput(elements,"md.smb.Gspini",SmbGspiniEnum);
66 iomodel->FetchDataToInput(elements,"md.smb.ECini",SmbECiniEnum);
67 iomodel->FetchDataToInput(elements,"md.smb.Wini",SmbWiniEnum);
68 iomodel->FetchDataToInput(elements,"md.smb.Aini",SmbAiniEnum);
69 iomodel->FetchDataToInput(elements,"md.smb.Tini",SmbTiniEnum);
70 iomodel->FetchDataToInput(elements,"md.smb.Sizeini",SmbSizeiniEnum);
[19554]71 break;
[19528]72 case SMBpddEnum:
[20690]73 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
74 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
75 iomodel->FetchDataToInput(elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
76 iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum);
77 iomodel->FetchDataToInput(elements,"md.smb.s0t",SmbS0tEnum);
[19528]78 if(isdelta18o || ismungsm){
[20690]79 iomodel->FetchDataToInput(elements,"md.smb.temperatures_lgm",SmbTemperaturesLgmEnum);
80 iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
81 iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
82 iomodel->FetchDataToInput(elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum);
[19528]83 }
84 else{
[20690]85 iomodel->FetchDataToInput(elements,"md.smb.precipitation",SmbPrecipitationEnum);
86 iomodel->FetchDataToInput(elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
[19528]87 }
88 break;
89 case SMBd18opddEnum:
[20690]90 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
91 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
[22448]92 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
[20690]93 iomodel->FetchDataToInput(elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
94 iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum);
95 iomodel->FetchDataToInput(elements,"md.smb.s0t",SmbS0tEnum);
96 if(isd18opd){
97 iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
98 iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
[19528]99 }
[22448]100 if(issetpddfac){
101 iomodel->FetchDataToInput(elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
102 iomodel->FetchDataToInput(elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
103 }
[19528]104 break;
105 case SMBgradientsEnum:
[20690]106 iomodel->FetchDataToInput(elements,"md.smb.href",SmbHrefEnum);
107 iomodel->FetchDataToInput(elements,"md.smb.smbref",SmbSmbrefEnum);
108 iomodel->FetchDataToInput(elements,"md.smb.b_pos",SmbBPosEnum);
109 iomodel->FetchDataToInput(elements,"md.smb.b_neg",SmbBNegEnum);
[19528]110 break;
[21469]111 case SMBgradientselaEnum:
112 iomodel->FetchDataToInput(elements,"md.smb.ela",SmbElaEnum);
113 iomodel->FetchDataToInput(elements,"md.smb.b_pos",SmbBPosEnum);
114 iomodel->FetchDataToInput(elements,"md.smb.b_neg",SmbBNegEnum);
115 iomodel->FetchDataToInput(elements,"md.smb.b_max",SmbBMaxEnum);
116 iomodel->FetchDataToInput(elements,"md.smb.b_min",SmbBMinEnum);
117 break;
[19528]118 case SMBhenningEnum:
[20690]119 iomodel->FetchDataToInput(elements,"md.smb.smbref",SmbSmbrefEnum,0.);
[19528]120 break;
121 case SMBcomponentsEnum:
[20690]122 iomodel->FetchDataToInput(elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
123 iomodel->FetchDataToInput(elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
124 iomodel->FetchDataToInput(elements,"md.smb.runoff",SmbRunoffEnum,0.);
[19528]125 break;
126 case SMBmeltcomponentsEnum:
[20690]127 iomodel->FetchDataToInput(elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
128 iomodel->FetchDataToInput(elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
129 iomodel->FetchDataToInput(elements,"md.smb.melt",SmbMeltEnum,0.);
130 iomodel->FetchDataToInput(elements,"md.smb.refreeze",SmbRefreezeEnum,0.);
[19528]131 break;
132 default:
133 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
134 }
[19554]135
136
[19528]137
138}/*}}}*/
139void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
140
141 int numoutputs;
142 char** requestedoutputs = NULL;
[22448]143 bool isdelta18o,ismungsm,isd18opd,issetpddfac,interp;
[19528]144 int smb_model;
145 IssmDouble *temp = NULL;
146 int N,M;
147
[20690]148 parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
[19528]149
[20690]150 iomodel->FindConstant(&smb_model,"md.smb.model");
151 iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
[19528]152
153 switch(smb_model){
154 case SMBforcingEnum:
155 /*Nothing to add to parameters*/
156 break;
[19554]157 case SMBgembEnum:
[20690]158 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIdx",SmbAIdxEnum));
159 parameters->AddObject(iomodel->CopyConstantObject("md.smb.swIdx",SmbSwIdxEnum));
160 parameters->AddObject(iomodel->CopyConstantObject("md.smb.denIdx",SmbDenIdxEnum));
161 parameters->AddObject(iomodel->CopyConstantObject("md.smb.outputFreq",SmbOutputFreqEnum));
162 parameters->AddObject(iomodel->CopyConstantObject("md.smb.cldFrac",SmbCldFracEnum));
163 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0wet",SmbT0wetEnum));
164 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0dry",SmbT0dryEnum));
165 parameters->AddObject(iomodel->CopyConstantObject("md.smb.K",SmbKEnum));
166 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aSnow",SmbASnowEnum));
167 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));
168 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dt",SmbDtEnum));
169 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isgraingrowth",SmbIsgraingrowthEnum));
170 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isalbedo",SmbIsalbedoEnum));
171 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isshortwave",SmbIsshortwaveEnum));
172 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isthermal",SmbIsthermalEnum));
173 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isaccumulation",SmbIsaccumulationEnum));
174 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismelt",SmbIsmeltEnum));
175 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdensification",SmbIsdensificationEnum));
176 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum));
177 parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
[19554]178 break;
[19528]179 case SMBpddEnum:
[20690]180 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdelta18o",SmbIsdelta18oEnum));
181 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
[22448]182 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
[20690]183 parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
184 parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
185 parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlapslgm",SmbRlapslgmEnum));
186 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
187 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
[19528]188
189 if(ismungsm){
[20690]190 iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2);
[19528]191 parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,M));
[20690]192 iomodel->DeleteData(temp,"md.smb.Pfac");
[19528]193
[20690]194 iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
[19528]195 parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,M));
[20690]196 iomodel->DeleteData(temp,"md.smb.Tdiff");
[19528]197
[20690]198 iomodel->FetchData(&temp,&N,&M,"md.smb.sealev"); _assert_(N==2);
[19528]199 parameters->AddObject(new TransientParam(SmbSealevEnum,&temp[0],&temp[M],interp,M));
[20690]200 iomodel->DeleteData(temp,"md.smb.sealev");
[19528]201 }
202 if(isdelta18o){
[20690]203 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
[19528]204 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
[20690]205 iomodel->DeleteData(temp,"md.smb.delta18o");
[19528]206
[20690]207 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o_surface"); _assert_(N==2);
[19528]208 parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,M));
[20690]209 iomodel->DeleteData(temp,"md.smb.delta18o_surface");
[19528]210 }
211 break;
212 case SMBd18opddEnum:
[20690]213 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
214 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isd18opd",SmbIsd18opdEnum));
[22448]215 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
[20690]216 parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
217 parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
218 parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlapslgm",SmbRlapslgmEnum));
219 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
220 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
[22448]221 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
[19528]222 if(isd18opd){
[20690]223 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
[19528]224 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
[20690]225 iomodel->DeleteData(temp,"md.smb.delta18o");
[19528]226
[20690]227 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dpermil",SmbDpermilEnum));
[21709]228 parameters->AddObject(iomodel->CopyConstantObject("md.smb.f",SmbFEnum));
[19528]229 }
230 break;
231 case SMBgradientsEnum:
232 /*Nothing to add to parameters*/
233 break;
[21469]234 case SMBgradientselaEnum:
235 /*Nothing to add to parameters*/
236 break;
[19528]237 case SMBhenningEnum:
238 /*Nothing to add to parameters*/
239 break;
240 case SMBcomponentsEnum:
241 /*Nothing to add to parameters*/
242 break;
243 case SMBmeltcomponentsEnum:
244 /*Nothing to add to parameters*/
245 break;
246 default:
247 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
248 }
249
[20690]250 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.smb.requested_outputs");
[19528]251 parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
252 if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
[20690]253 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.smb.requested_outputs");
[19528]254
255}/*}}}*/
256
257/*Finite Element Analysis*/
258void SmbAnalysis::Core(FemModel* femmodel){/*{{{*/
259
260 int smb_model;
261
262 /*Figure out smb model: */
263 femmodel->parameters->FindParam(&smb_model,SmbEnum);
264
265 /*branch to correct module*/
266 switch(smb_model){
267 case SMBforcingEnum:
268 /*Nothing to be done*/
269 break;
[19554]270 case SMBgembEnum:
271 Gembx(femmodel);
272 break;
[19528]273 case SMBpddEnum:
274 bool isdelta18o,ismungsm;
275 femmodel->parameters->FindParam(&isdelta18o,SmbIsdelta18oEnum);
276 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
277 if(isdelta18o){
278 if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n");
279 Delta18oParameterizationx(femmodel);
280 }
281 if(ismungsm){
282 if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n");
283 MungsmtpParameterizationx(femmodel);
284 }
285 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
286 PositiveDegreeDayx(femmodel);
287 break;
288 case SMBd18opddEnum:
289 bool isd18opd;
290 femmodel->parameters->FindParam(&isd18opd,SmbIsd18opdEnum);
291 if(isd18opd){
292 if(VerboseSolution()) _printf0_(" call Delta18opdParameterization module\n");
293 Delta18opdParameterizationx(femmodel);
294 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
295 PositiveDegreeDayx(femmodel);
296 }
297 break;
298 case SMBgradientsEnum:
299 if(VerboseSolution())_printf0_(" call smb gradients module\n");
300 SmbGradientsx(femmodel);
301 break;
[21469]302 case SMBgradientselaEnum:
303 if(VerboseSolution())_printf0_(" call smb gradients ela module\n");
304 SmbGradientsElax(femmodel);
305 break;
[19528]306 case SMBhenningEnum:
307 if(VerboseSolution())_printf0_(" call smb Henning module\n");
308 SmbHenningx(femmodel);
309 break;
310 case SMBcomponentsEnum:
311 if(VerboseSolution())_printf0_(" call smb Components module\n");
312 SmbComponentsx(femmodel);
313 break;
314 case SMBmeltcomponentsEnum:
315 if(VerboseSolution())_printf0_(" call smb Melt Components module\n");
316 SmbMeltComponentsx(femmodel);
317 break;
318 case SMBgcmEnum:
319 /*Nothing to be done*/
320 break;
321 default:
322 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
323 }
324
325}/*}}}*/
326ElementVector* SmbAnalysis::CreateDVector(Element* element){/*{{{*/
327 _error_("not implemented");
328}/*}}}*/
329ElementMatrix* SmbAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
330_error_("Not implemented");
331}/*}}}*/
332ElementMatrix* SmbAnalysis::CreateKMatrix(Element* element){/*{{{*/
333 _error_("not implemented yet");
334}/*}}}*/
335ElementVector* SmbAnalysis::CreatePVector(Element* element){/*{{{*/
336_error_("not implemented yet");
337}/*}}}*/
338void SmbAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
339 _error_("not implemented yet");
340}/*}}}*/
341void SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
342 _error_("Not implemented yet");
343}/*}}}*/
344void SmbAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
345 _error_("not implemented yet");
346}/*}}}*/
347void SmbAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
348 /*Default, do nothing*/
349 return;
350}/*}}}*/
Note: See TracBrowser for help on using the repository browser.