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

Last change on this file since 19565 was 19565, checked in by Eric.Larour, 10 years ago

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