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

Last change on this file since 25317 was 25317, checked in by Mathieu Morlighem, 5 years ago

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