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

Last change on this file since 25374 was 25374, checked in by schlegel, 5 years ago

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