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

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

CHG: add options for del018 pdd to take in reconstructed climate
File size: 15.9 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/*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;
23 bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled;
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
35 /*Figure out smb model: */
36 iomodel->FindConstant(&smb_model,"md.smb.model");
37
38 switch(smb_model){
39 case SMBforcingEnum:
40 iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
41 break;
42 case SMBgembEnum:
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);
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);
71 iomodel->FetchDataToInput(elements,"md.smb.aValue",SmbAValueEnum);
72 iomodel->FetchDataToInput(elements,"md.smb.teValue",SmbTeValueEnum);
73 break;
74 case SMBpddEnum:
75 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
76 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
77 iomodel->FetchDataToInput(elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
78 iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum);
79 iomodel->FetchDataToInput(elements,"md.smb.s0t",SmbS0tEnum);
80 if(isdelta18o || ismungsm){
81 iomodel->FetchDataToInput(elements,"md.smb.temperatures_lgm",SmbTemperaturesLgmEnum);
82 iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
83 iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
84 iomodel->FetchDataToInput(elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum);
85 }else{
86 iomodel->FetchDataToInput(elements,"md.smb.precipitation",SmbPrecipitationEnum);
87 iomodel->FetchDataToInput(elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
88 }
89 break;
90 case SMBd18opddEnum:
91 iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
92 iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
93 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
94 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
95 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
96 iomodel->FetchDataToInput(elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
97 iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum);
98 iomodel->FetchDataToInput(elements,"md.smb.s0t",SmbS0tEnum);
99 if(isd18opd){
100 iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
101 iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
102 if(!istemperaturescaled){
103 iomodel->FetchDataToInput(elements,"md.smb.temperatures_reconstructed",SmbTemperaturesReconstructedEnum);
104 }
105 if(!isprecipscaled){
106 iomodel->FetchDataToInput(elements,"md.smb.precipitations_reconstructed",SmbPrecipitationsReconstructedEnum);
107 }
108 }
109 if(issetpddfac){
110 iomodel->FetchDataToInput(elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
111 iomodel->FetchDataToInput(elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
112 }
113 break;
114 case SMBgradientsEnum:
115 iomodel->FetchDataToInput(elements,"md.smb.href",SmbHrefEnum);
116 iomodel->FetchDataToInput(elements,"md.smb.smbref",SmbSmbrefEnum);
117 iomodel->FetchDataToInput(elements,"md.smb.b_pos",SmbBPosEnum);
118 iomodel->FetchDataToInput(elements,"md.smb.b_neg",SmbBNegEnum);
119 break;
120 case SMBgradientselaEnum:
121 iomodel->FetchDataToInput(elements,"md.smb.ela",SmbElaEnum);
122 iomodel->FetchDataToInput(elements,"md.smb.b_pos",SmbBPosEnum);
123 iomodel->FetchDataToInput(elements,"md.smb.b_neg",SmbBNegEnum);
124 iomodel->FetchDataToInput(elements,"md.smb.b_max",SmbBMaxEnum);
125 iomodel->FetchDataToInput(elements,"md.smb.b_min",SmbBMinEnum);
126 break;
127 case SMBhenningEnum:
128 iomodel->FetchDataToInput(elements,"md.smb.smbref",SmbSmbrefEnum,0.);
129 break;
130 case SMBcomponentsEnum:
131 iomodel->FetchDataToInput(elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
132 iomodel->FetchDataToInput(elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
133 iomodel->FetchDataToInput(elements,"md.smb.runoff",SmbRunoffEnum,0.);
134 break;
135 case SMBmeltcomponentsEnum:
136 iomodel->FetchDataToInput(elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
137 iomodel->FetchDataToInput(elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
138 iomodel->FetchDataToInput(elements,"md.smb.melt",SmbMeltEnum,0.);
139 iomodel->FetchDataToInput(elements,"md.smb.refreeze",SmbRefreezeEnum,0.);
140 break;
141 default:
142 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
143 }
144
145
146
147}/*}}}*/
148void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
149
150 int numoutputs;
151 char** requestedoutputs = NULL;
152 bool isdelta18o,ismungsm,isd18opd,issetpddfac,interp;
153 int smb_model;
154 IssmDouble *temp = NULL;
155 int N,M;
156
157 parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
158
159 iomodel->FindConstant(&smb_model,"md.smb.model");
160 iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
161
162 switch(smb_model){
163 case SMBforcingEnum:
164 /*Nothing to add to parameters*/
165 break;
166 case SMBgembEnum:
167 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIdx",SmbAIdxEnum));
168 parameters->AddObject(iomodel->CopyConstantObject("md.smb.swIdx",SmbSwIdxEnum));
169 parameters->AddObject(iomodel->CopyConstantObject("md.smb.denIdx",SmbDenIdxEnum));
170 parameters->AddObject(iomodel->CopyConstantObject("md.smb.cldFrac",SmbCldFracEnum));
171 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0wet",SmbT0wetEnum));
172 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0dry",SmbT0dryEnum));
173 parameters->AddObject(iomodel->CopyConstantObject("md.smb.K",SmbKEnum));
174 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aSnow",SmbASnowEnum));
175 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));
176 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dt",SmbDtEnum));
177 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isgraingrowth",SmbIsgraingrowthEnum));
178 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isalbedo",SmbIsalbedoEnum));
179 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isshortwave",SmbIsshortwaveEnum));
180 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isthermal",SmbIsthermalEnum));
181 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isaccumulation",SmbIsaccumulationEnum));
182 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismelt",SmbIsmeltEnum));
183 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdensification",SmbIsdensificationEnum));
184 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum));
185 parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
186 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
187 parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));
188 break;
189 case SMBpddEnum:
190 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdelta18o",SmbIsdelta18oEnum));
191 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
192 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
193 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
194 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
195
196 if(ismungsm){
197 iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2);
198 parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,M));
199 iomodel->DeleteData(temp,"md.smb.Pfac");
200
201 iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
202 parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,M));
203 iomodel->DeleteData(temp,"md.smb.Tdiff");
204
205 iomodel->FetchData(&temp,&N,&M,"md.smb.sealev"); _assert_(N==2);
206 parameters->AddObject(new TransientParam(SmbSealevEnum,&temp[0],&temp[M],interp,M));
207 iomodel->DeleteData(temp,"md.smb.sealev");
208 }
209 if(isdelta18o){
210 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
211 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
212 iomodel->DeleteData(temp,"md.smb.delta18o");
213
214 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o_surface"); _assert_(N==2);
215 parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,M));
216 iomodel->DeleteData(temp,"md.smb.delta18o_surface");
217 }
218 break;
219 case SMBd18opddEnum:
220 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
221 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isd18opd",SmbIsd18opdEnum));
222 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
223 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
224 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
225 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
226 if(isd18opd){
227 parameters->AddObject(iomodel->CopyConstantObject("md.smb.f",SmbFEnum));
228 parameters->AddObject(iomodel->CopyConstantObject("md.smb.istemperaturescaled",SmbIstemperaturescaledEnum));
229 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipscaled",SmbIsprecipscaledEnum));
230 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
231 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
232 iomodel->DeleteData(temp,"md.smb.delta18o");
233 }
234 break;
235 case SMBgradientsEnum:
236 /*Nothing to add to parameters*/
237 break;
238 case SMBgradientselaEnum:
239 /*Nothing to add to parameters*/
240 break;
241 case SMBhenningEnum:
242 /*Nothing to add to parameters*/
243 break;
244 case SMBcomponentsEnum:
245 /*Nothing to add to parameters*/
246 break;
247 case SMBmeltcomponentsEnum:
248 /*Nothing to add to parameters*/
249 break;
250 default:
251 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
252 }
253
254 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.smb.requested_outputs");
255 parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
256 if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
257 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.smb.requested_outputs");
258
259}/*}}}*/
260
261/*Finite Element Analysis*/
262void SmbAnalysis::Core(FemModel* femmodel){/*{{{*/
263
264 int smb_model;
265
266 /*Figure out smb model: */
267 femmodel->parameters->FindParam(&smb_model,SmbEnum);
268
269 /*branch to correct module*/
270 switch(smb_model){
271 case SMBforcingEnum:
272 /*Nothing to be done*/
273 break;
274 case SMBgembEnum:
275 Gembx(femmodel);
276 break;
277 case SMBpddEnum:
278 bool isdelta18o,ismungsm;
279 femmodel->parameters->FindParam(&isdelta18o,SmbIsdelta18oEnum);
280 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
281 if(isdelta18o){
282 if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n");
283 Delta18oParameterizationx(femmodel);
284 }
285 if(ismungsm){
286 if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n");
287 MungsmtpParameterizationx(femmodel);
288 }
289 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
290 PositiveDegreeDayx(femmodel);
291 break;
292 case SMBd18opddEnum:
293 bool isd18opd;
294 femmodel->parameters->FindParam(&isd18opd,SmbIsd18opdEnum);
295 if(isd18opd){
296 if(VerboseSolution()) _printf0_(" call Delta18opdParameterization module\n");
297 Delta18opdParameterizationx(femmodel);
298 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
299 PositiveDegreeDayx(femmodel);
300 }
301 break;
302 case SMBgradientsEnum:
303 if(VerboseSolution())_printf0_(" call smb gradients module\n");
304 SmbGradientsx(femmodel);
305 break;
306 case SMBgradientselaEnum:
307 if(VerboseSolution())_printf0_(" call smb gradients ela module\n");
308 SmbGradientsElax(femmodel);
309 break;
310 case SMBhenningEnum:
311 if(VerboseSolution())_printf0_(" call smb Henning module\n");
312 SmbHenningx(femmodel);
313 break;
314 case SMBcomponentsEnum:
315 if(VerboseSolution())_printf0_(" call smb Components module\n");
316 SmbComponentsx(femmodel);
317 break;
318 case SMBmeltcomponentsEnum:
319 if(VerboseSolution())_printf0_(" call smb Melt Components module\n");
320 SmbMeltComponentsx(femmodel);
321 break;
322 case SMBgcmEnum:
323 /*Nothing to be done*/
324 break;
325 default:
326 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
327 }
328
329}/*}}}*/
330ElementVector* SmbAnalysis::CreateDVector(Element* element){/*{{{*/
331 _error_("not implemented");
332}/*}}}*/
333ElementMatrix* SmbAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
334_error_("Not implemented");
335}/*}}}*/
336ElementMatrix* SmbAnalysis::CreateKMatrix(Element* element){/*{{{*/
337 _error_("not implemented yet");
338}/*}}}*/
339ElementVector* SmbAnalysis::CreatePVector(Element* element){/*{{{*/
340_error_("not implemented yet");
341}/*}}}*/
342void SmbAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
343 _error_("not implemented yet");
344}/*}}}*/
345void SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
346 _error_("Not implemented yet");
347}/*}}}*/
348void SmbAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
349 _error_("not implemented yet");
350}/*}}}*/
351void SmbAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
352 /*Default, do nothing*/
353 return;
354}/*}}}*/
Note: See TracBrowser for help on using the repository browser.