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

Last change on this file since 27847 was 27847, checked in by rueckamp, 2 years ago

CHG: Several imprvements to Debris analysis
File size: 31.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;
[26750]26 bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled,isfirnwarming,isstochastic;
[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");
[26750]40 iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing");
[19528]41 switch(smb_model){
42 case SMBforcingEnum:
[25379]43 iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
[26750]44 if(isstochastic){
45 iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
46 iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",BaselineSmbMassBalanceEnum,0.);
47 }
[19528]48 break;
[19554]49 case SMBgembEnum:
[25379]50 iomodel->FetchDataToInput(inputs,elements,"md.smb.Ta",SmbTaEnum);
51 iomodel->FetchDataToInput(inputs,elements,"md.smb.V",SmbVEnum);
52 iomodel->FetchDataToInput(inputs,elements,"md.smb.dswrf",SmbDswrfEnum);
[25997]53 iomodel->FetchDataToInput(inputs,elements,"md.smb.dswdiffrf",SmbDswdiffrfEnum);
[25379]54 iomodel->FetchDataToInput(inputs,elements,"md.smb.dlwrf",SmbDlwrfEnum);
55 iomodel->FetchDataToInput(inputs,elements,"md.smb.P",SmbPEnum);
56 iomodel->FetchDataToInput(inputs,elements,"md.smb.eAir",SmbEAirEnum);
57 iomodel->FetchDataToInput(inputs,elements,"md.smb.pAir",SmbPAirEnum);
58 iomodel->FetchDataToInput(inputs,elements,"md.smb.zTop",SmbZTopEnum);
59 iomodel->FetchDataToInput(inputs,elements,"md.smb.dzTop",SmbDzTopEnum);
60 iomodel->FetchDataToInput(inputs,elements,"md.smb.dzMin",SmbDzMinEnum);
61 iomodel->FetchDataToInput(inputs,elements,"md.smb.zY",SmbZYEnum);
62 iomodel->FetchDataToInput(inputs,elements,"md.smb.zMax",SmbZMaxEnum);
63 iomodel->FetchDataToInput(inputs,elements,"md.smb.zMin",SmbZMinEnum);
64 iomodel->FetchDataToInput(inputs,elements,"md.smb.Tmean",SmbTmeanEnum);
65 iomodel->FetchDataToInput(inputs,elements,"md.smb.Vmean",SmbVmeanEnum);
66 iomodel->FetchDataToInput(inputs,elements,"md.smb.C",SmbCEnum);
67 iomodel->FetchDataToInput(inputs,elements,"md.smb.Tz",SmbTzEnum);
68 iomodel->FetchDataToInput(inputs,elements,"md.smb.Vz",SmbVzEnum);
69 InputUpdateFromConstantx(inputs,elements,false,SmbIsInitializedEnum);
70 iomodel->FetchDataToInput(inputs,elements,"md.smb.Dzini",SmbDziniEnum);
71 iomodel->FetchDataToInput(inputs,elements,"md.smb.Dini",SmbDiniEnum);
72 iomodel->FetchDataToInput(inputs,elements,"md.smb.Reini",SmbReiniEnum);
73 iomodel->FetchDataToInput(inputs,elements,"md.smb.Gdnini",SmbGdniniEnum);
74 iomodel->FetchDataToInput(inputs,elements,"md.smb.Gspini",SmbGspiniEnum);
75 iomodel->FetchDataToInput(inputs,elements,"md.smb.ECini",SmbECiniEnum);
76 iomodel->FetchDataToInput(inputs,elements,"md.smb.Wini",SmbWiniEnum);
77 iomodel->FetchDataToInput(inputs,elements,"md.smb.Aini",SmbAiniEnum);
[25997]78 iomodel->FetchDataToInput(inputs,elements,"md.smb.Adiffini",SmbAdiffiniEnum);
[25379]79 iomodel->FetchDataToInput(inputs,elements,"md.smb.Tini",SmbTiniEnum);
80 iomodel->FetchDataToInput(inputs,elements,"md.smb.Sizeini",SmbSizeiniEnum);
81 iomodel->FetchDataToInput(inputs,elements,"md.smb.aValue",SmbAValueEnum);
[26550]82 iomodel->FetchDataToInput(inputs,elements,"md.smb.dulwrfValue",SmbDulwrfValueEnum);
[25379]83 iomodel->FetchDataToInput(inputs,elements,"md.smb.teValue",SmbTeValueEnum);
[25997]84 iomodel->FetchDataToInput(inputs,elements,"md.smb.szaValue",SmbSzaValueEnum);
85 iomodel->FetchDataToInput(inputs,elements,"md.smb.cotValue",SmbCotValueEnum);
86 iomodel->FetchDataToInput(inputs,elements,"md.smb.ccsnowValue",SmbCcsnowValueEnum);
87 iomodel->FetchDataToInput(inputs,elements,"md.smb.cciceValue",SmbCciceValueEnum);
[19554]88 break;
[19528]89 case SMBpddEnum:
[20690]90 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
91 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
[25374]92 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
[25379]93 iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
94 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
95 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
[19528]96 if(isdelta18o || ismungsm){
[25379]97 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_lgm",SmbTemperaturesLgmEnum);
98 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
99 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
100 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum);
[22852]101 }else{
[25379]102 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitation",SmbPrecipitationEnum);
103 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
[19528]104 }
[25374]105 if(issetpddfac){
[25379]106 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
107 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
[25374]108 }
[19528]109 break;
[23317]110 case SMBpddSicopolisEnum:
[25379]111 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
112 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
[23317]113 iomodel->FindConstant(&isfirnwarming,"md.smb.isfirnwarming");
[25379]114 iomodel->FetchDataToInput(inputs,elements,"md.smb.smb_corr",SmbSmbCorrEnum);
115 iomodel->FetchDataToInput(inputs,elements,"md.smb.precipitation_anomaly",SmbPrecipitationsAnomalyEnum);
116 iomodel->FetchDataToInput(inputs,elements,"md.smb.temperature_anomaly",SmbTemperaturesAnomalyEnum);
117 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
118 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitation",SmbPrecipitationEnum);
[23317]119 break;
[19528]120 case SMBd18opddEnum:
[22852]121 iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
122 iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
[20690]123 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
124 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
[22448]125 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
[25379]126 iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
127 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
128 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
[20690]129 if(isd18opd){
[25379]130 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
131 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
[22852]132 if(!istemperaturescaled){
[24626]133 /*Fetch array*/
134 IssmDouble* doublearray = NULL;
135 int M,N;
136 iomodel->FetchData(&doublearray,&M,&N,"md.smb.temperatures_reconstructed");
137 if(M!=iomodel->numberofvertices+1) _error_("md.smb.temperatures_reconstructed should have nbv+1 rows");
138 if(N%12!=0) _error_("md.smb.temperatures_reconstructed should have a multiple of 12 columns (since it is monthly)");
139
140 /*Build indices*/
141 int* ids = xNew<int>(N); for(int i=0;i<N;i++) ids[i] = i;
142
[25539]143 for(Object* & object : elements->objects){
144 Element* element=xDynamicCast<Element*>(object);
[25379]145 element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs,iomodel,SmbTemperaturesReconstructedEnum);
[24626]146 }
147 xDelete<int>(ids);
148 iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
[22852]149 }
150 if(!isprecipscaled){
[24626]151 /*Fetch array*/
152 IssmDouble* doublearray = NULL;
153 int M,N;
154 iomodel->FetchData(&doublearray,&M,&N,"md.smb.precipitations_reconstructed");
155 if(M!=iomodel->numberofvertices+1) _error_("md.smb.precipitations_reconstructed should have nbv+1 rows");
156 if(N%12!=0) _error_("md.smb.precipitations_reconstructed should have a multiple of 12 columns (since it is monthly)");
157
158 /*Build indices*/
159 int* ids = xNew<int>(N); for(int i=0;i<N;i++) ids[i] = i;
160
[25539]161 for(Object* & object : elements->objects){
162 Element* element=xDynamicCast<Element*>(object);
[25379]163 element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs,iomodel,SmbPrecipitationsReconstructedEnum);
[24626]164 }
165 xDelete<int>(ids);
166 iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
[22852]167 }
[19528]168 }
[22448]169 if(issetpddfac){
[25379]170 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
171 iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
[22448]172 }
[19528]173 break;
174 case SMBgradientsEnum:
[25379]175 iomodel->FetchDataToInput(inputs,elements,"md.smb.href",SmbHrefEnum);
176 iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum);
177 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_pos",SmbBPosEnum);
178 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_neg",SmbBNegEnum);
[19528]179 break;
[21469]180 case SMBgradientselaEnum:
[25379]181 iomodel->FetchDataToInput(inputs,elements,"md.smb.ela",SmbElaEnum);
182 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_pos",SmbBPosEnum);
183 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_neg",SmbBNegEnum);
184 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_max",SmbBMaxEnum);
185 iomodel->FetchDataToInput(inputs,elements,"md.smb.b_min",SmbBMinEnum);
[21469]186 break;
[27250]187 case SMBarmaEnum:
[26477]188 iomodel->FetchDataToInput(inputs,elements,"md.smb.basin_id",SmbBasinsIdEnum);
189 break;
[19528]190 case SMBhenningEnum:
[25379]191 iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum,0.);
[19528]192 break;
193 case SMBcomponentsEnum:
[25379]194 iomodel->FetchDataToInput(inputs,elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
195 iomodel->FetchDataToInput(inputs,elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
196 iomodel->FetchDataToInput(inputs,elements,"md.smb.runoff",SmbRunoffEnum,0.);
[19528]197 break;
198 case SMBmeltcomponentsEnum:
[25379]199 iomodel->FetchDataToInput(inputs,elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
200 iomodel->FetchDataToInput(inputs,elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
201 iomodel->FetchDataToInput(inputs,elements,"md.smb.melt",SmbMeltEnum,0.);
202 iomodel->FetchDataToInput(inputs,elements,"md.smb.refreeze",SmbRefreezeEnum,0.);
[19528]203 break;
[23366]204 case SMBgradientscomponentsEnum:
[24385]205 /* Nothing to add to input */
[23366]206 break;
[23540]207 case SMBsemicEnum:
[27498]208 int ismethod;
209 //if(VerboseSolution()) _printf0_(" smb semic: UpdateElements.\n");
210 //iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
[25379]211 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0gcm",SmbS0gcmEnum);
212 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailysnowfall",SmbDailysnowfallEnum);
213 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyrainfall",SmbDailyrainfallEnum);
214 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailydsradiation",SmbDailydsradiationEnum);
215 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailydlradiation",SmbDailydlradiationEnum);
216 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailywindspeed",SmbDailywindspeedEnum);
217 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailypressure",SmbDailypressureEnum);
218 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyairdensity",SmbDailyairdensityEnum);
219 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyairhumidity",SmbDailyairhumidityEnum);
220 iomodel->FetchDataToInput(inputs,elements,"md.smb.dailytemperature",SmbDailytemperatureEnum);
[27498]221 // assign initial SEMIC temperature from initialization class.
222 if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - temperature.\n");
[27510]223 iomodel->FetchDataToInput(inputs,elements,"md.initialization.temperature",TemperatureEnum);
[27498]224
225 iomodel->FindConstant(&ismethod,"md.smb.ismethod");
226 if (ismethod == 1){
[27646]227 if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - albedo.\n");
[27510]228 iomodel->FetchDataToInput(inputs,elements,"md.smb.albedo",SmbAlbedoInitEnum);
229 iomodel->FetchDataToInput(inputs,elements,"md.smb.albedo_snow",SmbAlbedoSnowInitEnum);
[27498]230 if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - Hice/Hsnow.\n");
[27510]231 iomodel->FetchDataToInput(inputs,elements,"md.smb.hice",SmbHIceInitEnum);
232 iomodel->FetchDataToInput(inputs,elements,"md.smb.hsnow",SmbHSnowInitEnum);
[27498]233
234 // initial Temperature amplitude.
235 if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - Tamp.\n");
236 iomodel->FetchDataToInput(inputs,elements,"md.smb.Tamp",SmbTampEnum);
237
238 // assign masking
239 iomodel->FetchDataToInput(inputs,elements,"md.smb.mask",SmbMaskEnum);
[27646]240 iomodel->FetchDataToInput(inputs,elements,"md.smb.qmr",SmbSemicQmrInitEnum);
[27498]241 if(VerboseSolution()) _printf0_(" smb semic: UpdateElements - done.\n");
242 }
[23540]243 break;
[27847]244 case SMBdebrisEvattEnum:
245 iomodel->FetchDataToInput(inputs,elements,"md.initialization.debris",DebrisThicknessEnum);
246 iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
247 iomodel->FetchDataToInput(inputs,elements,"md.smb.snowheight",SmbSnowheightEnum);
248 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
249 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlyprecipitation",SmbPrecipitationEnum);
250 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlydsradiation",SmbMonthlydsradiationEnum);
251 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlydlradiation",SmbMonthlydlradiationEnum);
252 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlywindspeed",SmbMonthlywindspeedEnum);
253 iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlyairhumidity",SmbMonthlyairhumidityEnum);
254 iomodel->FetchDataToInput(inputs,elements,"md.smb.precipitation_anomaly",SmbPrecipitationsAnomalyEnum);
255 iomodel->FetchDataToInput(inputs,elements,"md.smb.temperature_anomaly",SmbTemperaturesAnomalyEnum);
256 iomodel->FetchDataToInput(inputs,elements,"md.smb.dsradiation_anomaly",SmbDsradiationAnomalyEnum);
257 iomodel->FetchDataToInput(inputs,elements,"md.smb.dlradiation_anomaly",SmbDlradiationAnomalyEnum);
258 iomodel->FetchDataToInput(inputs,elements,"md.smb.windspeed_anomaly",SmbWindspeedAnomalyEnum);
259 iomodel->FetchDataToInput(inputs,elements,"md.smb.airhumidity_anomaly",SmbAirhumidityAnomalyEnum);
260 break;
[19528]261 default:
262 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
263 }
264
265}/*}}}*/
266void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
267
268 int numoutputs;
269 char** requestedoutputs = NULL;
[26196]270 bool isdelta18o,ismungsm,isd18opd,issetpddfac,interp,cycle,isfirnwarming;
[24791]271 int smb_model, smbslices, averaging;
[19528]272 IssmDouble *temp = NULL;
273 int N,M;
[23066]274
[20690]275 parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
[23066]276
[20690]277 iomodel->FindConstant(&smb_model,"md.smb.model");
[26208]278 iomodel->FindConstant(&interp,"md.timestepping.interp_forcing");
[26196]279 iomodel->FindConstant(&cycle,"md.timestepping.cycle_forcing");
[23066]280
[24791]281 iomodel->FindConstant(&smbslices,"md.smb.steps_per_step");
[24240]282 parameters->AddObject(new IntParam(SmbStepsPerStepEnum,smbslices));
283
[24791]284 parameters->AddObject(iomodel->CopyConstantObject("md.smb.averaging",SmbAveragingEnum));
285
[19528]286 switch(smb_model){
287 case SMBforcingEnum:
288 break;
[19554]289 case SMBgembEnum:
[20690]290 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIdx",SmbAIdxEnum));
[26550]291 parameters->AddObject(iomodel->CopyConstantObject("md.smb.eIdx",SmbEIdxEnum));
[27229]292 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tcIdx",SmbTcIdxEnum));
[20690]293 parameters->AddObject(iomodel->CopyConstantObject("md.smb.swIdx",SmbSwIdxEnum));
294 parameters->AddObject(iomodel->CopyConstantObject("md.smb.denIdx",SmbDenIdxEnum));
[23468]295 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dsnowIdx",SmbDsnowIdxEnum));
[20690]296 parameters->AddObject(iomodel->CopyConstantObject("md.smb.cldFrac",SmbCldFracEnum));
297 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0wet",SmbT0wetEnum));
298 parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0dry",SmbT0dryEnum));
299 parameters->AddObject(iomodel->CopyConstantObject("md.smb.K",SmbKEnum));
300 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aSnow",SmbASnowEnum));
301 parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));
302 parameters->AddObject(iomodel->CopyConstantObject("md.smb.dt",SmbDtEnum));
303 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isgraingrowth",SmbIsgraingrowthEnum));
304 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isalbedo",SmbIsalbedoEnum));
305 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isshortwave",SmbIsshortwaveEnum));
306 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isthermal",SmbIsthermalEnum));
307 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isaccumulation",SmbIsaccumulationEnum));
308 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismelt",SmbIsmeltEnum));
309 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdensification",SmbIsdensificationEnum));
310 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum));
[24735]311 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isconstrainsurfaceT",SmbIsconstrainsurfaceTEnum));
[26550]312 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdeltaLWup",SmbIsdeltaLWupEnum));
[20690]313 parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
[22475]314 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
[22482]315 parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));
[26550]316 parameters->AddObject(iomodel->CopyConstantObject("md.smb.teThresh",SmbTeThreshEnum));
[19554]317 break;
[19528]318 case SMBpddEnum:
[20690]319 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdelta18o",SmbIsdelta18oEnum));
320 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
[22448]321 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
[20690]322 iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
323 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
[19528]324
325 if(ismungsm){
[20690]326 iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2);
[26196]327 parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,cycle,M));
[20690]328 iomodel->DeleteData(temp,"md.smb.Pfac");
[23066]329
[20690]330 iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
[26196]331 parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,cycle,M));
[20690]332 iomodel->DeleteData(temp,"md.smb.Tdiff");
[19528]333
[20690]334 iomodel->FetchData(&temp,&N,&M,"md.smb.sealev"); _assert_(N==2);
[26196]335 parameters->AddObject(new TransientParam(SmbSealevEnum,&temp[0],&temp[M],interp,cycle,M));
[20690]336 iomodel->DeleteData(temp,"md.smb.sealev");
[19528]337 }
338 if(isdelta18o){
[20690]339 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
[26196]340 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,cycle,M));
[20690]341 iomodel->DeleteData(temp,"md.smb.delta18o");
[19528]342
[20690]343 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o_surface"); _assert_(N==2);
[26196]344 parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,cycle,M));
[20690]345 iomodel->DeleteData(temp,"md.smb.delta18o_surface");
[19528]346 }
[24626]347
[19528]348 break;
[23317]349 case SMBpddSicopolisEnum:
[23328]350 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isfirnwarming",SmbIsfirnwarmingEnum));
[23317]351 break;
[19528]352 case SMBd18opddEnum:
[20690]353 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
354 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isd18opd",SmbIsd18opdEnum));
[22448]355 parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
[20690]356 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
357 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
[22448]358 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
[19528]359 if(isd18opd){
[22608]360 parameters->AddObject(iomodel->CopyConstantObject("md.smb.f",SmbFEnum));
[22495]361 parameters->AddObject(iomodel->CopyConstantObject("md.smb.istemperaturescaled",SmbIstemperaturescaledEnum));
[22852]362 parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipscaled",SmbIsprecipscaledEnum));
[20690]363 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
[26196]364 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,cycle,M));
[20690]365 iomodel->DeleteData(temp,"md.smb.delta18o");
[24626]366
367 IssmDouble yts;
368 bool istemperaturescaled,isprecipscaled;
369 iomodel->FindConstant(&yts,"md.constants.yts");
370 iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
371 iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
372 if(!istemperaturescaled){
373 /*Fetch array*/
374 IssmDouble* doublearray = NULL;
375 int M,N;
376 iomodel->FetchData(&doublearray,&M,&N,"md.smb.temperatures_reconstructed");
377 if(M!=iomodel->numberofvertices+1) _error_("md.smb.temperatures_reconstructed should have nbv+1 rows");
378 if(N%12!=0) _error_("md.smb.temperatures_reconstructed should have a multiple of 12 columns (since it is monthly)");
379 int numyears = N/12; _assert_(numyears*12==N);
380
381 /*Check times*/
382 #ifdef _ISSM_DEBUG_
383 for(int i=0;i<numyears;i++){
384 for(int j=1;j<12;j++){
385 //_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
386 _assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
387 }
388 }
389 #endif
390
391 /*Build time*/
392 IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
393 parameters->AddObject(new DoubleVecParam(SmbTemperaturesReconstructedYearsEnum,times,numyears));
394 xDelete<IssmDouble>(times);
395 iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
396 }
397 if(!isprecipscaled){
398 /*Fetch array*/
399 IssmDouble* doublearray = NULL;
400 int M,N;
401 iomodel->FetchData(&doublearray,&M,&N,"md.smb.precipitations_reconstructed");
402 if(M!=iomodel->numberofvertices+1) _error_("md.smb.precipitations_reconstructed should have nbv+1 rows");
403 if(N%12!=0) _error_("md.smb.precipitations_reconstructed should have a multiple of 12 columns (since it is monthly)");
404 int numyears = N/12; _assert_(numyears*12==N);
405
406 /*Check times*/
407 #ifdef _ISSM_DEBUG_
408 for(int i=0;i<numyears;i++){
409 for(int j=1;j<12;j++){
410 //_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
411 _assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
412 }
413 }
414 #endif
415
416 /*Build time*/
417 IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
418 parameters->AddObject(new DoubleVecParam(SmbPrecipitationsReconstructedYearsEnum,times,numyears));
419 xDelete<IssmDouble>(times);
420 iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
421 }
[19528]422 }
[24626]423
[19528]424 break;
425 case SMBgradientsEnum:
426 /*Nothing to add to parameters*/
427 break;
[21469]428 case SMBgradientselaEnum:
429 /*Nothing to add to parameters*/
430 break;
[27250]431 case SMBarmaEnum:
[26477]432 /*Nothing to add to parameters*/
433 break;
[19528]434 case SMBhenningEnum:
435 /*Nothing to add to parameters*/
436 break;
437 case SMBcomponentsEnum:
438 break;
439 case SMBmeltcomponentsEnum:
440 break;
[23366]441 case SMBgradientscomponentsEnum:
442 parameters->AddObject(iomodel->CopyConstantObject("md.smb.accualti",SmbAccualtiEnum));
443 parameters->AddObject(iomodel->CopyConstantObject("md.smb.runoffalti",SmbRunoffaltiEnum));
444
[25209]445 iomodel->FetchData(&temp,&N,&M,"md.smb.accugrad"); _assert_(N==2);
[26196]446 parameters->AddObject(new TransientParam(SmbAccugradEnum,&temp[0],&temp[M],interp,cycle,M));
[25209]447 iomodel->DeleteData(temp,"md.smb.accugrad");
448 iomodel->FetchData(&temp,&N,&M,"md.smb.runoffgrad"); _assert_(N==2);
[26196]449 parameters->AddObject(new TransientParam(SmbRunoffgradEnum,&temp[0],&temp[M],interp,cycle,M));
[25209]450 iomodel->DeleteData(temp,"md.smb.runoffgrad");
[23366]451
[25209]452 iomodel->FetchData(&temp,&N,&M,"md.smb.accuref"); _assert_(N==2);
[26196]453 parameters->AddObject(new TransientParam(SmbAccurefEnum,&temp[0],&temp[M],interp,cycle,M));
[25209]454 iomodel->DeleteData(temp,"md.smb.accuref");
455 iomodel->FetchData(&temp,&N,&M,"md.smb.runoffref"); _assert_(N==2);
[26196]456 parameters->AddObject(new TransientParam(SmbRunoffrefEnum,&temp[0],&temp[M],interp,cycle,M));
[25209]457 iomodel->DeleteData(temp,"md.smb.runoffref");
[23366]458 break;
[23540]459 case SMBsemicEnum:
[27498]460 int ismethod;
461 parameters->FindParam(&ismethod,SmbSemicMethodEnum);
462 if (ismethod == 1){
463 parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfacElevation",SmbDesfacElevEnum));
464 parameters->AddObject(iomodel->CopyConstantObject("md.smb.hcrit",SmbSemicHcritEnum));
465 parameters->AddObject(iomodel->CopyConstantObject("md.smb.rcrit",SmbSemicRcritEnum));
466 /*Define albedo parameters.*/
467 parameters->AddObject(iomodel->CopyConstantObject("md.smb.albedo_scheme",SmbAlbedoSchemeEnum));
468 parameters->AddObject(iomodel->CopyConstantObject("md.smb.alb_smax",SmbAlbedoSnowMaxEnum));
469 parameters->AddObject(iomodel->CopyConstantObject("md.smb.alb_smin",SmbAlbedoSnowMinEnum));
470 parameters->AddObject(iomodel->CopyConstantObject("md.smb.albi",SmbAlbedoIceEnum));
471 parameters->AddObject(iomodel->CopyConstantObject("md.smb.albl",SmbAlbedoLandEnum));
[27510]472
473 //albedo parameter - slatter
474 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tmin",SmbSemicTmaxEnum));
475 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tmax",SmbSemicTminEnum));
476
477 //albedo parameter - isba & denby
478 parameters->AddObject(iomodel->CopyConstantObject("md.smb.mcrit",SmbSemicMcritEnum));
479 //albedo parameter - isba
480 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tau_a",SmbSemicTauAEnum));
481 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tau_f",SmbSemicTauFEnum));
482 parameters->AddObject(iomodel->CopyConstantObject("md.smb.wcrit",SmbSemicWcritEnum));
483 //albedo parameter - alex
484 parameters->AddObject(iomodel->CopyConstantObject("md.smb.tmid",SmbSemicTmidEnum));
485 parameters->AddObject(iomodel->CopyConstantObject("md.smb.afac",SmbSemicAfacEnum));
[27498]486 }
[23540]487 /*Nothing to add to parameters*/
488 break;
[27847]489 case SMBdebrisEvattEnum:
490 /*Nothing to add to parameters*/
491 break;
[19528]492 default:
493 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
494 }
495
[20690]496 iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.smb.requested_outputs");
[19528]497 parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
498 if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
[20690]499 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.smb.requested_outputs");
[19528]500
501}/*}}}*/
502
503/*Finite Element Analysis*/
504void SmbAnalysis::Core(FemModel* femmodel){/*{{{*/
505
506 int smb_model;
507
508 /*Figure out smb model: */
509 femmodel->parameters->FindParam(&smb_model,SmbEnum);
[23066]510
[19528]511 /*branch to correct module*/
512 switch(smb_model){
513 case SMBforcingEnum:
[23814]514 SmbForcingx(femmodel);
[19528]515 break;
[19554]516 case SMBgembEnum:
517 Gembx(femmodel);
518 break;
[19528]519 case SMBpddEnum:
520 bool isdelta18o,ismungsm;
521 femmodel->parameters->FindParam(&isdelta18o,SmbIsdelta18oEnum);
522 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
523 if(isdelta18o){
524 if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n");
525 Delta18oParameterizationx(femmodel);
[23366]526 }
[19528]527 if(ismungsm){
528 if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n");
529 MungsmtpParameterizationx(femmodel);
[23366]530 }
[19528]531 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
532 PositiveDegreeDayx(femmodel);
533 break;
[23317]534 case SMBpddSicopolisEnum:
535 if(VerboseSolution()) _printf0_(" call SICOPOLIS positive degree day module\n");
536 PositiveDegreeDaySicopolisx(femmodel);
537 break;
[19528]538 case SMBd18opddEnum:
539 bool isd18opd;
540 femmodel->parameters->FindParam(&isd18opd,SmbIsd18opdEnum);
541 if(isd18opd){
542 if(VerboseSolution()) _printf0_(" call Delta18opdParameterization module\n");
543 Delta18opdParameterizationx(femmodel);
544 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
545 PositiveDegreeDayx(femmodel);
[23317]546 }
[19528]547 break;
548 case SMBgradientsEnum:
549 if(VerboseSolution())_printf0_(" call smb gradients module\n");
550 SmbGradientsx(femmodel);
551 break;
[21469]552 case SMBgradientselaEnum:
553 if(VerboseSolution())_printf0_(" call smb gradients ela module\n");
554 SmbGradientsElax(femmodel);
555 break;
[27250]556 case SMBarmaEnum:
[27297]557 if(VerboseSolution())_printf0_(" call smb arma module\n");
558 Smbarmax(femmodel);
559 break;
[19528]560 case SMBhenningEnum:
561 if(VerboseSolution())_printf0_(" call smb Henning module\n");
562 SmbHenningx(femmodel);
563 break;
564 case SMBcomponentsEnum:
565 if(VerboseSolution())_printf0_(" call smb Components module\n");
566 SmbComponentsx(femmodel);
567 break;
568 case SMBmeltcomponentsEnum:
569 if(VerboseSolution())_printf0_(" call smb Melt Components module\n");
570 SmbMeltComponentsx(femmodel);
571 break;
572 case SMBgcmEnum:
573 /*Nothing to be done*/
574 break;
[23366]575 case SMBgradientscomponentsEnum:
576 if(VerboseSolution())_printf0_(" call smb gradients components module\n");
577 SmbGradientsComponentsx(femmodel);
578 break;
[23540]579 case SMBsemicEnum:
580 #ifdef _HAVE_SEMIC_
[27498]581 if(VerboseSolution())_printf0_(" call smb SEMIC module\n");
582 int ismethod;
583 femmodel->parameters->FindParam(&ismethod,SmbSemicMethodEnum);
584 SmbSemicx(femmodel,ismethod);
[23540]585 #else
586 _error_("SEMIC not installed");
587 #endif //_HAVE_SEMIC_
588 break;
[27847]589 case SMBdebrisEvattEnum:
590 if(VerboseSolution())_printf0_(" call smb Evatt debris module\n");
591 SmbDebrisEvattx(femmodel);
592 break;
[19528]593 default:
594 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
595 }
596
597}/*}}}*/
[26047]598void SmbAnalysis::PreCore(FemModel* femmodel){/*{{{*/
599 _error_("not implemented");
600}/*}}}*/
[19528]601ElementVector* SmbAnalysis::CreateDVector(Element* element){/*{{{*/
602 _error_("not implemented");
603}/*}}}*/
604ElementMatrix* SmbAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
605_error_("Not implemented");
606}/*}}}*/
607ElementMatrix* SmbAnalysis::CreateKMatrix(Element* element){/*{{{*/
608 _error_("not implemented yet");
609}/*}}}*/
610ElementVector* SmbAnalysis::CreatePVector(Element* element){/*{{{*/
611_error_("not implemented yet");
612}/*}}}*/
613void SmbAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
614 _error_("not implemented yet");
615}/*}}}*/
[25317]616void SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/
[19528]617 _error_("Not implemented yet");
618}/*}}}*/
619void SmbAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
620 _error_("not implemented yet");
621}/*}}}*/
622void SmbAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
623 /*Default, do nothing*/
624 return;
625}/*}}}*/
Note: See TracBrowser for help on using the repository browser.