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

Last change on this file since 26477 was 26477, checked in by vverjans, 3 years ago

autoregression SMB module, Random utilities, Cholesky decomposition
File size: 25.7 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// FIX
8#include "./shared/io/Print/Print.h"
9
10/*Model processing*/
11void SmbAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
12 /*No constraints*/
13}/*}}}*/
14void SmbAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
15 /*No loads*/
16}/*}}}*/
17void SmbAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
18 ::CreateNodes(nodes,iomodel,SmbAnalysisEnum,P1Enum);
19}/*}}}*/
20int SmbAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
21 return 1;
22}/*}}}*/
23void SmbAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
24
25 int smb_model;
26 bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled,isfirnwarming;
27
28 /*Update elements: */
29 int counter=0;
30 for(int i=0;i<iomodel->numberofelements;i++){
31 if(iomodel->my_elements[i]){
32 Element* element=(Element*)elements->GetObjectByOffset(counter);
33 element->Update(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum);
34 counter++;
35 }
36 }
37
38 /*Figure out smb model: */
39 iomodel->FindConstant(&smb_model,"md.smb.model");
40 switch(smb_model){
41 case SMBforcingEnum:
42 iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
43 break;
44 case SMBgembEnum:
45 iomodel->FetchDataToInput(inputs,elements,"md.smb.Ta",SmbTaEnum);
46 iomodel->FetchDataToInput(inputs,elements,"md.smb.V",SmbVEnum);
47 iomodel->FetchDataToInput(inputs,elements,"md.smb.dswrf",SmbDswrfEnum);
48 iomodel->FetchDataToInput(inputs,elements,"md.smb.dswdiffrf",SmbDswdiffrfEnum);
49 iomodel->FetchDataToInput(inputs,elements,"md.smb.dlwrf",SmbDlwrfEnum);
50 iomodel->FetchDataToInput(inputs,elements,"md.smb.P",SmbPEnum);
51 iomodel->FetchDataToInput(inputs,elements,"md.smb.eAir",SmbEAirEnum);
52 iomodel->FetchDataToInput(inputs,elements,"md.smb.pAir",SmbPAirEnum);
53 iomodel->FetchDataToInput(inputs,elements,"md.smb.zTop",SmbZTopEnum);
54 iomodel->FetchDataToInput(inputs,elements,"md.smb.dzTop",SmbDzTopEnum);
55 iomodel->FetchDataToInput(inputs,elements,"md.smb.dzMin",SmbDzMinEnum);
56 iomodel->FetchDataToInput(inputs,elements,"md.smb.zY",SmbZYEnum);
57 iomodel->FetchDataToInput(inputs,elements,"md.smb.zMax",SmbZMaxEnum);
58 iomodel->FetchDataToInput(inputs,elements,"md.smb.zMin",SmbZMinEnum);
59 iomodel->FetchDataToInput(inputs,elements,"md.smb.Tmean",SmbTmeanEnum);
60 iomodel->FetchDataToInput(inputs,elements,"md.smb.Vmean",SmbVmeanEnum);
61 iomodel->FetchDataToInput(inputs,elements,"md.smb.C",SmbCEnum);
62 iomodel->FetchDataToInput(inputs,elements,"md.smb.Tz",SmbTzEnum);
63 iomodel->FetchDataToInput(inputs,elements,"md.smb.Vz",SmbVzEnum);
64 InputUpdateFromConstantx(inputs,elements,false,SmbIsInitializedEnum);
65 iomodel->FetchDataToInput(inputs,elements,"md.smb.Dzini",SmbDziniEnum);
66 iomodel->FetchDataToInput(inputs,elements,"md.smb.Dini",SmbDiniEnum);
67 iomodel->FetchDataToInput(inputs,elements,"md.smb.Reini",SmbReiniEnum);
68 iomodel->FetchDataToInput(inputs,elements,"md.smb.Gdnini",SmbGdniniEnum);
69 iomodel->FetchDataToInput(inputs,elements,"md.smb.Gspini",SmbGspiniEnum);
70 iomodel->FetchDataToInput(inputs,elements,"md.smb.ECini",SmbECiniEnum);
71 iomodel->FetchDataToInput(inputs,elements,"md.smb.Wini",SmbWiniEnum);
72 iomodel->FetchDataToInput(inputs,elements,"md.smb.Aini",SmbAiniEnum);
73 iomodel->FetchDataToInput(inputs,elements,"md.smb.Adiffini",SmbAdiffiniEnum);
74 iomodel->FetchDataToInput(inputs,elements,"md.smb.Tini",SmbTiniEnum);
75 iomodel->FetchDataToInput(inputs,elements,"md.smb.Sizeini",SmbSizeiniEnum);
76 iomodel->FetchDataToInput(inputs,elements,"md.smb.aValue",SmbAValueEnum);
77 iomodel->FetchDataToInput(inputs,elements,"md.smb.teValue",SmbTeValueEnum);
78 iomodel->FetchDataToInput(inputs,elements,"md.smb.szaValue",SmbSzaValueEnum);
79 iomodel->FetchDataToInput(inputs,elements,"md.smb.cotValue",SmbCotValueEnum);
80 iomodel->FetchDataToInput(inputs,elements,"md.smb.ccsnowValue",SmbCcsnowValueEnum);
81 iomodel->FetchDataToInput(inputs,elements,"md.smb.cciceValue",SmbCciceValueEnum);
82 break;
83 case SMBpddEnum:
84 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
85 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
86 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
87 iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
88 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
89 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
90 if(isdelta18o || ismungsm){
91 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_lgm",SmbTemperaturesLgmEnum);
92 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
93 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
94 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum);
95 }else{
96 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitation",SmbPrecipitationEnum);
97 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
98 }
99 if(issetpddfac){
100 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
101 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
102 }
103 break;
104 case SMBpddSicopolisEnum:
105 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
106 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
107 iomodel->FindConstant(&isfirnwarming,"md.smb.isfirnwarming");
108 iomodel->FetchDataToInput(inputs,elements,"md.smb.smb_corr",SmbSmbCorrEnum);
109 iomodel->FetchDataToInput(inputs,elements,"md.smb.precipitation_anomaly",SmbPrecipitationsAnomalyEnum);
110 iomodel->FetchDataToInput(inputs,elements,"md.smb.temperature_anomaly",SmbTemperaturesAnomalyEnum);
111 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
112 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitation",SmbPrecipitationEnum);
113 break;
114 case SMBd18opddEnum:
115 iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
116 iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
117 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
118 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
119 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
120 iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
121 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
122 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
123 if(isd18opd){
124 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
125 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
126 if(!istemperaturescaled){
127 /*Fetch array*/
128 IssmDouble* doublearray = NULL;
129 int M,N;
130 iomodel->FetchData(&doublearray,&M,&N,"md.smb.temperatures_reconstructed");
131 if(M!=iomodel->numberofvertices+1) _error_("md.smb.temperatures_reconstructed should have nbv+1 rows");
132 if(N%12!=0) _error_("md.smb.temperatures_reconstructed should have a multiple of 12 columns (since it is monthly)");
133
134 /*Build indices*/
135 int* ids = xNew<int>(N); for(int i=0;i<N;i++) ids[i] = i;
136
137 for(Object* & object : elements->objects){
138 Element* element=xDynamicCast<Element*>(object);
139 element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs,iomodel,SmbTemperaturesReconstructedEnum);
140 }
141 xDelete<int>(ids);
142 iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
143 }
144 if(!isprecipscaled){
145 /*Fetch array*/
146 IssmDouble* doublearray = NULL;
147 int M,N;
148 iomodel->FetchData(&doublearray,&M,&N,"md.smb.precipitations_reconstructed");
149 if(M!=iomodel->numberofvertices+1) _error_("md.smb.precipitations_reconstructed should have nbv+1 rows");
150 if(N%12!=0) _error_("md.smb.precipitations_reconstructed should have a multiple of 12 columns (since it is monthly)");
151
152 /*Build indices*/
153 int* ids = xNew<int>(N); for(int i=0;i<N;i++) ids[i] = i;
154
155 for(Object* & object : elements->objects){
156 Element* element=xDynamicCast<Element*>(object);
157 element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs,iomodel,SmbPrecipitationsReconstructedEnum);
158 }
159 xDelete<int>(ids);
160 iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
161 }
162 }
163 if(issetpddfac){
164 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
165 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
166 }
167 break;
168 case SMBgradientsEnum:
169 iomodel->FetchDataToInput(inputs,elements,"md.smb.href",SmbHrefEnum);
170 iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum);
171 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_pos",SmbBPosEnum);
172 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_neg",SmbBNegEnum);
173 break;
174 case SMBgradientselaEnum:
175 iomodel->FetchDataToInput(inputs,elements,"md.smb.ela",SmbElaEnum);
176 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_pos",SmbBPosEnum);
177 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_neg",SmbBNegEnum);
178 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_max",SmbBMaxEnum);
179 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_min",SmbBMinEnum);
180 break;
181 case SMBautoregressionEnum:
182 iomodel->FetchDataToInput(inputs,elements,"md.smb.basin_id",SmbBasinsIdEnum);
183 break;
184 case SMBhenningEnum:
185 iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum,0.);
186 break;
187 case SMBcomponentsEnum:
188 iomodel->FetchDataToInput(inputs,elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
189 iomodel->FetchDataToInput(inputs,elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
190 iomodel->FetchDataToInput(inputs,elements,"md.smb.runoff",SmbRunoffEnum,0.);
191 break;
192 case SMBmeltcomponentsEnum:
193 iomodel->FetchDataToInput(inputs,elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
194 iomodel->FetchDataToInput(inputs,elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
195 iomodel->FetchDataToInput(inputs,elements,"md.smb.melt",SmbMeltEnum,0.);
196 iomodel->FetchDataToInput(inputs,elements,"md.smb.refreeze",SmbRefreezeEnum,0.);
197 break;
198 case SMBgradientscomponentsEnum:
199 /* Nothing to add to input */
200 break;
201 case SMBsemicEnum:
202 iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
203 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0gcm",SmbS0gcmEnum);
204 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailysnowfall",SmbDailysnowfallEnum);
205 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyrainfall",SmbDailyrainfallEnum);
206 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailydsradiation",SmbDailydsradiationEnum);
207 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailydlradiation",SmbDailydlradiationEnum);
208 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailywindspeed",SmbDailywindspeedEnum);
209 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailypressure",SmbDailypressureEnum);
210 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyairdensity",SmbDailyairdensityEnum);
211 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyairhumidity",SmbDailyairhumidityEnum);
212 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailytemperature",SmbDailytemperatureEnum);
213 break;
214 default:
215 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
216 }
217
218}/*}}}*/
219void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
220
221 int numoutputs;
222 char** requestedoutputs = NULL;
223 bool isdelta18o,ismungsm,isd18opd,issetpddfac,interp,cycle,isfirnwarming;
224 int smb_model, smbslices, averaging;
225 IssmDouble *temp = NULL;
226 int N,M;
227
228 parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
229
230 iomodel->FindConstant(&smb_model,"md.smb.model");
231 iomodel->FindConstant(&interp,"md.timestepping.interp_forcing");
232 iomodel->FindConstant(&cycle,"md.timestepping.cycle_forcing");
233
234 iomodel->FindConstant(&smbslices,"md.smb.steps_per_step");
235 parameters->AddObject(new IntParam(SmbStepsPerStepEnum,smbslices));
236
237 parameters->AddObject(iomodel->CopyConstantObject("md.smb.averaging",SmbAveragingEnum));
238
239 switch(smb_model){
240 case SMBforcingEnum:
241 break;
242 case SMBgembEnum:
243 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIdx",SmbAIdxEnum));
244 parameters->AddObject(iomodel->CopyConstantObject("md.smb.swIdx",SmbSwIdxEnum));
245 parameters->AddObject(iomodel->CopyConstantObject("md.smb.denIdx",SmbDenIdxEnum));
246 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dsnowIdx",SmbDsnowIdxEnum));
247 parameters->AddObject(iomodel->CopyConstantObject("md.smb.cldFrac",SmbCldFracEnum));
248 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0wet",SmbT0wetEnum));
249 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0dry",SmbT0dryEnum));
250 parameters->AddObject(iomodel->CopyConstantObject("md.smb.K",SmbKEnum));
251 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aSnow",SmbASnowEnum));
252 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));
253 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dt",SmbDtEnum));
254 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isgraingrowth",SmbIsgraingrowthEnum));
255 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isalbedo",SmbIsalbedoEnum));
256 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isshortwave",SmbIsshortwaveEnum));
257 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isthermal",SmbIsthermalEnum));
258 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isaccumulation",SmbIsaccumulationEnum));
259 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismelt",SmbIsmeltEnum));
260 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdensification",SmbIsdensificationEnum));
261 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum));
262 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isconstrainsurfaceT",SmbIsconstrainsurfaceTEnum));
263 parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
264 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
265 parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));
266 break;
267 case SMBpddEnum:
268 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdelta18o",SmbIsdelta18oEnum));
269 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
270 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
271 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
272 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
273
274 if(ismungsm){
275 iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2);
276 parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,cycle,M));
277 iomodel->DeleteData(temp,"md.smb.Pfac");
278
279 iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
280 parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,cycle,M));
281 iomodel->DeleteData(temp,"md.smb.Tdiff");
282
283 iomodel->FetchData(&temp,&N,&M,"md.smb.sealev"); _assert_(N==2);
284 parameters->AddObject(new TransientParam(SmbSealevEnum,&temp[0],&temp[M],interp,cycle,M));
285 iomodel->DeleteData(temp,"md.smb.sealev");
286 }
287 if(isdelta18o){
288 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
289 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,cycle,M));
290 iomodel->DeleteData(temp,"md.smb.delta18o");
291
292 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o_surface"); _assert_(N==2);
293 parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,cycle,M));
294 iomodel->DeleteData(temp,"md.smb.delta18o_surface");
295 }
296
297 break;
298 case SMBpddSicopolisEnum:
299 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isfirnwarming",SmbIsfirnwarmingEnum));
300 break;
301 case SMBd18opddEnum:
302 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
303 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isd18opd",SmbIsd18opdEnum));
304 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
305 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
306 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
307 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
308 if(isd18opd){
309 parameters->AddObject(iomodel->CopyConstantObject("md.smb.f",SmbFEnum));
310 parameters->AddObject(iomodel->CopyConstantObject("md.smb.istemperaturescaled",SmbIstemperaturescaledEnum));
311 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipscaled",SmbIsprecipscaledEnum));
312 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
313 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,cycle,M));
314 iomodel->DeleteData(temp,"md.smb.delta18o");
315
316 IssmDouble yts;
317 bool istemperaturescaled,isprecipscaled;
318 iomodel->FindConstant(&yts,"md.constants.yts");
319 iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
320 iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
321 if(!istemperaturescaled){
322 /*Fetch array*/
323 IssmDouble* doublearray = NULL;
324 int M,N;
325 iomodel->FetchData(&doublearray,&M,&N,"md.smb.temperatures_reconstructed");
326 if(M!=iomodel->numberofvertices+1) _error_("md.smb.temperatures_reconstructed should have nbv+1 rows");
327 if(N%12!=0) _error_("md.smb.temperatures_reconstructed should have a multiple of 12 columns (since it is monthly)");
328 int numyears = N/12; _assert_(numyears*12==N);
329
330 /*Check times*/
331 #ifdef _ISSM_DEBUG_
332 for(int i=0;i<numyears;i++){
333 for(int j=1;j<12;j++){
334 //_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
335 _assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
336 }
337 }
338 #endif
339
340 /*Build time*/
341 IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
342 parameters->AddObject(new DoubleVecParam(SmbTemperaturesReconstructedYearsEnum,times,numyears));
343 xDelete<IssmDouble>(times);
344 iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
345 }
346 if(!isprecipscaled){
347 /*Fetch array*/
348 IssmDouble* doublearray = NULL;
349 int M,N;
350 iomodel->FetchData(&doublearray,&M,&N,"md.smb.precipitations_reconstructed");
351 if(M!=iomodel->numberofvertices+1) _error_("md.smb.precipitations_reconstructed should have nbv+1 rows");
352 if(N%12!=0) _error_("md.smb.precipitations_reconstructed should have a multiple of 12 columns (since it is monthly)");
353 int numyears = N/12; _assert_(numyears*12==N);
354
355 /*Check times*/
356 #ifdef _ISSM_DEBUG_
357 for(int i=0;i<numyears;i++){
358 for(int j=1;j<12;j++){
359 //_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
360 _assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
361 }
362 }
363 #endif
364
365 /*Build time*/
366 IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
367 parameters->AddObject(new DoubleVecParam(SmbPrecipitationsReconstructedYearsEnum,times,numyears));
368 xDelete<IssmDouble>(times);
369 iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
370 }
371 }
372
373 break;
374 case SMBgradientsEnum:
375 /*Nothing to add to parameters*/
376 break;
377 case SMBgradientselaEnum:
378 /*Nothing to add to parameters*/
379 break;
380 case SMBautoregressionEnum:
381 /*Nothing to add to parameters*/
382 break;
383 case SMBhenningEnum:
384 /*Nothing to add to parameters*/
385 break;
386 case SMBcomponentsEnum:
387 break;
388 case SMBmeltcomponentsEnum:
389 break;
390 case SMBgradientscomponentsEnum:
391 parameters->AddObject(iomodel->CopyConstantObject("md.smb.accualti",SmbAccualtiEnum));
392 parameters->AddObject(iomodel->CopyConstantObject("md.smb.runoffalti",SmbRunoffaltiEnum));
393
394 iomodel->FetchData(&temp,&N,&M,"md.smb.accugrad"); _assert_(N==2);
395 parameters->AddObject(new TransientParam(SmbAccugradEnum,&temp[0],&temp[M],interp,cycle,M));
396 iomodel->DeleteData(temp,"md.smb.accugrad");
397 iomodel->FetchData(&temp,&N,&M,"md.smb.runoffgrad"); _assert_(N==2);
398 parameters->AddObject(new TransientParam(SmbRunoffgradEnum,&temp[0],&temp[M],interp,cycle,M));
399 iomodel->DeleteData(temp,"md.smb.runoffgrad");
400
401 iomodel->FetchData(&temp,&N,&M,"md.smb.accuref"); _assert_(N==2);
402 parameters->AddObject(new TransientParam(SmbAccurefEnum,&temp[0],&temp[M],interp,cycle,M));
403 iomodel->DeleteData(temp,"md.smb.accuref");
404 iomodel->FetchData(&temp,&N,&M,"md.smb.runoffref"); _assert_(N==2);
405 parameters->AddObject(new TransientParam(SmbRunoffrefEnum,&temp[0],&temp[M],interp,cycle,M));
406 iomodel->DeleteData(temp,"md.smb.runoffref");
407 break;
408 case SMBsemicEnum:
409 /*Nothing to add to parameters*/
410 break;
411 default:
412 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
413 }
414
415 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.smb.requested_outputs");
416 parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
417 if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
418 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.smb.requested_outputs");
419
420}/*}}}*/
421
422/*Finite Element Analysis*/
423void SmbAnalysis::Core(FemModel* femmodel){/*{{{*/
424
425 int smb_model;
426
427 /*Figure out smb model: */
428 femmodel->parameters->FindParam(&smb_model,SmbEnum);
429
430 /*branch to correct module*/
431 switch(smb_model){
432 case SMBforcingEnum:
433 SmbForcingx(femmodel);
434 break;
435 case SMBgembEnum:
436 Gembx(femmodel);
437 break;
438 case SMBpddEnum:
439 bool isdelta18o,ismungsm;
440 femmodel->parameters->FindParam(&isdelta18o,SmbIsdelta18oEnum);
441 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
442 if(isdelta18o){
443 if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n");
444 Delta18oParameterizationx(femmodel);
445 }
446 if(ismungsm){
447 if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n");
448 MungsmtpParameterizationx(femmodel);
449 }
450 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
451 PositiveDegreeDayx(femmodel);
452 break;
453 case SMBpddSicopolisEnum:
454 if(VerboseSolution()) _printf0_(" call SICOPOLIS positive degree day module\n");
455 PositiveDegreeDaySicopolisx(femmodel);
456 break;
457 case SMBd18opddEnum:
458 bool isd18opd;
459 femmodel->parameters->FindParam(&isd18opd,SmbIsd18opdEnum);
460 if(isd18opd){
461 if(VerboseSolution()) _printf0_(" call Delta18opdParameterization module\n");
462 Delta18opdParameterizationx(femmodel);
463 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
464 PositiveDegreeDayx(femmodel);
465 }
466 break;
467 case SMBgradientsEnum:
468 if(VerboseSolution())_printf0_(" call smb gradients module\n");
469 SmbGradientsx(femmodel);
470 break;
471 case SMBgradientselaEnum:
472 if(VerboseSolution())_printf0_(" call smb gradients ela module\n");
473 SmbGradientsElax(femmodel);
474 break;
475 case SMBautoregressionEnum:
476 if(VerboseSolution())_printf0_(" call smb autoregression module\n");
477 Smbautoregressionx(femmodel);
478 break;
479 case SMBhenningEnum:
480 if(VerboseSolution())_printf0_(" call smb Henning module\n");
481 SmbHenningx(femmodel);
482 break;
483 case SMBcomponentsEnum:
484 if(VerboseSolution())_printf0_(" call smb Components module\n");
485 SmbComponentsx(femmodel);
486 break;
487 case SMBmeltcomponentsEnum:
488 if(VerboseSolution())_printf0_(" call smb Melt Components module\n");
489 SmbMeltComponentsx(femmodel);
490 break;
491 case SMBgcmEnum:
492 /*Nothing to be done*/
493 break;
494 case SMBgradientscomponentsEnum:
495 if(VerboseSolution())_printf0_(" call smb gradients components module\n");
496 SmbGradientsComponentsx(femmodel);
497 break;
498 case SMBsemicEnum:
499 #ifdef _HAVE_SEMIC_
500 if(VerboseSolution())_printf0_(" call smb SEMIC module\n");
501 SmbSemicx(femmodel);
502 #else
503 _error_("SEMIC not installed");
504 #endif //_HAVE_SEMIC_
505 break;
506 default:
507 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
508 }
509
510}/*}}}*/
511void SmbAnalysis::PreCore(FemModel* femmodel){/*{{{*/
512 _error_("not implemented");
513}/*}}}*/
514ElementVector* SmbAnalysis::CreateDVector(Element* element){/*{{{*/
515 _error_("not implemented");
516}/*}}}*/
517ElementMatrix* SmbAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
518_error_("Not implemented");
519}/*}}}*/
520ElementMatrix* SmbAnalysis::CreateKMatrix(Element* element){/*{{{*/
521 _error_("not implemented yet");
522}/*}}}*/
523ElementVector* SmbAnalysis::CreatePVector(Element* element){/*{{{*/
524_error_("not implemented yet");
525}/*}}}*/
526void SmbAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
527 _error_("not implemented yet");
528}/*}}}*/
529void SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
530 _error_("Not implemented yet");
531}/*}}}*/
532void SmbAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
533 _error_("not implemented yet");
534}/*}}}*/
535void SmbAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
536 /*Default, do nothing*/
537 return;
538}/*}}}*/
Note: See TracBrowser for help on using the repository browser.