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

Last change on this file since 21232 was 21232, checked in by langchar, 8 years ago

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