Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 24582) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 24583) @@ -82,11 +82,11 @@ iomodel->FetchDataToInput(inputs2,elements,"md.smb.s0t",SmbS0tEnum); if(isdelta18o || ismungsm){ - iomodel->FetchDataToInput(inputs2,elements,"md.smb.temperatures_lgm",SmbTemperaturesLgmEnum); - iomodel->FetchDataToInput(inputs2,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum); - iomodel->FetchDataToInput(inputs2,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum); - iomodel->FetchDataToInput(inputs2,elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum); + iomodel->FetchDataToDatasetInput(inputs2,elements,"md.smb.temperatures_lgm",SmbTemperaturesLgmEnum); + iomodel->FetchDataToDatasetInput(inputs2,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum); + iomodel->FetchDataToDatasetInput(inputs2,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum); + iomodel->FetchDataToDatasetInput(inputs2,elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum); }else{ - iomodel->FetchDataToInput(inputs2,elements,"md.smb.precipitation",SmbPrecipitationEnum); - iomodel->FetchDataToInput(inputs2,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum); + iomodel->FetchDataToDatasetInput(inputs2,elements,"md.smb.precipitation",SmbPrecipitationEnum); + iomodel->FetchDataToDatasetInput(inputs2,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum); } break; @@ -111,6 +111,6 @@ iomodel->FetchDataToInput(inputs2,elements,"md.smb.s0t",SmbS0tEnum); if(isd18opd){ - iomodel->FetchDataToInput(inputs2,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum); - iomodel->FetchDataToInput(inputs2,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum); + iomodel->FetchDataToDatasetInput(inputs2,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum); + iomodel->FetchDataToDatasetInput(inputs2,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum); if(!istemperaturescaled){ iomodel->FetchDataToInput(inputs2,elements,"md.smb.temperatures_reconstructed",SmbTemperaturesReconstructedEnum); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24582) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24583) @@ -455,18 +455,17 @@ IssmDouble time_yr=floor(time/yts)*yts; + /*Recover present day temperature and precipitation*/ + DatasetInput2* dinput1=this->GetDatasetInput2(SmbTemperaturesPresentdayEnum); _assert_(dinput1); + DatasetInput2* dinput2=this->GetDatasetInput2(SmbTemperaturesLgmEnum); _assert_(dinput2); + DatasetInput2* dinput3=this->GetDatasetInput2(SmbPrecipitationsPresentdayEnum); _assert_(dinput3); + /*loop over vertices: */ Gauss* gauss=this->NewGauss(); for(int month=0;month<12;month++){ - - /*Recover present day temperature and precipitation*/ - Input2* input1=this->GetInput2(SmbTemperaturesPresentdayEnum,month/12.*yts); _assert_(input1); - Input2* input2=this->GetInput2(SmbTemperaturesLgmEnum,month/12.*yts); _assert_(input2); - Input2* input3=this->GetInput2(SmbPrecipitationsPresentdayEnum,month/12.*yts); _assert_(input3); - for(int iv=0;ivGaussVertex(iv); - input1->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss); - input2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss); - input3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss); + dinput1->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month); + dinput2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss,month); + dinput3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month); PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts; @@ -537,5 +536,5 @@ const int NUM_VERTICES_MONTHS_PER_YEAR = NUM_VERTICES * 12; - int i; + int i,offset; int* vertexlids=xNew(NUM_VERTICES); IssmDouble* monthlytemperatures=xNew(NUM_VERTICES_MONTHS_PER_YEAR); @@ -547,8 +546,6 @@ IssmDouble* tmp=xNew(NUM_VERTICES); IssmDouble Delta18oTime; - IssmDouble dpermil,f; - IssmDouble time,yts,time_yr,month,time_clim,time_climt,time_climp,del_clim; - bool isTemperatureScaled=true; - bool isPrecipScaled=true; + IssmDouble f; + IssmDouble time,yts,time_yr,month,time_climt,time_climp,del_clim; this->parameters->FindParam(&time,TimeEnum); this->parameters->FindParam(&yts,ConstantsYtsEnum); @@ -556,24 +553,16 @@ this->GetVerticesLidList(vertexlids); time_yr=floor(time/yts)*yts; - time_clim =ceil(time/yts + 1e-10)*yts; time_climt=ceil(time/yts + 1e-10)*yts; time_climp=ceil(time/yts + 1e-10)*yts; /*Get some pdd parameters*/ - dpermil=this->FindParam(SmbDpermilEnum); - + bool isTemperatureScaled,isPrecipScaled; + IssmDouble dpermil=this->FindParam(SmbDpermilEnum); this->parameters->FindParam(&isTemperatureScaled,SmbIstemperaturescaledEnum); this->parameters->FindParam(&isPrecipScaled,SmbIsprecipscaledEnum); /*Recover present day temperature and precipitation*/ - TransientInput2 *tinput = this->inputs2->GetTransientInput(SmbTemperaturesPresentdayEnum); _assert_(tinput); - TransientInput2 *tinput2 = this->inputs2->GetTransientInput(SmbPrecipitationsPresentdayEnum); _assert_(tinput2); TransientInput2 *tinput3 = NULL; TransientInput2 *tinput4 = NULL; - int offset; - - offset = tinput->GetTimeInputOffset(time_clim); - time_clim = tinput->GetTimeByOffset(offset-11-fmod(offset-11,12.)); - if(!isTemperatureScaled){ tinput3=this->inputs2->GetTransientInput(SmbTemperaturesReconstructedEnum); _assert_(tinput3); @@ -587,9 +576,12 @@ } + /*Get present day temp and precip (monthly)*/ + DatasetInput2 *dinput1 = this->GetDatasetInput2(SmbTemperaturesPresentdayEnum); _assert_(dinput1); + DatasetInput2 *dinput2 = this->GetDatasetInput2(SmbPrecipitationsPresentdayEnum); _assert_(dinput2); + /*loop over vertices: */ Gauss* gauss=this->NewGauss(); for(int month=0;month<12;month++) { - Input2 *input = this->GetInput2(SmbTemperaturesPresentdayEnum,time_clim+month/12.*yts); _assert_(input); - Input2 *input2 = this->GetInput2(SmbPrecipitationsPresentdayEnum,time_clim+month/12.*yts); _assert_(input2); + Input2 *input3 = NULL; Input2 *input4 = NULL; @@ -602,6 +594,6 @@ for(int iv=0;ivGaussVertex(iv); - input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss); - input2->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss); + dinput1->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month); + dinput2->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month); PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts; @@ -2439,18 +2431,19 @@ time_yr=floor(time/yts)*yts; + /*Recover present day temperature and precipitation*/ + DatasetInput2* dinput1=this->GetDatasetInput2(SmbTemperaturesPresentdayEnum); _assert_(dinput1); + DatasetInput2* dinput2=this->GetDatasetInput2(SmbTemperaturesLgmEnum); _assert_(dinput2); + DatasetInput2* dinput3=this->GetDatasetInput2(SmbPrecipitationsPresentdayEnum); _assert_(dinput3); + DatasetInput2* dinput4=this->GetDatasetInput2(SmbPrecipitationsLgmEnum); _assert_(dinput4); + /*loop over vertices: */ Gauss* gauss=this->NewGauss(); for(int month=0;month<12;month++) { - /*Recover present day temperature and precipitation*/ - Input2* input =this->GetInput2(SmbTemperaturesPresentdayEnum,month/12.*yts); _assert_(input); - Input2* input2=this->GetInput2(SmbTemperaturesLgmEnum,month/12.*yts); _assert_(input2); - Input2* input3=this->GetInput2(SmbPrecipitationsPresentdayEnum,month/12.*yts); _assert_(input3); - Input2* input4=this->GetInput2(SmbPrecipitationsLgmEnum,month/12.*yts); _assert_(input4); for(int iv=0;ivGaussVertex(iv); - input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss); - input2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss); - input3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss); - input4->GetInputValue(&PrecipitationsLgm[iv*12+month],gauss); + dinput1->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month); + dinput2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss,month); + dinput3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month); + dinput4->GetInputValue(&PrecipitationsLgm[iv*12+month],gauss,month); PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 24582) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 24583) @@ -228,5 +228,26 @@ void Penta::DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,Inputs2* inputs2,IoModel* iomodel,int input_enum){/*{{{*/ - _error_("Look at Tria :)"); + /*Intermediaries*/ + int vertexsids[NUMVERTICES]; + int vertexlids[NUMVERTICES]; + IssmDouble nodeinputs[NUMVERTICES]; + + /*Some sanity checks*/ + if(num_inputs<1) _error_("Cannot create a DatasetInput of size <1"); + if(M!=iomodel->numberofvertices) _error_("Input size not supported yet"); + if(N!=num_inputs) _error_("Sizes are not consistent"); + + /*Get indices*/ + _assert_(iomodel->elements); + for(int i=0;i(iomodel->elements[NUMVERTICES*this->Sid()+i])-1; + vertexlids[i] = iomodel->my_vertices_lids[vertexsids[i]]; + } + + /*Create inputs and add to DataSetInput*/ + for(int i=0;iSetPentaDatasetInput(input_enum,individual_enums[i],P1Enum,NUMVERTICES,vertexlids,nodeinputs); + } } /*}}}*/ Index: /issm/trunk-jpl/src/m/classes/SMBd18opdd.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBd18opdd.m (revision 24582) +++ /issm/trunk-jpl/src/m/classes/SMBd18opdd.m (revision 24583) @@ -93,16 +93,6 @@ if(self.isd18opd==1) - md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1); - %cross check that time stamp values are between 0 and 1: - i=find(self.temperatures_presentday(end,:)-self.precipitations_presentday(end,:)~=0); - if length(i), - error('Smb checkconsistency fail: timestamp (last row of matrix) for present day temp and precip must be equal!'); - end - j=find(self.temperatures_presentday(end,:)<0 | self.temperatures_presentday(end,:)>1); - if length(j), - error('Smb checkconsistency fail: timestamp (last row of matrix) for present day temp and precip must be between 0 and 1!'); - end - + md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices 12],'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices 12],'NaN',1,'Inf',1,'timeseries',1); md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,NaN],'singletimeseries',1); md = checkfield(md,'fieldname','smb.dpermil','>=',0,'numel',1); Index: /issm/trunk-jpl/src/m/classes/SMBd18opdd.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBd18opdd.py (revision 24582) +++ /issm/trunk-jpl/src/m/classes/SMBd18opdd.py (revision 24583) @@ -141,6 +141,6 @@ multt = np.ceil(lent / 12.) * 12. multp = np.ceil(lenp / 12.) * 12. - md = checkfield(md, 'fieldname', 'smb.temperatures_presentday', 'size', [md.mesh.numberofvertices + 1, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'smb.precipitations_presentday', 'size', [md.mesh.numberofvertices + 1, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'smb.temperatures_presentday', 'size', [md.mesh.numberofvertices, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'smb.precipitations_presentday', 'size', [md.mesh.numberofvertices, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) if self.istemperaturescaled == 0: Index: /issm/trunk-jpl/src/m/classes/SMBpdd.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBpdd.m (revision 24582) +++ /issm/trunk-jpl/src/m/classes/SMBpdd.m (revision 24583) @@ -90,15 +90,15 @@ md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,NaN],'singletimeseries',1); md = checkfield(md,'fieldname','smb.delta18o_surface','NaN',1,'Inf',1,'size',[2,NaN],'singletimeseries',1); - md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','smb.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','smb.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices 12],'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','smb.temperatures_lgm','size',[md.mesh.numberofvertices 12],'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices 12],'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','smb.precipitations_lgm','size',[md.mesh.numberofvertices 12],'NaN',1,'Inf',1,'timeseries',1); md = checkfield(md,'fieldname','smb.Tdiff','NaN',1,'Inf',1,'size',[2,NaN],'singletimeseries',1); md = checkfield(md,'fieldname','smb.sealev','NaN',1,'Inf',1,'size',[2,NaN],'singletimeseries',1); elseif(self.ismungsm==1) - md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','smb.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','smb.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices 12],'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','smb.temperatures_lgm','size',[md.mesh.numberofvertices 12],'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices 12],'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','smb.precipitations_lgm','size',[md.mesh.numberofvertices 12],'NaN',1,'Inf',1,'timeseries',1); md = checkfield(md,'fieldname','smb.Pfac','NaN',1,'Inf',1,'size',[2,NaN],'singletimeseries',1); md = checkfield(md,'fieldname','smb.Tdiff','NaN',1,'Inf',1,'size',[2,NaN],'singletimeseries',1); Index: /issm/trunk-jpl/src/m/classes/SMBpdd.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBpdd.py (revision 24582) +++ /issm/trunk-jpl/src/m/classes/SMBpdd.py (revision 24583) @@ -143,18 +143,18 @@ md = checkfield(md, 'fieldname', 'smb.delta18o', 'NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1) md = checkfield(md, 'fieldname', 'smb.delta18o_surface', 'NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1) - md = checkfield(md, 'fieldname', 'smb.temperatures_presentday', 'size', [md.mesh.numberofvertices + 1, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'smb.temperatures_lgm', 'size', [md.mesh.numberofvertices + 1, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'smb.precipitations_presentday', 'size', [md.mesh.numberofvertices + 1, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'smb.precipitations_lgm', 'size', [md.mesh.numberofvertices + 1, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'smb.temperatures_presentday', 'size', [md.mesh.numberofvertices, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'smb.temperatures_lgm', 'size', [md.mesh.numberofvertices, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'smb.precipitations_presentday', 'size', [md.mesh.numberofvertices, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'smb.precipitations_lgm', 'size', [md.mesh.numberofvertices, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) md = checkfield(md, 'fieldname', 'smb.Tdiff', 'NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1) md = checkfield(md, 'fieldname', 'smb.sealev', 'NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1) elif self.ismungsm: - md = checkfield(md, 'fieldname', 'smb.temperatures_presentday', 'size', [md.mesh.numberofvertices + 1, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'smb.temperatures_lgm', 'size', [md.mesh.numberofvertices + 1, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'smb.precipitations_presentday', 'size', [md.mesh.numberofvertices + 1, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) - md = checkfield(md, 'fieldname', 'smb.precipitations_lgm', 'size', [md.mesh.numberofvertices + 1, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'smb.temperatures_presentday', 'size', [md.mesh.numberofvertices, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'smb.temperatures_lgm', 'size', [md.mesh.numberofvertices, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'smb.precipitations_presentday', 'size', [md.mesh.numberofvertices, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'smb.precipitations_lgm', 'size', [md.mesh.numberofvertices, 12], 'NaN', 1, 'Inf', 1, 'timeseries', 1) md = checkfield(md, 'fieldname', 'smb.Pfac', 'NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1) - md = checkfield(md, 'fieldname', 'smb.Tdiff', 'NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1) - md = checkfield(md, 'fieldname', 'smb.sealev', 'NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1) + md = checkfield(md, 'fieldname', 'smb.Tdiff','NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1) + md = checkfield(md, 'fieldname', 'smb.sealev','NaN',1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1) md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])