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

Last change on this file since 25997 was 25997, checked in by schlegel, 4 years ago

CHG: implement inputs for gardner albedo scheme
File size: 25.5 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}/*}}}*/
[25379]23void SmbAnalysis::UpdateElements(Elements* elements,Inputs* inputs,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);
[25379]33 element->Update(inputs,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:
[25379]42 iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
[19528]43 break;
[19554]44 case SMBgembEnum:
[25379]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);
[25997]48 iomodel->FetchDataToInput(inputs,elements,"md.smb.dswdiffrf",SmbDswdiffrfEnum);
[25379]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);
[25997]73 iomodel->FetchDataToInput(inputs,elements,"md.smb.Adiffini",SmbAdiffiniEnum);
[25379]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);
[25997]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);
[19554]82 break;
[19528]83 case SMBpddEnum:
[20690]84 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
85 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
[25374]86 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
[25379]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);
[19528]90 if(isdelta18o || ismungsm){
[25379]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);
[22852]95 }else{
[25379]96 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitation",SmbPrecipitationEnum);
97 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
[19528]98 }
[25374]99 if(issetpddfac){
[25379]100 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
101 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
[25374]102 }
[19528]103 break;
[23317]104 case SMBpddSicopolisEnum:
[25379]105 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
106 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
[23317]107 iomodel->FindConstant(&isfirnwarming,"md.smb.isfirnwarming");
[25379]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);
[23317]113 break;
[19528]114 case SMBd18opddEnum:
[22852]115 iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
116 iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
[20690]117 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
118 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
[22448]119 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
[25379]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);
[20690]123 if(isd18opd){
[25379]124 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
125 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
[22852]126 if(!istemperaturescaled){
[24626]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
[25539]137 for(Object* & object : elements->objects){
138 Element* element=xDynamicCast<Element*>(object);
[25379]139 element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs,iomodel,SmbTemperaturesReconstructedEnum);
[24626]140 }
141 xDelete<int>(ids);
142 iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
[22852]143 }
144 if(!isprecipscaled){
[24626]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
[25539]155 for(Object* & object : elements->objects){
156 Element* element=xDynamicCast<Element*>(object);
[25379]157 element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs,iomodel,SmbPrecipitationsReconstructedEnum);
[24626]158 }
159 xDelete<int>(ids);
160 iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
[22852]161 }
[19528]162 }
[22448]163 if(issetpddfac){
[25379]164 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
165 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
[22448]166 }
[19528]167 break;
168 case SMBgradientsEnum:
[25379]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);
[19528]173 break;
[21469]174 case SMBgradientselaEnum:
[25379]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);
[21469]180 break;
[19528]181 case SMBhenningEnum:
[25379]182 iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum,0.);
[19528]183 break;
184 case SMBcomponentsEnum:
[25379]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.);
[19528]188 break;
189 case SMBmeltcomponentsEnum:
[25379]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.);
[19528]194 break;
[23366]195 case SMBgradientscomponentsEnum:
[24385]196 /* Nothing to add to input */
[23366]197 break;
[23540]198 case SMBsemicEnum:
[25379]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);
[23540]210 break;
[19528]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;
[23317]220 bool isdelta18o,ismungsm,isd18opd,issetpddfac,interp,isfirnwarming;
[24791]221 int smb_model, smbslices, averaging;
[19528]222 IssmDouble *temp = NULL;
223 int N,M;
[23066]224
[20690]225 parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
[23066]226
[20690]227 iomodel->FindConstant(&smb_model,"md.smb.model");
228 iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
[23066]229
[24791]230 iomodel->FindConstant(&smbslices,"md.smb.steps_per_step");
[24240]231 parameters->AddObject(new IntParam(SmbStepsPerStepEnum,smbslices));
232
[24791]233
234 parameters->AddObject(iomodel->CopyConstantObject("md.smb.averaging",SmbAveragingEnum));
235
[19528]236 switch(smb_model){
237 case SMBforcingEnum:
[23814]238 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
[19528]239 break;
[19554]240 case SMBgembEnum:
[20690]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));
[23468]244 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dsnowIdx",SmbDsnowIdxEnum));
[20690]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));
[23808]260 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
[24735]261 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isconstrainsurfaceT",SmbIsconstrainsurfaceTEnum));
[20690]262 parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
[22475]263 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
[22482]264 parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));
[19554]265 break;
[19528]266 case SMBpddEnum:
[20690]267 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdelta18o",SmbIsdelta18oEnum));
268 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
[22448]269 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
[20690]270 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
271 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
[19528]272
273 if(ismungsm){
[20690]274 iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2);
[19528]275 parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,M));
[20690]276 iomodel->DeleteData(temp,"md.smb.Pfac");
[23066]277
[20690]278 iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
[19528]279 parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,M));
[20690]280 iomodel->DeleteData(temp,"md.smb.Tdiff");
[19528]281
[20690]282 iomodel->FetchData(&temp,&N,&M,"md.smb.sealev"); _assert_(N==2);
[19528]283 parameters->AddObject(new TransientParam(SmbSealevEnum,&temp[0],&temp[M],interp,M));
[20690]284 iomodel->DeleteData(temp,"md.smb.sealev");
[19528]285 }
286 if(isdelta18o){
[20690]287 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
[19528]288 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
[20690]289 iomodel->DeleteData(temp,"md.smb.delta18o");
[19528]290
[20690]291 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o_surface"); _assert_(N==2);
[19528]292 parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,M));
[20690]293 iomodel->DeleteData(temp,"md.smb.delta18o_surface");
[19528]294 }
[24626]295
[19528]296 break;
[23317]297 case SMBpddSicopolisEnum:
[23328]298 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isfirnwarming",SmbIsfirnwarmingEnum));
[23317]299 break;
[19528]300 case SMBd18opddEnum:
[20690]301 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
302 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isd18opd",SmbIsd18opdEnum));
[22448]303 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
[20690]304 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
305 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
[22448]306 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
[19528]307 if(isd18opd){
[22608]308 parameters->AddObject(iomodel->CopyConstantObject("md.smb.f",SmbFEnum));
[22495]309 parameters->AddObject(iomodel->CopyConstantObject("md.smb.istemperaturescaled",SmbIstemperaturescaledEnum));
[22852]310 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipscaled",SmbIsprecipscaledEnum));
[20690]311 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
[19528]312 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
[20690]313 iomodel->DeleteData(temp,"md.smb.delta18o");
[24626]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 }
[19528]370 }
[24626]371
[19528]372 break;
373 case SMBgradientsEnum:
374 /*Nothing to add to parameters*/
375 break;
[21469]376 case SMBgradientselaEnum:
377 /*Nothing to add to parameters*/
378 break;
[19528]379 case SMBhenningEnum:
380 /*Nothing to add to parameters*/
381 break;
382 case SMBcomponentsEnum:
[23814]383 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
[19528]384 break;
385 case SMBmeltcomponentsEnum:
[23814]386 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
[19528]387 break;
[23366]388 case SMBgradientscomponentsEnum:
389 parameters->AddObject(iomodel->CopyConstantObject("md.smb.accualti",SmbAccualtiEnum));
390 parameters->AddObject(iomodel->CopyConstantObject("md.smb.runoffalti",SmbRunoffaltiEnum));
391
[25209]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");
[23366]398
[25209]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");
[23366]405 break;
[23540]406 case SMBsemicEnum:
407 /*Nothing to add to parameters*/
408 break;
[19528]409 default:
410 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
411 }
412
[20690]413 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.smb.requested_outputs");
[19528]414 parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
415 if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
[20690]416 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.smb.requested_outputs");
[19528]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);
[23066]427
[19528]428 /*branch to correct module*/
429 switch(smb_model){
430 case SMBforcingEnum:
[23814]431 SmbForcingx(femmodel);
[19528]432 break;
[19554]433 case SMBgembEnum:
434 Gembx(femmodel);
435 break;
[19528]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);
[23366]443 }
[19528]444 if(ismungsm){
445 if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n");
446 MungsmtpParameterizationx(femmodel);
[23366]447 }
[19528]448 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
449 PositiveDegreeDayx(femmodel);
450 break;
[23317]451 case SMBpddSicopolisEnum:
452 if(VerboseSolution()) _printf0_(" call SICOPOLIS positive degree day module\n");
453 PositiveDegreeDaySicopolisx(femmodel);
454 break;
[19528]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);
[23317]463 }
[19528]464 break;
465 case SMBgradientsEnum:
466 if(VerboseSolution())_printf0_(" call smb gradients module\n");
467 SmbGradientsx(femmodel);
468 break;
[21469]469 case SMBgradientselaEnum:
470 if(VerboseSolution())_printf0_(" call smb gradients ela module\n");
471 SmbGradientsElax(femmodel);
472 break;
[19528]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;
[23366]488 case SMBgradientscomponentsEnum:
489 if(VerboseSolution())_printf0_(" call smb gradients components module\n");
490 SmbGradientsComponentsx(femmodel);
491 break;
[23540]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;
[19528]500 default:
501 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
502 }
503
504}/*}}}*/
505ElementVector* SmbAnalysis::CreateDVector(Element* element){/*{{{*/
506 _error_("not implemented");
507}/*}}}*/
508ElementMatrix* SmbAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
509_error_("Not implemented");
510}/*}}}*/
511ElementMatrix* SmbAnalysis::CreateKMatrix(Element* element){/*{{{*/
512 _error_("not implemented yet");
513}/*}}}*/
514ElementVector* SmbAnalysis::CreatePVector(Element* element){/*{{{*/
515_error_("not implemented yet");
516}/*}}}*/
517void SmbAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
518 _error_("not implemented yet");
519}/*}}}*/
[25317]520void SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
[19528]521 _error_("Not implemented yet");
522}/*}}}*/
523void SmbAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
524 _error_("not implemented yet");
525}/*}}}*/
526void SmbAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
527 /*Default, do nothing*/
528 return;
529}/*}}}*/
Note: See TracBrowser for help on using the repository browser.