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

Last change on this file since 19644 was 19644, checked in by Mathieu Morlighem, 9 years ago

CHG: cleaning up unneeded code
File size: 12.2 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;
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
[19554]35 /*Figure out smb model: */
36 iomodel->Constant(&smb_model,SmbEnum);
37
38
[19528]39 switch(smb_model){
40 case SMBforcingEnum:
41 iomodel->FetchDataToInput(elements,SmbMassBalanceEnum,0.);
42 break;
[19554]43 case SMBgembEnum:
44 iomodel->FetchDataToInput(elements,SmbTaEnum);
45 iomodel->FetchDataToInput(elements,SmbVEnum);
46 iomodel->FetchDataToInput(elements,SmbDswrfEnum);
47 iomodel->FetchDataToInput(elements,SmbDlwrfEnum);
48 iomodel->FetchDataToInput(elements,SmbPEnum);
49 iomodel->FetchDataToInput(elements,SmbEAirEnum);
50 iomodel->FetchDataToInput(elements,SmbPAirEnum);
51 iomodel->FetchDataToInput(elements,SmbZTopEnum);
52 iomodel->FetchDataToInput(elements,SmbDzTopEnum);
53 iomodel->FetchDataToInput(elements,SmbDzMinEnum);
54 iomodel->FetchDataToInput(elements,SmbZYEnum);
55 iomodel->FetchDataToInput(elements,SmbZMaxEnum);
56 iomodel->FetchDataToInput(elements,SmbZMinEnum);
57 iomodel->FetchDataToInput(elements,SmbTmeanEnum);
58 iomodel->FetchDataToInput(elements,SmbCEnum);
59 iomodel->FetchDataToInput(elements,SmbTzEnum);
60 iomodel->FetchDataToInput(elements,SmbVzEnum);
61 break;
[19528]62 case SMBpddEnum:
63 iomodel->Constant(&isdelta18o,SmbIsdelta18oEnum);
64 iomodel->Constant(&ismungsm,SmbIsmungsmEnum);
65 iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
66 iomodel->FetchDataToInput(elements,SmbS0pEnum);
67 iomodel->FetchDataToInput(elements,SmbS0tEnum);
68 if(isdelta18o || ismungsm){
69 iomodel->FetchDataToInput(elements,SmbTemperaturesLgmEnum);
70 iomodel->FetchDataToInput(elements,SmbTemperaturesPresentdayEnum);
71 iomodel->FetchDataToInput(elements,SmbPrecipitationsPresentdayEnum);
72 iomodel->FetchDataToInput(elements,SmbPrecipitationsLgmEnum);
73 }
74 else{
75 iomodel->FetchDataToInput(elements,SmbPrecipitationEnum);
76 iomodel->FetchDataToInput(elements,SmbMonthlytemperaturesEnum);
77 }
78 break;
79 case SMBd18opddEnum:
80 iomodel->Constant(&ismungsm,SmbIsmungsmEnum);
81 iomodel->Constant(&isd18opd,SmbIsd18opdEnum);
82 iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
83 iomodel->FetchDataToInput(elements,SmbS0pEnum);
84 iomodel->FetchDataToInput(elements,SmbS0tEnum);
85 if (isd18opd){
86 iomodel->FetchDataToInput(elements,SmbTemperaturesPresentdayEnum);
87 iomodel->FetchDataToInput(elements,SmbPrecipitationsPresentdayEnum);
88 }
89
90 break;
91 case SMBgradientsEnum:
92 iomodel->FetchDataToInput(elements,SmbHrefEnum);
93 iomodel->FetchDataToInput(elements,SmbSmbrefEnum);
94 iomodel->FetchDataToInput(elements,SmbBPosEnum);
95 iomodel->FetchDataToInput(elements,SmbBNegEnum);
96 break;
97 case SMBhenningEnum:
98 iomodel->FetchDataToInput(elements,SmbSmbrefEnum,0.);
99 break;
100 case SMBcomponentsEnum:
101 iomodel->FetchDataToInput(elements,SmbAccumulationEnum,0.);
102 iomodel->FetchDataToInput(elements,SmbEvaporationEnum,0.);
103 iomodel->FetchDataToInput(elements,SmbRunoffEnum,0.);
104 break;
105 case SMBmeltcomponentsEnum:
106 iomodel->FetchDataToInput(elements,SmbAccumulationEnum,0.);
107 iomodel->FetchDataToInput(elements,SmbEvaporationEnum,0.);
108 iomodel->FetchDataToInput(elements,SmbMeltEnum,0.);
109 iomodel->FetchDataToInput(elements,SmbRefreezeEnum,0.);
110 break;
111 default:
112 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
113 }
[19554]114
115
[19528]116
117}/*}}}*/
118void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
119
120 int numoutputs;
121 char** requestedoutputs = NULL;
122 bool isdelta18o,ismungsm,isd18opd,interp;
123 int smb_model;
124 IssmDouble *temp = NULL;
125 int N,M;
126
127 parameters->AddObject(iomodel->CopyConstantObject(SmbEnum));
128
129 iomodel->Constant(&smb_model,SmbEnum);
130 iomodel->Constant(&interp,TimesteppingInterpForcingsEnum);
131
132 switch(smb_model){
133 case SMBforcingEnum:
134 /*Nothing to add to parameters*/
135 break;
[19554]136 case SMBgembEnum:
137 parameters->AddObject(iomodel->CopyConstantObject(SmbAIdxEnum));
138 parameters->AddObject(iomodel->CopyConstantObject(SmbSwIdxEnum));
139 parameters->AddObject(iomodel->CopyConstantObject(SmbDenIdxEnum));
140 parameters->AddObject(iomodel->CopyConstantObject(SmbOutputFreqEnum));
141 parameters->AddObject(iomodel->CopyConstantObject(SmbCldFracEnum));
142 parameters->AddObject(iomodel->CopyConstantObject(SmbT0wetEnum));
143 parameters->AddObject(iomodel->CopyConstantObject(SmbT0dryEnum));
144 parameters->AddObject(iomodel->CopyConstantObject(SmbKEnum));
145 parameters->AddObject(iomodel->CopyConstantObject(SmbASnowEnum));
146 parameters->AddObject(iomodel->CopyConstantObject(SmbAIceEnum));
147 parameters->AddObject(iomodel->CopyConstantObject(SmbDtEnum));
148 parameters->AddObject(iomodel->CopyConstantObject(SmbIsgraingrowthEnum));
149 parameters->AddObject(iomodel->CopyConstantObject(SmbIsalbedoEnum));
150 parameters->AddObject(iomodel->CopyConstantObject(SmbIsshortwaveEnum));
151 parameters->AddObject(iomodel->CopyConstantObject(SmbIsthermalEnum));
152 parameters->AddObject(iomodel->CopyConstantObject(SmbIsaccumulationEnum));
153 parameters->AddObject(iomodel->CopyConstantObject(SmbIsmeltEnum));
154 parameters->AddObject(iomodel->CopyConstantObject(SmbIsdensificationEnum));
155 parameters->AddObject(iomodel->CopyConstantObject(SmbIsturbulentfluxEnum));
[19565]156 parameters->AddObject(iomodel->CopyConstantObject(SmbInitDensityScalingEnum));
[19554]157 break;
[19528]158 case SMBpddEnum:
159 parameters->AddObject(iomodel->CopyConstantObject(SmbIsdelta18oEnum));
160 parameters->AddObject(iomodel->CopyConstantObject(SmbIsmungsmEnum));
161 parameters->AddObject(iomodel->CopyConstantObject(SmbDesfacEnum));
162 parameters->AddObject(iomodel->CopyConstantObject(SmbRlapsEnum));
163 parameters->AddObject(iomodel->CopyConstantObject(SmbRlapslgmEnum));
164 iomodel->Constant(&isdelta18o,SmbIsdelta18oEnum);
165 iomodel->Constant(&ismungsm,SmbIsmungsmEnum);
166
167 if(ismungsm){
168 iomodel->FetchData(&temp,&N,&M,SmbPfacEnum); _assert_(N==2);
169 parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,M));
170 iomodel->DeleteData(temp,SmbPfacEnum);
171
172 iomodel->FetchData(&temp,&N,&M,SmbTdiffEnum); _assert_(N==2);
173 parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,M));
174 iomodel->DeleteData(temp,SmbTdiffEnum);
175
176 iomodel->FetchData(&temp,&N,&M,SmbSealevEnum); _assert_(N==2);
177 parameters->AddObject(new TransientParam(SmbSealevEnum,&temp[0],&temp[M],interp,M));
178 iomodel->DeleteData(temp,SmbSealevEnum);
179 }
180 if(isdelta18o){
181 iomodel->FetchData(&temp,&N,&M,SmbDelta18oEnum); _assert_(N==2);
182 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
183 iomodel->DeleteData(temp,SmbDelta18oEnum);
184
185 iomodel->FetchData(&temp,&N,&M,SmbDelta18oSurfaceEnum); _assert_(N==2);
186 parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,M));
187 iomodel->DeleteData(temp,SmbDelta18oSurfaceEnum);
188 }
189 break;
190 case SMBd18opddEnum:
191 parameters->AddObject(iomodel->CopyConstantObject(SmbIsmungsmEnum));
192 parameters->AddObject(iomodel->CopyConstantObject(SmbIsd18opdEnum));
193 parameters->AddObject(iomodel->CopyConstantObject(SmbDesfacEnum));
194 parameters->AddObject(iomodel->CopyConstantObject(SmbRlapsEnum));
195 parameters->AddObject(iomodel->CopyConstantObject(SmbRlapslgmEnum));
196 iomodel->Constant(&ismungsm,SmbIsmungsmEnum);
197 iomodel->Constant(&isd18opd,SmbIsd18opdEnum);
198 if(isd18opd){
199 iomodel->FetchData(&temp,&N,&M,SmbDelta18oEnum); _assert_(N==2);
200 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
201 iomodel->DeleteData(temp,SmbDelta18oEnum);
202
203 parameters->AddObject(iomodel->CopyConstantObject(SmbDpermilEnum));
204 }
205 break;
206 case SMBgradientsEnum:
207 /*Nothing to add to parameters*/
208 break;
209 case SMBhenningEnum:
210 /*Nothing to add to parameters*/
211 break;
212 case SMBcomponentsEnum:
213 /*Nothing to add to parameters*/
214 break;
215 case SMBmeltcomponentsEnum:
216 /*Nothing to add to parameters*/
217 break;
218 default:
219 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
220 }
221
222 iomodel->FetchData(&requestedoutputs,&numoutputs,SmbRequestedOutputsEnum);
223 parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
224 if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
225 iomodel->DeleteData(&requestedoutputs,numoutputs,SmbRequestedOutputsEnum);
226
227}/*}}}*/
228
229/*Finite Element Analysis*/
230void SmbAnalysis::Core(FemModel* femmodel){/*{{{*/
231
232 int smb_model;
233
234 /*Figure out smb model: */
235 femmodel->parameters->FindParam(&smb_model,SmbEnum);
236
237 /*branch to correct module*/
238 switch(smb_model){
239 case SMBforcingEnum:
240 /*Nothing to be done*/
241 break;
[19554]242 case SMBgembEnum:
243 Gembx(femmodel);
244 break;
[19528]245 case SMBpddEnum:
246 bool isdelta18o,ismungsm;
247 femmodel->parameters->FindParam(&isdelta18o,SmbIsdelta18oEnum);
248 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
249 if(isdelta18o){
250 if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n");
251 Delta18oParameterizationx(femmodel);
252 }
253 if(ismungsm){
254 if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n");
255 MungsmtpParameterizationx(femmodel);
256 }
257 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
258 PositiveDegreeDayx(femmodel);
259 break;
260 case SMBd18opddEnum:
261 bool isd18opd;
262 femmodel->parameters->FindParam(&isd18opd,SmbIsd18opdEnum);
263 if(isd18opd){
264 if(VerboseSolution()) _printf0_(" call Delta18opdParameterization module\n");
265 Delta18opdParameterizationx(femmodel);
266 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
267 PositiveDegreeDayx(femmodel);
268 }
269 break;
270 case SMBgradientsEnum:
271 if(VerboseSolution())_printf0_(" call smb gradients module\n");
272 SmbGradientsx(femmodel);
273 break;
274 case SMBhenningEnum:
275 if(VerboseSolution())_printf0_(" call smb Henning module\n");
276 SmbHenningx(femmodel);
277 break;
278 case SMBcomponentsEnum:
279 if(VerboseSolution())_printf0_(" call smb Components module\n");
280 SmbComponentsx(femmodel);
281 break;
282 case SMBmeltcomponentsEnum:
283 if(VerboseSolution())_printf0_(" call smb Melt Components module\n");
284 SmbMeltComponentsx(femmodel);
285 break;
286 case SMBgcmEnum:
287 /*Nothing to be done*/
288 break;
289 default:
290 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
291 }
292
293}/*}}}*/
294ElementVector* SmbAnalysis::CreateDVector(Element* element){/*{{{*/
295 _error_("not implemented");
296}/*}}}*/
297ElementMatrix* SmbAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
298_error_("Not implemented");
299}/*}}}*/
300ElementMatrix* SmbAnalysis::CreateKMatrix(Element* element){/*{{{*/
301 _error_("not implemented yet");
302}/*}}}*/
303ElementVector* SmbAnalysis::CreatePVector(Element* element){/*{{{*/
304_error_("not implemented yet");
305}/*}}}*/
306void SmbAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
307 _error_("not implemented yet");
308}/*}}}*/
309void SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
310 _error_("Not implemented yet");
311}/*}}}*/
312void SmbAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
313 _error_("not implemented yet");
314}/*}}}*/
315void SmbAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
316 /*Default, do nothing*/
317 return;
318}/*}}}*/
Note: See TracBrowser for help on using the repository browser.