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

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

CHG: huge commit on solid earth capability rewrite. Complete cleanup of the sea level core.
New mass transport capabilities for ocean and tws. No more giacore. GiaIvins folded into
the sea level core. Debugging of Materials.

File size: 25.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// 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 SMBhenningEnum:
182 iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum,0.);
183 break;
184 case SMBcomponentsEnum:
185 iomodel->FetchDataToInput(inputs,elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
186 iomodel->FetchDataToInput(inputs,elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
187 iomodel->FetchDataToInput(inputs,elements,"md.smb.runoff",SmbRunoffEnum,0.);
188 break;
189 case SMBmeltcomponentsEnum:
190 iomodel->FetchDataToInput(inputs,elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
191 iomodel->FetchDataToInput(inputs,elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
192 iomodel->FetchDataToInput(inputs,elements,"md.smb.melt",SmbMeltEnum,0.);
193 iomodel->FetchDataToInput(inputs,elements,"md.smb.refreeze",SmbRefreezeEnum,0.);
194 break;
195 case SMBgradientscomponentsEnum:
196 /* Nothing to add to input */
197 break;
198 case SMBsemicEnum:
199 iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
200 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0gcm",SmbS0gcmEnum);
201 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailysnowfall",SmbDailysnowfallEnum);
202 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyrainfall",SmbDailyrainfallEnum);
203 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailydsradiation",SmbDailydsradiationEnum);
204 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailydlradiation",SmbDailydlradiationEnum);
205 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailywindspeed",SmbDailywindspeedEnum);
206 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailypressure",SmbDailypressureEnum);
207 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyairdensity",SmbDailyairdensityEnum);
208 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyairhumidity",SmbDailyairhumidityEnum);
209 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailytemperature",SmbDailytemperatureEnum);
210 break;
211 default:
212 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
213 }
214
215}/*}}}*/
216void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
217
218 int numoutputs;
219 char** requestedoutputs = NULL;
220 bool isdelta18o,ismungsm,isd18opd,issetpddfac,interp,isfirnwarming;
221 int smb_model, smbslices, averaging;
222 IssmDouble *temp = NULL;
223 int N,M;
224
225 parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
226
227 iomodel->FindConstant(&smb_model,"md.smb.model");
228 iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
229
230 iomodel->FindConstant(&smbslices,"md.smb.steps_per_step");
231 parameters->AddObject(new IntParam(SmbStepsPerStepEnum,smbslices));
232
233
234 parameters->AddObject(iomodel->CopyConstantObject("md.smb.averaging",SmbAveragingEnum));
235
236 switch(smb_model){
237 case SMBforcingEnum:
238 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
239 break;
240 case SMBgembEnum:
241 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIdx",SmbAIdxEnum));
242 parameters->AddObject(iomodel->CopyConstantObject("md.smb.swIdx",SmbSwIdxEnum));
243 parameters->AddObject(iomodel->CopyConstantObject("md.smb.denIdx",SmbDenIdxEnum));
244 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dsnowIdx",SmbDsnowIdxEnum));
245 parameters->AddObject(iomodel->CopyConstantObject("md.smb.cldFrac",SmbCldFracEnum));
246 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0wet",SmbT0wetEnum));
247 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0dry",SmbT0dryEnum));
248 parameters->AddObject(iomodel->CopyConstantObject("md.smb.K",SmbKEnum));
249 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aSnow",SmbASnowEnum));
250 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));
251 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dt",SmbDtEnum));
252 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isgraingrowth",SmbIsgraingrowthEnum));
253 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isalbedo",SmbIsalbedoEnum));
254 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isshortwave",SmbIsshortwaveEnum));
255 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isthermal",SmbIsthermalEnum));
256 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isaccumulation",SmbIsaccumulationEnum));
257 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismelt",SmbIsmeltEnum));
258 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdensification",SmbIsdensificationEnum));
259 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum));
260 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
261 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isconstrainsurfaceT",SmbIsconstrainsurfaceTEnum));
262 parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
263 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
264 parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));
265 break;
266 case SMBpddEnum:
267 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdelta18o",SmbIsdelta18oEnum));
268 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
269 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
270 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
271 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
272
273 if(ismungsm){
274 iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2);
275 parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,M));
276 iomodel->DeleteData(temp,"md.smb.Pfac");
277
278 iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
279 parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,M));
280 iomodel->DeleteData(temp,"md.smb.Tdiff");
281
282 iomodel->FetchData(&temp,&N,&M,"md.smb.sealev"); _assert_(N==2);
283 parameters->AddObject(new TransientParam(SmbSealevEnum,&temp[0],&temp[M],interp,M));
284 iomodel->DeleteData(temp,"md.smb.sealev");
285 }
286 if(isdelta18o){
287 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
288 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
289 iomodel->DeleteData(temp,"md.smb.delta18o");
290
291 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o_surface"); _assert_(N==2);
292 parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,M));
293 iomodel->DeleteData(temp,"md.smb.delta18o_surface");
294 }
295
296 break;
297 case SMBpddSicopolisEnum:
298 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isfirnwarming",SmbIsfirnwarmingEnum));
299 break;
300 case SMBd18opddEnum:
301 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
302 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isd18opd",SmbIsd18opdEnum));
303 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
304 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
305 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
306 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
307 if(isd18opd){
308 parameters->AddObject(iomodel->CopyConstantObject("md.smb.f",SmbFEnum));
309 parameters->AddObject(iomodel->CopyConstantObject("md.smb.istemperaturescaled",SmbIstemperaturescaledEnum));
310 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipscaled",SmbIsprecipscaledEnum));
311 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
312 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
313 iomodel->DeleteData(temp,"md.smb.delta18o");
314
315 IssmDouble yts;
316 bool istemperaturescaled,isprecipscaled;
317 iomodel->FindConstant(&yts,"md.constants.yts");
318 iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
319 iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
320 if(!istemperaturescaled){
321 /*Fetch array*/
322 IssmDouble* doublearray = NULL;
323 int M,N;
324 iomodel->FetchData(&doublearray,&M,&N,"md.smb.temperatures_reconstructed");
325 if(M!=iomodel->numberofvertices+1) _error_("md.smb.temperatures_reconstructed should have nbv+1 rows");
326 if(N%12!=0) _error_("md.smb.temperatures_reconstructed should have a multiple of 12 columns (since it is monthly)");
327 int numyears = N/12; _assert_(numyears*12==N);
328
329 /*Check times*/
330 #ifdef _ISSM_DEBUG_
331 for(int i=0;i<numyears;i++){
332 for(int j=1;j<12;j++){
333 //_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
334 _assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
335 }
336 }
337 #endif
338
339 /*Build time*/
340 IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
341 parameters->AddObject(new DoubleVecParam(SmbTemperaturesReconstructedYearsEnum,times,numyears));
342 xDelete<IssmDouble>(times);
343 iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
344 }
345 if(!isprecipscaled){
346 /*Fetch array*/
347 IssmDouble* doublearray = NULL;
348 int M,N;
349 iomodel->FetchData(&doublearray,&M,&N,"md.smb.precipitations_reconstructed");
350 if(M!=iomodel->numberofvertices+1) _error_("md.smb.precipitations_reconstructed should have nbv+1 rows");
351 if(N%12!=0) _error_("md.smb.precipitations_reconstructed should have a multiple of 12 columns (since it is monthly)");
352 int numyears = N/12; _assert_(numyears*12==N);
353
354 /*Check times*/
355 #ifdef _ISSM_DEBUG_
356 for(int i=0;i<numyears;i++){
357 for(int j=1;j<12;j++){
358 //_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
359 _assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
360 }
361 }
362 #endif
363
364 /*Build time*/
365 IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
366 parameters->AddObject(new DoubleVecParam(SmbPrecipitationsReconstructedYearsEnum,times,numyears));
367 xDelete<IssmDouble>(times);
368 iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
369 }
370 }
371
372 break;
373 case SMBgradientsEnum:
374 /*Nothing to add to parameters*/
375 break;
376 case SMBgradientselaEnum:
377 /*Nothing to add to parameters*/
378 break;
379 case SMBhenningEnum:
380 /*Nothing to add to parameters*/
381 break;
382 case SMBcomponentsEnum:
383 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
384 break;
385 case SMBmeltcomponentsEnum:
386 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
387 break;
388 case SMBgradientscomponentsEnum:
389 parameters->AddObject(iomodel->CopyConstantObject("md.smb.accualti",SmbAccualtiEnum));
390 parameters->AddObject(iomodel->CopyConstantObject("md.smb.runoffalti",SmbRunoffaltiEnum));
391
392 iomodel->FetchData(&temp,&N,&M,"md.smb.accugrad"); _assert_(N==2);
393 parameters->AddObject(new TransientParam(SmbAccugradEnum,&temp[0],&temp[M],interp,M));
394 iomodel->DeleteData(temp,"md.smb.accugrad");
395 iomodel->FetchData(&temp,&N,&M,"md.smb.runoffgrad"); _assert_(N==2);
396 parameters->AddObject(new TransientParam(SmbRunoffgradEnum,&temp[0],&temp[M],interp,M));
397 iomodel->DeleteData(temp,"md.smb.runoffgrad");
398
399 iomodel->FetchData(&temp,&N,&M,"md.smb.accuref"); _assert_(N==2);
400 parameters->AddObject(new TransientParam(SmbAccurefEnum,&temp[0],&temp[M],interp,M));
401 iomodel->DeleteData(temp,"md.smb.accuref");
402 iomodel->FetchData(&temp,&N,&M,"md.smb.runoffref"); _assert_(N==2);
403 parameters->AddObject(new TransientParam(SmbRunoffrefEnum,&temp[0],&temp[M],interp,M));
404 iomodel->DeleteData(temp,"md.smb.runoffref");
405 break;
406 case SMBsemicEnum:
407 /*Nothing to add to parameters*/
408 break;
409 default:
410 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
411 }
412
413 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.smb.requested_outputs");
414 parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
415 if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
416 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.smb.requested_outputs");
417
418}/*}}}*/
419
420/*Finite Element Analysis*/
421void SmbAnalysis::Core(FemModel* femmodel){/*{{{*/
422
423 int smb_model;
424
425 /*Figure out smb model: */
426 femmodel->parameters->FindParam(&smb_model,SmbEnum);
427
428 /*branch to correct module*/
429 switch(smb_model){
430 case SMBforcingEnum:
431 SmbForcingx(femmodel);
432 break;
433 case SMBgembEnum:
434 Gembx(femmodel);
435 break;
436 case SMBpddEnum:
437 bool isdelta18o,ismungsm;
438 femmodel->parameters->FindParam(&isdelta18o,SmbIsdelta18oEnum);
439 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
440 if(isdelta18o){
441 if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n");
442 Delta18oParameterizationx(femmodel);
443 }
444 if(ismungsm){
445 if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n");
446 MungsmtpParameterizationx(femmodel);
447 }
448 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
449 PositiveDegreeDayx(femmodel);
450 break;
451 case SMBpddSicopolisEnum:
452 if(VerboseSolution()) _printf0_(" call SICOPOLIS positive degree day module\n");
453 PositiveDegreeDaySicopolisx(femmodel);
454 break;
455 case SMBd18opddEnum:
456 bool isd18opd;
457 femmodel->parameters->FindParam(&isd18opd,SmbIsd18opdEnum);
458 if(isd18opd){
459 if(VerboseSolution()) _printf0_(" call Delta18opdParameterization module\n");
460 Delta18opdParameterizationx(femmodel);
461 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
462 PositiveDegreeDayx(femmodel);
463 }
464 break;
465 case SMBgradientsEnum:
466 if(VerboseSolution())_printf0_(" call smb gradients module\n");
467 SmbGradientsx(femmodel);
468 break;
469 case SMBgradientselaEnum:
470 if(VerboseSolution())_printf0_(" call smb gradients ela module\n");
471 SmbGradientsElax(femmodel);
472 break;
473 case SMBhenningEnum:
474 if(VerboseSolution())_printf0_(" call smb Henning module\n");
475 SmbHenningx(femmodel);
476 break;
477 case SMBcomponentsEnum:
478 if(VerboseSolution())_printf0_(" call smb Components module\n");
479 SmbComponentsx(femmodel);
480 break;
481 case SMBmeltcomponentsEnum:
482 if(VerboseSolution())_printf0_(" call smb Melt Components module\n");
483 SmbMeltComponentsx(femmodel);
484 break;
485 case SMBgcmEnum:
486 /*Nothing to be done*/
487 break;
488 case SMBgradientscomponentsEnum:
489 if(VerboseSolution())_printf0_(" call smb gradients components module\n");
490 SmbGradientsComponentsx(femmodel);
491 break;
492 case SMBsemicEnum:
493 #ifdef _HAVE_SEMIC_
494 if(VerboseSolution())_printf0_(" call smb SEMIC module\n");
495 SmbSemicx(femmodel);
496 #else
497 _error_("SEMIC not installed");
498 #endif //_HAVE_SEMIC_
499 break;
500 default:
501 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
502 }
503
504}/*}}}*/
505void SmbAnalysis::PreCore(FemModel* femmodel){/*{{{*/
506 _error_("not implemented");
507}/*}}}*/
508ElementVector* SmbAnalysis::CreateDVector(Element* element){/*{{{*/
509 _error_("not implemented");
510}/*}}}*/
511ElementMatrix* SmbAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
512_error_("Not implemented");
513}/*}}}*/
514ElementMatrix* SmbAnalysis::CreateKMatrix(Element* element){/*{{{*/
515 _error_("not implemented yet");
516}/*}}}*/
517ElementVector* SmbAnalysis::CreatePVector(Element* element){/*{{{*/
518_error_("not implemented yet");
519}/*}}}*/
520void SmbAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
521 _error_("not implemented yet");
522}/*}}}*/
523void SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
524 _error_("Not implemented yet");
525}/*}}}*/
526void SmbAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
527 _error_("not implemented yet");
528}/*}}}*/
529void SmbAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
530 /*Default, do nothing*/
531 return;
532}/*}}}*/
Note: See TracBrowser for help on using the repository browser.