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

Last change on this file since 21216 was 21216, checked in by langchar, 9 years ago

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