Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27317) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27318) @@ -246,13 +246,14 @@ case FrontalForcingsRignotarmaEnum: /*Retrieve autoregressive parameters*/ + parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_params",FrontalForcingsNumberofParamsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_breaks",FrontalForcingsNumberofBreaksEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_order",FrontalForcingsARMAarOrderEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ma_order",FrontalForcingsARMAmaOrderEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.arma_initialtime",FrontalForcingsARMAInitialTimeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.arma_timestep",FrontalForcingsARMATimestepEnum)); - iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.const"); - parameters->AddObject(new DoubleVecParam(FrontalForcingsARMAconstEnum,transparam,N)); + iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.polynomialparams"); + parameters->AddObject(new DoubleMatParam(FrontalForcingsARMApolyparamsEnum,transparam,M,N)); xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.trend"); - parameters->AddObject(new DoubleVecParam(FrontalForcingsARMAtrendEnum,transparam,N)); + iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.datebreaks"); + parameters->AddObject(new DoubleMatParam(FrontalForcingsARMAdatebreaksEnum,transparam,M,N)); xDelete(transparam); iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.arlag_coefs"); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27317) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27318) @@ -61,11 +61,15 @@ return false; }/*}}}*/ -void Element::ArmaProcess(bool isstepforarma,int arorder,int maorder,IssmDouble telapsed_arma,IssmDouble tstep_arma,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,bool isfieldstochastic,int enum_type){/*{{{*/ +void Element::ArmaProcess(bool isstepforarma,int arorder,int maorder,int numparams,int numbreaks,IssmDouble tstep_arma,IssmDouble* polyparams,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,IssmDouble* datebreaks,bool isfieldstochastic,int enum_type){/*{{{*/ const int numvertices = this->GetNumberOfVertices(); - int basinid,M,N,arenum_type,maenum_type,basinenum_type,noiseenum_type,outenum_type; - IssmDouble time,dt,starttime,termconstant_basin,trend_basin,noiseterm; + int numperiods = numbreaks+1; + int basinid,M,N,arenum_type,maenum_type,basinenum_type,noiseenum_type,outenum_type,indperiod; + IssmDouble time,dt,starttime,noiseterm; IssmDouble* arlagcoefs_basin = xNew(arorder); IssmDouble* malagcoefs_basin = xNew(maorder); + IssmDouble* datebreaks_basin = xNew(numbreaks); + IssmDouble* polyparams_basin = xNew(numperiods*numparams); IssmDouble* varlist = xNew(numvertices); + IssmDouble* sumpoly = xNewZeroInit(arorder+1); IssmDouble* valuesautoregression = NULL; IssmDouble* valuesmovingaverage = NULL; @@ -104,14 +108,36 @@ /*Get basin coefficients*/ this->GetInputValue(&basinid,basinenum_type); - for(int ii=0;ii0){ + for(int i=0;i0){ + /*Find index of tatstep compared to the breakpoints*/ + indperiod = 0; + for(int i=0;idatebreaks_basin[i]) indperiod = i+1; + } + /*Compute polynomial with parameters of indperiod*/ + if(indperiod==0) telapsed_break = tatstep; + else telapsed_break = tatstep-datebreaks_basin[indperiod]; + for(int j=0;j(numvertices*arorder); IssmDouble* initvaluesmovingaverage = xNewZeroInit(numvertices*maorder); - for(int i=0;iinputs->SetArrayInput(arenum_type,this->lid,initvaluesautoregression,numvertices*arorder); this->inputs->SetArrayInput(maenum_type,this->lid,initvaluesmovingaverage,numvertices*maorder); @@ -141,11 +167,11 @@ /*Compute autoregressive term*/ IssmDouble autoregressionterm=0.; - for(int ii=0;ii(arlagcoefs_basin); xDelete(malagcoefs_basin); + xDelete(datebreaks_basin); + xDelete(polyparams_basin); + xDelete(sumpoly); xDelete(varlist); xDelete(valuesautoregression); xDelete(valuesmovingaverage); -}/*}}}*/ -void Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble tstep_ar,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* lagcoefs,bool isfieldstochastic,int enum_type){/*{{{*/ - - const int numvertices = this->GetNumberOfVertices(); - int basinid,M,N,arenum_type,basinenum_type,noiseenum_type,outenum_type; - IssmDouble time,dt,starttime,termconstant_basin,trend_basin,noiseterm; - IssmDouble* lagcoefs_basin = xNew(arorder); - IssmDouble* varlist = xNew(numvertices); - IssmDouble* valuesautoregression = NULL; - Input* noiseterm_input = NULL; - - /*Get field-specific enums*/ - switch(enum_type){ - case(SMBarmaEnum): - arenum_type = SmbValuesAutoregressionEnum; - basinenum_type = SmbBasinsIdEnum; - noiseenum_type = SmbARMANoiseEnum; - outenum_type = SmbMassBalanceEnum; - break; - case(FrontalForcingsRignotarmaEnum): - arenum_type = ThermalforcingValuesAutoregressionEnum; - basinenum_type = FrontalForcingsBasinIdEnum; - noiseenum_type = ThermalforcingARMANoiseEnum; - outenum_type = ThermalForcingEnum; - break; - case(BasalforcingsDeepwaterMeltingRatearmaEnum): - arenum_type = BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum; - basinenum_type = BasalforcingsLinearBasinIdEnum; - noiseenum_type = BasalforcingsDeepwaterMeltingRateNoiseEnum; - outenum_type = BasalforcingsSpatialDeepwaterMeltingRateEnum; - break; - } - - /*Get time parameters*/ - this->parameters->FindParam(&time,TimeEnum); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); - - /*Get basin coefficients*/ - this->GetInputValue(&basinid,basinenum_type); - for(int ii=0;ii(numvertices*arorder); - for(int i=0;iinputs->SetArrayInput(arenum_type,this->lid,initvaluesautoregression,numvertices*arorder); - xDelete(initvaluesautoregression); - } - - /*Get noise and autoregressive terms*/ - if(isfieldstochastic){ - noiseterm_input = this->GetInput(noiseenum_type); - Gauss* gauss = this->NewGauss(); - noiseterm_input->GetInputValue(&noiseterm,gauss); - delete gauss; - } - else noiseterm = 0.0; - this->inputs->GetArray(arenum_type,this->lid,&valuesautoregression,&M); - - /*If not AR model timestep: take the old values of variable*/ - if(isstepforar==false){ - for(int i=0;i(numvertices*arorder); - /*Assign newest values and shift older values*/ - for(int i=0;iinputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder); - xDelete(temparray); - } - - /*Add input to element*/ - this->AddInput(outenum_type,varlist,P1Enum); - - /*Cleanup*/ - xDelete(lagcoefs_basin); - xDelete(varlist); - xDelete(valuesautoregression); -}/*}}}*/ +} + +/*}}}*/ void Element::BasinLinearFloatingiceMeltingRate(IssmDouble* deepwaterel,IssmDouble* upperwatermelt,IssmDouble* upperwaterel,IssmDouble* perturbation){/*{{{*/ @@ -3889,5 +3827,5 @@ /* calculate melt rates */ meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s] - } + } if(xIsNan(meltrates[iv])) _error_("NaN found in vector"); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27317) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27318) @@ -68,6 +68,6 @@ /*bool AnyActive(void);*/ bool AnyFSet(void); - void ArmaProcess(bool isstepforarma,int arorder,int maorder,IssmDouble telapsed_arma,IssmDouble tstep_arma,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,bool isfieldstochastic,int enum_type); - void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble tstep_ar,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* lagcoefs,bool isfieldstochastic,int enum_type); + void ArmaProcess_pre18Oct2022(bool isstepforarma,int arorder,int maorder,IssmDouble telapsed_arma,IssmDouble tstep_arma,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,bool isfieldstochastic,int enum_type); + void ArmaProcess(bool isstepforarma,int arorder,int maorder,int numparams,int numbreaks,IssmDouble tstep_arma,IssmDouble* polyparams,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,IssmDouble* datebreaks,bool isfieldstochastic,int enum_type); void BasinLinearFloatingiceMeltingRate(IssmDouble* deepwaterel,IssmDouble* upperwatermelt,IssmDouble* upperwaterel,IssmDouble* perturbation); void CalvingSetZeroRate(void); Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 27317) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 27318) @@ -231,14 +231,15 @@ bool isstochastic; bool isdeepmeltingstochastic = false; - int M,N,Narlagcoefs,Nmalagcoefs,arorder,maorder,numbasins,my_rank; + int M,N,arorder,maorder,numbasins,numparams,numbreaks,my_rank; femmodel->parameters->FindParam(&numbasins,BasalforcingsLinearNumBasinsEnum); femmodel->parameters->FindParam(&arorder,BasalforcingsARMAarOrderEnum); femmodel->parameters->FindParam(&maorder,BasalforcingsARMAmaOrderEnum); - IssmDouble tinit_arma; - IssmDouble* termconstant = NULL; - IssmDouble* trend = NULL; - IssmDouble* arlagcoefs = NULL; + femmodel->parameters->FindParam(&numparams,BasalforcingsLinearNumParamsEnum); + femmodel->parameters->FindParam(&numbreaks,BasalforcingsLinearNumBreaksEnum); + IssmDouble* datebreaks = NULL; + IssmDouble* arlagcoefs = NULL; IssmDouble* malagcoefs = NULL; - IssmDouble* deepwaterel = NULL; + IssmDouble* polyparams = NULL; + IssmDouble* deepwaterel = NULL; IssmDouble* upperwaterel = NULL; IssmDouble* upperwatermelt = NULL; @@ -246,9 +247,8 @@ /*Get autoregressive parameters*/ - femmodel->parameters->FindParam(&tinit_arma,BasalforcingsARMAInitialTimeEnum); - femmodel->parameters->FindParam(&termconstant,&M,BasalforcingsARMAconstEnum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&trend,&M,BasalforcingsARMAtrendEnum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&arlagcoefs,&M,&Narlagcoefs,BasalforcingsARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(Narlagcoefs==arorder); - femmodel->parameters->FindParam(&malagcoefs,&M,&Nmalagcoefs,BasalforcingsARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(Nmalagcoefs==maorder); + femmodel->parameters->FindParam(&datebreaks,&M,&N,BasalforcingsARMAdatebreaksEnum); _assert_(M==numbasins); _assert_(N==numbreaks); + femmodel->parameters->FindParam(&polyparams,&M,&N,BasalforcingsARMApolyparamsEnum); _assert_(M==numbasins); _assert_(N==numbreaks*numparams); + femmodel->parameters->FindParam(&arlagcoefs,&M,&N,BasalforcingsARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(N==arorder); + femmodel->parameters->FindParam(&malagcoefs,&M,&N,BasalforcingsARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(N==maorder); /*Get basin-specific parameters*/ @@ -269,6 +269,4 @@ xDelete(stochasticfields); } - /*Time elapsed with respect to ARMA model initial time*/ - IssmDouble telapsed_arma = time-tinit_arma; /*Loop over each element to compute FloatingiceMeltingRate at vertices*/ @@ -276,13 +274,13 @@ Element* element = xDynamicCast(object); /*Compute ARMA*/ - element->ArmaProcess(isstepforarma,arorder,maorder,telapsed_arma,tstep_arma,termconstant,trend,arlagcoefs,malagcoefs,isdeepmeltingstochastic,BasalforcingsDeepwaterMeltingRatearmaEnum); + element->ArmaProcess(isstepforarma,arorder,maorder,numparams,numbreaks,tstep_arma,polyparams,arlagcoefs,malagcoefs,datebreaks,isdeepmeltingstochastic,BasalforcingsDeepwaterMeltingRatearmaEnum); element->BasinLinearFloatingiceMeltingRate(deepwaterel,upperwatermelt,upperwaterel,perturbation); } /*Cleanup*/ - xDelete(termconstant); - xDelete(trend); xDelete(arlagcoefs); xDelete(malagcoefs); + xDelete(polyparams); + xDelete(datebreaks); xDelete(deepwaterel); xDelete(upperwaterel); Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 27317) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 27318) @@ -50,21 +50,21 @@ bool isstochastic; bool istfstochastic = false; - int M,N,Narlagcoefs,Nmalagcoefs,arorder,maorder,numbasins,my_rank; + int M,N,arorder,maorder,numbasins,numparams,numbreaks,my_rank; femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); + femmodel->parameters->FindParam(&numparams,FrontalForcingsNumberofParamsEnum); + femmodel->parameters->FindParam(&numbreaks,FrontalForcingsNumberofBreaksEnum); femmodel->parameters->FindParam(&arorder,FrontalForcingsARMAarOrderEnum); femmodel->parameters->FindParam(&maorder,FrontalForcingsARMAmaOrderEnum); - IssmDouble tinit_arma; - IssmDouble* termconstant = NULL; - IssmDouble* trend = NULL; + IssmDouble* datebreaks = NULL; IssmDouble* arlagcoefs = NULL; IssmDouble* malagcoefs = NULL; IssmDouble* monthlyeff = NULL; + IssmDouble* polyparams = NULL; - femmodel->parameters->FindParam(&tinit_arma,FrontalForcingsARMAInitialTimeEnum); - femmodel->parameters->FindParam(&termconstant,&M,FrontalForcingsARMAconstEnum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&trend,&M,FrontalForcingsARMAtrendEnum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&arlagcoefs,&M,&Narlagcoefs,FrontalForcingsARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(Narlagcoefs==arorder); - femmodel->parameters->FindParam(&malagcoefs,&M,&Nmalagcoefs,FrontalForcingsARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(Nmalagcoefs==maorder); - femmodel->parameters->FindParam(&monthlyeff,&M,&N,ThermalForcingMonthlyEffectsEnum); _assert_(M==numbasins); _assert_(N==12); + femmodel->parameters->FindParam(&datebreaks,&M,&N,FrontalForcingsARMAdatebreaksEnum); _assert_(M==numbasins); _assert_(N==numbreaks); + femmodel->parameters->FindParam(&polyparams,&M,&N,FrontalForcingsARMApolyparamsEnum); _assert_(M==numbasins); _assert_(N==numbreaks*numparams); + femmodel->parameters->FindParam(&arlagcoefs,&M,&N,FrontalForcingsARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(N==arorder); + femmodel->parameters->FindParam(&malagcoefs,&M,&N,FrontalForcingsARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(N==maorder); + femmodel->parameters->FindParam(&monthlyeff,&M,&N,ThermalForcingMonthlyEffectsEnum); _assert_(M==numbasins); _assert_(N==12); femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); @@ -79,6 +79,4 @@ xDelete(stochasticfields); } - /*Time elapsed with respect to ARMA model initial time*/ - IssmDouble telapsed_arma = time-tinit_arma; /*Loop over each element to compute Thermal Forcing at vertices*/ @@ -86,5 +84,5 @@ Element* element = xDynamicCast(object); /*Compute ARMA*/ - element->ArmaProcess(isstepforarma,arorder,maorder,telapsed_arma,tstep_arma,termconstant,trend,arlagcoefs,malagcoefs,istfstochastic,FrontalForcingsRignotarmaEnum); + element->ArmaProcess(isstepforarma,arorder,maorder,numparams,numbreaks,tstep_arma,polyparams,arlagcoefs,malagcoefs,datebreaks,istfstochastic,FrontalForcingsRignotarmaEnum); /*Compute monthly effects*/ element->MonthlyEffectBasin(monthlyeff,FrontalForcingsRignotarmaEnum); @@ -92,8 +90,8 @@ /*Cleanup*/ - xDelete(termconstant); - xDelete(trend); xDelete(arlagcoefs); xDelete(malagcoefs); xDelete(monthlyeff); + xDelete(polyparams); + xDelete(datebreaks); }/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27317) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27318) @@ -259,13 +259,14 @@ /*Add parameters that are not in standard nbvertices format*/ parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsLinearNumBasinsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_breaks",BasalforcingsLinearNumBreaksEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_params",BasalforcingsLinearNumParamsEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.ar_order",BasalforcingsARMAarOrderEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.ma_order",BasalforcingsARMAmaOrderEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.arma_initialtime",BasalforcingsARMAInitialTimeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.arma_timestep",BasalforcingsARMATimestepEnum)); - iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.const"); - parameters->AddObject(new DoubleVecParam(BasalforcingsARMAconstEnum,transparam,N)); - xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.trend"); - parameters->AddObject(new DoubleVecParam(BasalforcingsARMAtrendEnum,transparam,N)); + iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.datebreaks"); + parameters->AddObject(new DoubleMatParam(BasalforcingsARMAdatebreaksEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.polynomialparams"); + parameters->AddObject(new DoubleMatParam(BasalforcingsARMApolyparamsEnum,transparam,M,N)); xDelete(transparam); iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.arlag_coefs"); @@ -435,16 +436,17 @@ /*Add parameters that are not in standard nbvertices format*/ parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_basins",SmbNumBasinsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_params",SmbNumParamsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_breaks",SmbNumBreaksEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_order",SmbARMAarOrderEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.ma_order",SmbARMAmaOrderEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.smb.arma_initialtime",SmbARMAInitialTimeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.arma_timestep",SmbARMATimestepEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_bins",SmbNumElevationBinsEnum)); - iomodel->FetchData(&transparam,&M,&N,"md.smb.const"); - parameters->AddObject(new DoubleVecParam(SmbARMAconstEnum,transparam,N)); - xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.smb.trend"); - parameters->AddObject(new DoubleVecParam(SmbARMAtrendEnum,transparam,N)); - xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.smb.arlag_coefs"); + iomodel->FetchData(&transparam,&M,&N,"md.smb.datebreaks"); + parameters->AddObject(new DoubleMatParam(SmbARMAdatebreaksEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.smb.polynomialparams"); + parameters->AddObject(new DoubleMatParam(SmbARMApolyparamsEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.smb.arlag_coefs"); parameters->AddObject(new DoubleMatParam(SmbARMAarlagcoefsEnum,transparam,M,N)); xDelete(transparam); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 27317) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 27318) @@ -168,23 +168,23 @@ bool isstochastic; bool issmbstochastic = false; - int M,N,Narlagcoefs,Nmalagcoefs,arorder,maorder,numbasins,numelevbins,my_rank; + int M,N,arorder,maorder,numbasins,numparams,numbreaks,numelevbins,my_rank; femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum); - femmodel->parameters->FindParam(&arorder,SmbARMAarOrderEnum); + femmodel->parameters->FindParam(&numparams,SmbNumParamsEnum); + femmodel->parameters->FindParam(&numbreaks,SmbNumBreaksEnum); + femmodel->parameters->FindParam(&arorder,SmbARMAarOrderEnum); femmodel->parameters->FindParam(&maorder,SmbARMAmaOrderEnum); femmodel->parameters->FindParam(&numelevbins,SmbNumElevationBinsEnum); - IssmDouble tinit_arma; - IssmDouble* termconstant = NULL; - IssmDouble* trend = NULL; - IssmDouble* arlagcoefs = NULL; + IssmDouble* datebreaks = NULL; + IssmDouble* arlagcoefs = NULL; IssmDouble* malagcoefs = NULL; + IssmDouble* polyparams = NULL; IssmDouble* lapserates = NULL; IssmDouble* elevbins = NULL; IssmDouble* refelevation = NULL; - femmodel->parameters->FindParam(&tinit_arma,SmbARMAInitialTimeEnum); - femmodel->parameters->FindParam(&termconstant,&M,SmbARMAconstEnum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&trend,&M,SmbARMAtrendEnum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&arlagcoefs,&M,&Narlagcoefs,SmbARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(Narlagcoefs==arorder); - femmodel->parameters->FindParam(&malagcoefs,&M,&Nmalagcoefs,SmbARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(Nmalagcoefs==maorder); + femmodel->parameters->FindParam(&datebreaks,&M,&N,SmbARMAdatebreaksEnum); _assert_(M==numbasins); _assert_(N==numbreaks); + femmodel->parameters->FindParam(&polyparams,&M,&N,SmbARMApolyparamsEnum); _assert_(M==numbasins); _assert_(N==numbreaks*numparams); + femmodel->parameters->FindParam(&arlagcoefs,&M,&N,SmbARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(N==arorder); + femmodel->parameters->FindParam(&malagcoefs,&M,&N,SmbARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(N==maorder); femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum); _assert_(M==numbasins); _assert_(N==numelevbins); femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum); _assert_(M==numbasins); _assert_(N==numelevbins-1); @@ -202,6 +202,4 @@ xDelete(stochasticfields); } - /*Time elapsed with respect to ARMA model initial time*/ - IssmDouble telapsed_arma = time-tinit_arma; /*Loop over each element to compute SMB at vertices*/ @@ -209,5 +207,5 @@ Element* element = xDynamicCast(object); /*Compute ARMA*/ - element->ArmaProcess(isstepforarma,arorder,maorder,telapsed_arma,tstep_arma,termconstant,trend,arlagcoefs,malagcoefs,issmbstochastic,SMBarmaEnum); + element->ArmaProcess(isstepforarma,arorder,maorder,numparams,numbreaks,tstep_arma,polyparams,arlagcoefs,malagcoefs,datebreaks,issmbstochastic,SMBarmaEnum); /*Compute lapse rate adjustment*/ element->LapseRateBasinSMB(numelevbins,lapserates,elevbins,refelevation); @@ -215,8 +213,8 @@ /*Cleanup*/ - xDelete(termconstant); - xDelete(trend); xDelete(arlagcoefs); xDelete(malagcoefs); + xDelete(polyparams); + xDelete(datebreaks); xDelete(lapserates); xDelete(elevbins); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27317) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27318) @@ -67,10 +67,7 @@ syn keyword cConstant BalancethicknessStabilizationEnum syn keyword cConstant BarystaticContributionsEnum -syn keyword cConstant BasalforcingsARMAInitialTimeEnum syn keyword cConstant BasalforcingsARMATimestepEnum syn keyword cConstant BasalforcingsARMAarOrderEnum syn keyword cConstant BasalforcingsARMAmaOrderEnum -syn keyword cConstant BasalforcingsARMAconstEnum -syn keyword cConstant BasalforcingsARMAtrendEnum syn keyword cConstant BasalforcingsBottomplumedepthEnum syn keyword cConstant BasalforcingsCrustthicknessEnum @@ -87,10 +84,14 @@ syn keyword cConstant BasalforcingsIsmip6TfDepthsEnum syn keyword cConstant BasalforcingsLinearNumBasinsEnum -syn keyword cConstant BasalforcingsLowercrustheatEnum +syn keyword cConstant BasalforcingsLinearNumBreaksEnum +syn keyword cConstant BasalforcingsLinearNumParamsEnum syn keyword cConstant BasalforcingsMantleconductivityEnum syn keyword cConstant BasalforcingsNusseltEnum syn keyword cConstant BasalforcingsARMAarlagcoefsEnum +syn keyword cConstant BasalforcingsARMAdatebreaksEnum syn keyword cConstant BasalforcingsARMAmalagcoefsEnum +syn keyword cConstant BasalforcingsARMApolyparamsEnum syn keyword cConstant BasalforcingsIsThermalForcingEnum +syn keyword cConstant BasalforcingsLowercrustheatEnum syn keyword cConstant BasalforcingsPicoAverageOverturningEnum syn keyword cConstant BasalforcingsPicoAverageSalinityEnum @@ -213,11 +214,12 @@ syn keyword cConstant FrictionVoidRatioEnum syn keyword cConstant FrontalForcingsBasinIcefrontAreaEnum -syn keyword cConstant FrontalForcingsARMAInitialTimeEnum syn keyword cConstant FrontalForcingsARMATimestepEnum syn keyword cConstant FrontalForcingsARMAarOrderEnum syn keyword cConstant FrontalForcingsARMAmaOrderEnum -syn keyword cConstant FrontalForcingsARMAconstEnum -syn keyword cConstant FrontalForcingsARMAtrendEnum +syn keyword cConstant FrontalForcingsARMAdatebreaksEnum +syn keyword cConstant FrontalForcingsARMApolyparamsEnum syn keyword cConstant FrontalForcingsNumberofBasinsEnum +syn keyword cConstant FrontalForcingsNumberofBreaksEnum +syn keyword cConstant FrontalForcingsNumberofParamsEnum syn keyword cConstant FrontalForcingsParamEnum syn keyword cConstant FrontalForcingsARMAarlagcoefsEnum @@ -318,4 +320,5 @@ syn keyword cConstant LoveMinIntegrationStepsEnum syn keyword cConstant LoveMaxIntegrationdrEnum +syn keyword cConstant LoveIntegrationSchemeEnum syn keyword cConstant LoveKernelsEnum syn keyword cConstant LoveMu0Enum @@ -493,11 +496,8 @@ syn keyword cConstant SmbAccurefEnum syn keyword cConstant SmbAdThreshEnum -syn keyword cConstant SmbARMAInitialTimeEnum syn keyword cConstant SmbARMATimestepEnum syn keyword cConstant SmbARMAarOrderEnum syn keyword cConstant SmbARMAmaOrderEnum syn keyword cConstant SmbAveragingEnum -syn keyword cConstant SmbARMAconstEnum -syn keyword cConstant SmbARMAtrendEnum syn keyword cConstant SmbDesfacEnum syn keyword cConstant SmbDpermilEnum @@ -533,9 +533,13 @@ syn keyword cConstant SmbLapseRatesEnum syn keyword cConstant SmbNumBasinsEnum +syn keyword cConstant SmbNumBreaksEnum syn keyword cConstant SmbNumElevationBinsEnum +syn keyword cConstant SmbNumParamsEnum syn keyword cConstant SmbNumRequestedOutputsEnum syn keyword cConstant SmbPfacEnum syn keyword cConstant SmbARMAarlagcoefsEnum +syn keyword cConstant SmbARMAdatebreaksEnum syn keyword cConstant SmbARMAmalagcoefsEnum +syn keyword cConstant SmbARMApolyparamsEnum syn keyword cConstant SmbRdlEnum syn keyword cConstant SmbRefElevationEnum @@ -939,4 +943,5 @@ syn keyword cConstant SealevelUNorthEsaEnum syn keyword cConstant SealevelchangeIndicesEnum +syn keyword cConstant SealevelchangeConvolutionVerticesEnum syn keyword cConstant SealevelchangeAlphaIndexEnum syn keyword cConstant SealevelchangeAzimuthIndexEnum @@ -957,4 +962,5 @@ syn keyword cConstant SealevelchangeViscousNEnum syn keyword cConstant SealevelchangeViscousEEnum +syn keyword cConstant CouplingTransferCountEnum syn keyword cConstant SedimentHeadEnum syn keyword cConstant SedimentHeadOldEnum @@ -1680,4 +1686,5 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel +syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1686,6 +1693,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DatasetInput syn keyword cType DataSetParam -syn keyword cType DatasetInput syn keyword cType Definition syn keyword cType DependentObject @@ -1700,6 +1707,6 @@ syn keyword cType ElementInput syn keyword cType ElementMatrix +syn keyword cType Elements syn keyword cType ElementVector -syn keyword cType Elements syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1708,9 +1715,10 @@ syn keyword cType Friction syn keyword cType Gauss +syn keyword cType GaussianVariogram +syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria -syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1729,4 +1737,5 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface +syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1739,4 +1748,5 @@ syn keyword cType Matice syn keyword cType Matlitho +syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1751,6 +1761,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType Options syn keyword cType OptionUtilities -syn keyword cType Options syn keyword cType Param syn keyword cType Parameters @@ -1766,11 +1776,11 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType Riftfront syn keyword cType RiftStruct -syn keyword cType Riftfront syn keyword cType SealevelGeometry syn keyword cType Seg syn keyword cType SegInput +syn keyword cType Segment syn keyword cType SegRef -syn keyword cType Segment syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1791,8 +1801,4 @@ syn keyword cType Vertex syn keyword cType Vertices -syn keyword cType classes -syn keyword cType gaussobjects -syn keyword cType krigingobjects -syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis @@ -1805,4 +1811,5 @@ syn keyword cType BalancevelocityAnalysis syn keyword cType DamageEvolutionAnalysis +syn keyword cType DebrisAnalysis syn keyword cType DepthAverageAnalysis syn keyword cType EnthalpyAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27317) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27318) @@ -61,10 +61,7 @@ BalancethicknessStabilizationEnum, BarystaticContributionsEnum, - BasalforcingsARMAInitialTimeEnum, BasalforcingsARMATimestepEnum, BasalforcingsARMAarOrderEnum, BasalforcingsARMAmaOrderEnum, - BasalforcingsARMAconstEnum, - BasalforcingsARMAtrendEnum, BasalforcingsBottomplumedepthEnum, BasalforcingsCrustthicknessEnum, @@ -81,10 +78,14 @@ BasalforcingsIsmip6TfDepthsEnum, BasalforcingsLinearNumBasinsEnum, - BasalforcingsLowercrustheatEnum, + BasalforcingsLinearNumBreaksEnum, + BasalforcingsLinearNumParamsEnum, BasalforcingsMantleconductivityEnum, BasalforcingsNusseltEnum, BasalforcingsARMAarlagcoefsEnum, + BasalforcingsARMAdatebreaksEnum, BasalforcingsARMAmalagcoefsEnum, + BasalforcingsARMApolyparamsEnum, BasalforcingsIsThermalForcingEnum, + BasalforcingsLowercrustheatEnum, BasalforcingsPicoAverageOverturningEnum, BasalforcingsPicoAverageSalinityEnum, @@ -207,11 +208,12 @@ FrictionVoidRatioEnum, FrontalForcingsBasinIcefrontAreaEnum, - FrontalForcingsARMAInitialTimeEnum, FrontalForcingsARMATimestepEnum, FrontalForcingsARMAarOrderEnum, FrontalForcingsARMAmaOrderEnum, - FrontalForcingsARMAconstEnum, - FrontalForcingsARMAtrendEnum, + FrontalForcingsARMAdatebreaksEnum, + FrontalForcingsARMApolyparamsEnum, FrontalForcingsNumberofBasinsEnum, + FrontalForcingsNumberofBreaksEnum, + FrontalForcingsNumberofParamsEnum, FrontalForcingsParamEnum, FrontalForcingsARMAarlagcoefsEnum, @@ -488,11 +490,8 @@ SmbAccurefEnum, SmbAdThreshEnum, - SmbARMAInitialTimeEnum, SmbARMATimestepEnum, SmbARMAarOrderEnum, SmbARMAmaOrderEnum, SmbAveragingEnum, - SmbARMAconstEnum, - SmbARMAtrendEnum, SmbDesfacEnum, SmbDpermilEnum, @@ -528,9 +527,13 @@ SmbLapseRatesEnum, SmbNumBasinsEnum, + SmbNumBreaksEnum, SmbNumElevationBinsEnum, + SmbNumParamsEnum, SmbNumRequestedOutputsEnum, SmbPfacEnum, SmbARMAarlagcoefsEnum, + SmbARMAdatebreaksEnum, SmbARMAmalagcoefsEnum, + SmbARMApolyparamsEnum, SmbRdlEnum, SmbRefElevationEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27317) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27318) @@ -69,10 +69,7 @@ case BalancethicknessStabilizationEnum : return "BalancethicknessStabilization"; case BarystaticContributionsEnum : return "BarystaticContributions"; - case BasalforcingsARMAInitialTimeEnum : return "BasalforcingsARMAInitialTime"; case BasalforcingsARMATimestepEnum : return "BasalforcingsARMATimestep"; case BasalforcingsARMAarOrderEnum : return "BasalforcingsARMAarOrder"; case BasalforcingsARMAmaOrderEnum : return "BasalforcingsARMAmaOrder"; - case BasalforcingsARMAconstEnum : return "BasalforcingsARMAconst"; - case BasalforcingsARMAtrendEnum : return "BasalforcingsARMAtrend"; case BasalforcingsBottomplumedepthEnum : return "BasalforcingsBottomplumedepth"; case BasalforcingsCrustthicknessEnum : return "BasalforcingsCrustthickness"; @@ -89,10 +86,14 @@ case BasalforcingsIsmip6TfDepthsEnum : return "BasalforcingsIsmip6TfDepths"; case BasalforcingsLinearNumBasinsEnum : return "BasalforcingsLinearNumBasins"; - case BasalforcingsLowercrustheatEnum : return "BasalforcingsLowercrustheat"; + case BasalforcingsLinearNumBreaksEnum : return "BasalforcingsLinearNumBreaks"; + case BasalforcingsLinearNumParamsEnum : return "BasalforcingsLinearNumParams"; case BasalforcingsMantleconductivityEnum : return "BasalforcingsMantleconductivity"; case BasalforcingsNusseltEnum : return "BasalforcingsNusselt"; case BasalforcingsARMAarlagcoefsEnum : return "BasalforcingsARMAarlagcoefs"; + case BasalforcingsARMAdatebreaksEnum : return "BasalforcingsARMAdatebreaks"; case BasalforcingsARMAmalagcoefsEnum : return "BasalforcingsARMAmalagcoefs"; + case BasalforcingsARMApolyparamsEnum : return "BasalforcingsARMApolyparams"; case BasalforcingsIsThermalForcingEnum : return "BasalforcingsIsThermalForcing"; + case BasalforcingsLowercrustheatEnum : return "BasalforcingsLowercrustheat"; case BasalforcingsPicoAverageOverturningEnum : return "BasalforcingsPicoAverageOverturning"; case BasalforcingsPicoAverageSalinityEnum : return "BasalforcingsPicoAverageSalinity"; @@ -215,11 +216,12 @@ case FrictionVoidRatioEnum : return "FrictionVoidRatio"; case FrontalForcingsBasinIcefrontAreaEnum : return "FrontalForcingsBasinIcefrontArea"; - case FrontalForcingsARMAInitialTimeEnum : return "FrontalForcingsARMAInitialTime"; case FrontalForcingsARMATimestepEnum : return "FrontalForcingsARMATimestep"; case FrontalForcingsARMAarOrderEnum : return "FrontalForcingsARMAarOrder"; case FrontalForcingsARMAmaOrderEnum : return "FrontalForcingsARMAmaOrder"; - case FrontalForcingsARMAconstEnum : return "FrontalForcingsARMAconst"; - case FrontalForcingsARMAtrendEnum : return "FrontalForcingsARMAtrend"; + case FrontalForcingsARMAdatebreaksEnum : return "FrontalForcingsARMAdatebreaks"; + case FrontalForcingsARMApolyparamsEnum : return "FrontalForcingsARMApolyparams"; case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins"; + case FrontalForcingsNumberofBreaksEnum : return "FrontalForcingsNumberofBreaks"; + case FrontalForcingsNumberofParamsEnum : return "FrontalForcingsNumberofParams"; case FrontalForcingsParamEnum : return "FrontalForcingsParam"; case FrontalForcingsARMAarlagcoefsEnum : return "FrontalForcingsARMAarlagcoefs"; @@ -496,11 +498,8 @@ case SmbAccurefEnum : return "SmbAccuref"; case SmbAdThreshEnum : return "SmbAdThresh"; - case SmbARMAInitialTimeEnum : return "SmbARMAInitialTime"; case SmbARMATimestepEnum : return "SmbARMATimestep"; case SmbARMAarOrderEnum : return "SmbARMAarOrder"; case SmbARMAmaOrderEnum : return "SmbARMAmaOrder"; case SmbAveragingEnum : return "SmbAveraging"; - case SmbARMAconstEnum : return "SmbARMAconst"; - case SmbARMAtrendEnum : return "SmbARMAtrend"; case SmbDesfacEnum : return "SmbDesfac"; case SmbDpermilEnum : return "SmbDpermil"; @@ -536,9 +535,13 @@ case SmbLapseRatesEnum : return "SmbLapseRates"; case SmbNumBasinsEnum : return "SmbNumBasins"; + case SmbNumBreaksEnum : return "SmbNumBreaks"; case SmbNumElevationBinsEnum : return "SmbNumElevationBins"; + case SmbNumParamsEnum : return "SmbNumParams"; case SmbNumRequestedOutputsEnum : return "SmbNumRequestedOutputs"; case SmbPfacEnum : return "SmbPfac"; case SmbARMAarlagcoefsEnum : return "SmbARMAarlagcoefs"; + case SmbARMAdatebreaksEnum : return "SmbARMAdatebreaks"; case SmbARMAmalagcoefsEnum : return "SmbARMAmalagcoefs"; + case SmbARMApolyparamsEnum : return "SmbARMApolyparams"; case SmbRdlEnum : return "SmbRdl"; case SmbRefElevationEnum : return "SmbRefElevation"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27317) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27318) @@ -60,10 +60,7 @@ syn keyword juliaConstC BalancethicknessStabilizationEnum syn keyword juliaConstC BarystaticContributionsEnum -syn keyword juliaConstC BasalforcingsARMAInitialTimeEnum syn keyword juliaConstC BasalforcingsARMATimestepEnum syn keyword juliaConstC BasalforcingsARMAarOrderEnum syn keyword juliaConstC BasalforcingsARMAmaOrderEnum -syn keyword juliaConstC BasalforcingsARMAconstEnum -syn keyword juliaConstC BasalforcingsARMAtrendEnum syn keyword juliaConstC BasalforcingsBottomplumedepthEnum syn keyword juliaConstC BasalforcingsCrustthicknessEnum @@ -80,10 +77,14 @@ syn keyword juliaConstC BasalforcingsIsmip6TfDepthsEnum syn keyword juliaConstC BasalforcingsLinearNumBasinsEnum -syn keyword juliaConstC BasalforcingsLowercrustheatEnum +syn keyword juliaConstC BasalforcingsLinearNumBreaksEnum +syn keyword juliaConstC BasalforcingsLinearNumParamsEnum syn keyword juliaConstC BasalforcingsMantleconductivityEnum syn keyword juliaConstC BasalforcingsNusseltEnum syn keyword juliaConstC BasalforcingsARMAarlagcoefsEnum +syn keyword juliaConstC BasalforcingsARMAdatebreaksEnum syn keyword juliaConstC BasalforcingsARMAmalagcoefsEnum +syn keyword juliaConstC BasalforcingsARMApolyparamsEnum syn keyword juliaConstC BasalforcingsIsThermalForcingEnum +syn keyword juliaConstC BasalforcingsLowercrustheatEnum syn keyword juliaConstC BasalforcingsPicoAverageOverturningEnum syn keyword juliaConstC BasalforcingsPicoAverageSalinityEnum @@ -206,11 +207,12 @@ syn keyword juliaConstC FrictionVoidRatioEnum syn keyword juliaConstC FrontalForcingsBasinIcefrontAreaEnum -syn keyword juliaConstC FrontalForcingsARMAInitialTimeEnum syn keyword juliaConstC FrontalForcingsARMATimestepEnum syn keyword juliaConstC FrontalForcingsARMAarOrderEnum syn keyword juliaConstC FrontalForcingsARMAmaOrderEnum -syn keyword juliaConstC FrontalForcingsARMAconstEnum -syn keyword juliaConstC FrontalForcingsARMAtrendEnum +syn keyword juliaConstC FrontalForcingsARMAdatebreaksEnum +syn keyword juliaConstC FrontalForcingsARMApolyparamsEnum syn keyword juliaConstC FrontalForcingsNumberofBasinsEnum +syn keyword juliaConstC FrontalForcingsNumberofBreaksEnum +syn keyword juliaConstC FrontalForcingsNumberofParamsEnum syn keyword juliaConstC FrontalForcingsParamEnum syn keyword juliaConstC FrontalForcingsARMAarlagcoefsEnum @@ -311,4 +313,5 @@ syn keyword juliaConstC LoveMinIntegrationStepsEnum syn keyword juliaConstC LoveMaxIntegrationdrEnum +syn keyword juliaConstC LoveIntegrationSchemeEnum syn keyword juliaConstC LoveKernelsEnum syn keyword juliaConstC LoveMu0Enum @@ -486,11 +489,8 @@ syn keyword juliaConstC SmbAccurefEnum syn keyword juliaConstC SmbAdThreshEnum -syn keyword juliaConstC SmbARMAInitialTimeEnum syn keyword juliaConstC SmbARMATimestepEnum syn keyword juliaConstC SmbARMAarOrderEnum syn keyword juliaConstC SmbARMAmaOrderEnum syn keyword juliaConstC SmbAveragingEnum -syn keyword juliaConstC SmbARMAconstEnum -syn keyword juliaConstC SmbARMAtrendEnum syn keyword juliaConstC SmbDesfacEnum syn keyword juliaConstC SmbDpermilEnum @@ -526,9 +526,13 @@ syn keyword juliaConstC SmbLapseRatesEnum syn keyword juliaConstC SmbNumBasinsEnum +syn keyword juliaConstC SmbNumBreaksEnum syn keyword juliaConstC SmbNumElevationBinsEnum +syn keyword juliaConstC SmbNumParamsEnum syn keyword juliaConstC SmbNumRequestedOutputsEnum syn keyword juliaConstC SmbPfacEnum syn keyword juliaConstC SmbARMAarlagcoefsEnum +syn keyword juliaConstC SmbARMAdatebreaksEnum syn keyword juliaConstC SmbARMAmalagcoefsEnum +syn keyword juliaConstC SmbARMApolyparamsEnum syn keyword juliaConstC SmbRdlEnum syn keyword juliaConstC SmbRefElevationEnum @@ -932,4 +936,5 @@ syn keyword juliaConstC SealevelUNorthEsaEnum syn keyword juliaConstC SealevelchangeIndicesEnum +syn keyword juliaConstC SealevelchangeConvolutionVerticesEnum syn keyword juliaConstC SealevelchangeAlphaIndexEnum syn keyword juliaConstC SealevelchangeAzimuthIndexEnum @@ -950,4 +955,5 @@ syn keyword juliaConstC SealevelchangeViscousNEnum syn keyword juliaConstC SealevelchangeViscousEEnum +syn keyword juliaConstC CouplingTransferCountEnum syn keyword juliaConstC SedimentHeadEnum syn keyword juliaConstC SedimentHeadOldEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27317) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27318) @@ -69,10 +69,7 @@ else if (strcmp(name,"BalancethicknessStabilization")==0) return BalancethicknessStabilizationEnum; else if (strcmp(name,"BarystaticContributions")==0) return BarystaticContributionsEnum; - else if (strcmp(name,"BasalforcingsARMAInitialTime")==0) return BasalforcingsARMAInitialTimeEnum; else if (strcmp(name,"BasalforcingsARMATimestep")==0) return BasalforcingsARMATimestepEnum; else if (strcmp(name,"BasalforcingsARMAarOrder")==0) return BasalforcingsARMAarOrderEnum; else if (strcmp(name,"BasalforcingsARMAmaOrder")==0) return BasalforcingsARMAmaOrderEnum; - else if (strcmp(name,"BasalforcingsARMAconst")==0) return BasalforcingsARMAconstEnum; - else if (strcmp(name,"BasalforcingsARMAtrend")==0) return BasalforcingsARMAtrendEnum; else if (strcmp(name,"BasalforcingsBottomplumedepth")==0) return BasalforcingsBottomplumedepthEnum; else if (strcmp(name,"BasalforcingsCrustthickness")==0) return BasalforcingsCrustthicknessEnum; @@ -89,10 +86,14 @@ else if (strcmp(name,"BasalforcingsIsmip6TfDepths")==0) return BasalforcingsIsmip6TfDepthsEnum; else if (strcmp(name,"BasalforcingsLinearNumBasins")==0) return BasalforcingsLinearNumBasinsEnum; - else if (strcmp(name,"BasalforcingsLowercrustheat")==0) return BasalforcingsLowercrustheatEnum; + else if (strcmp(name,"BasalforcingsLinearNumBreaks")==0) return BasalforcingsLinearNumBreaksEnum; + else if (strcmp(name,"BasalforcingsLinearNumParams")==0) return BasalforcingsLinearNumParamsEnum; else if (strcmp(name,"BasalforcingsMantleconductivity")==0) return BasalforcingsMantleconductivityEnum; else if (strcmp(name,"BasalforcingsNusselt")==0) return BasalforcingsNusseltEnum; else if (strcmp(name,"BasalforcingsARMAarlagcoefs")==0) return BasalforcingsARMAarlagcoefsEnum; + else if (strcmp(name,"BasalforcingsARMAdatebreaks")==0) return BasalforcingsARMAdatebreaksEnum; else if (strcmp(name,"BasalforcingsARMAmalagcoefs")==0) return BasalforcingsARMAmalagcoefsEnum; + else if (strcmp(name,"BasalforcingsARMApolyparams")==0) return BasalforcingsARMApolyparamsEnum; else if (strcmp(name,"BasalforcingsIsThermalForcing")==0) return BasalforcingsIsThermalForcingEnum; + else if (strcmp(name,"BasalforcingsLowercrustheat")==0) return BasalforcingsLowercrustheatEnum; else if (strcmp(name,"BasalforcingsPicoAverageOverturning")==0) return BasalforcingsPicoAverageOverturningEnum; else if (strcmp(name,"BasalforcingsPicoAverageSalinity")==0) return BasalforcingsPicoAverageSalinityEnum; @@ -136,9 +137,9 @@ else if (strcmp(name,"ConstantsNewtonGravity")==0) return ConstantsNewtonGravityEnum; else if (strcmp(name,"ConstantsReferencetemperature")==0) return ConstantsReferencetemperatureEnum; - else if (strcmp(name,"ConstantsYts")==0) return ConstantsYtsEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"ControlInputSizeM")==0) return ControlInputSizeMEnum; + if (strcmp(name,"ConstantsYts")==0) return ConstantsYtsEnum; + else if (strcmp(name,"ControlInputSizeM")==0) return ControlInputSizeMEnum; else if (strcmp(name,"ControlInputSizeN")==0) return ControlInputSizeNEnum; else if (strcmp(name,"ControlInputInterpolation")==0) return ControlInputInterpolationEnum; @@ -218,11 +219,12 @@ else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum; else if (strcmp(name,"FrontalForcingsBasinIcefrontArea")==0) return FrontalForcingsBasinIcefrontAreaEnum; - else if (strcmp(name,"FrontalForcingsARMAInitialTime")==0) return FrontalForcingsARMAInitialTimeEnum; else if (strcmp(name,"FrontalForcingsARMATimestep")==0) return FrontalForcingsARMATimestepEnum; else if (strcmp(name,"FrontalForcingsARMAarOrder")==0) return FrontalForcingsARMAarOrderEnum; else if (strcmp(name,"FrontalForcingsARMAmaOrder")==0) return FrontalForcingsARMAmaOrderEnum; - else if (strcmp(name,"FrontalForcingsARMAconst")==0) return FrontalForcingsARMAconstEnum; - else if (strcmp(name,"FrontalForcingsARMAtrend")==0) return FrontalForcingsARMAtrendEnum; + else if (strcmp(name,"FrontalForcingsARMAdatebreaks")==0) return FrontalForcingsARMAdatebreaksEnum; + else if (strcmp(name,"FrontalForcingsARMApolyparams")==0) return FrontalForcingsARMApolyparamsEnum; else if (strcmp(name,"FrontalForcingsNumberofBasins")==0) return FrontalForcingsNumberofBasinsEnum; + else if (strcmp(name,"FrontalForcingsNumberofBreaks")==0) return FrontalForcingsNumberofBreaksEnum; + else if (strcmp(name,"FrontalForcingsNumberofParams")==0) return FrontalForcingsNumberofParamsEnum; else if (strcmp(name,"FrontalForcingsParam")==0) return FrontalForcingsParamEnum; else if (strcmp(name,"FrontalForcingsARMAarlagcoefs")==0) return FrontalForcingsARMAarlagcoefsEnum; @@ -258,10 +260,10 @@ else if (strcmp(name,"HydrologydcEplflipLock")==0) return HydrologydcEplflipLockEnum; else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum; - else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum; - else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum; + if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum; + else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum; + else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum; else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum; else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum; @@ -381,10 +383,10 @@ else if (strcmp(name,"OceanGridY")==0) return OceanGridYEnum; else if (strcmp(name,"OutputBufferPointer")==0) return OutputBufferPointerEnum; - else if (strcmp(name,"OutputBufferSizePointer")==0) return OutputBufferSizePointerEnum; - else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; + if (strcmp(name,"OutputBufferSizePointer")==0) return OutputBufferSizePointerEnum; + else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum; + else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum; @@ -504,15 +506,12 @@ else if (strcmp(name,"SmbAccugrad")==0) return SmbAccugradEnum; else if (strcmp(name,"SmbAccuref")==0) return SmbAccurefEnum; - else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum; - else if (strcmp(name,"SmbARMAInitialTime")==0) return SmbARMAInitialTimeEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum; + if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum; + else if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum; else if (strcmp(name,"SmbARMAarOrder")==0) return SmbARMAarOrderEnum; else if (strcmp(name,"SmbARMAmaOrder")==0) return SmbARMAmaOrderEnum; else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum; - else if (strcmp(name,"SmbARMAconst")==0) return SmbARMAconstEnum; - else if (strcmp(name,"SmbARMAtrend")==0) return SmbARMAtrendEnum; else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; @@ -548,9 +547,13 @@ else if (strcmp(name,"SmbLapseRates")==0) return SmbLapseRatesEnum; else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum; + else if (strcmp(name,"SmbNumBreaks")==0) return SmbNumBreaksEnum; else if (strcmp(name,"SmbNumElevationBins")==0) return SmbNumElevationBinsEnum; + else if (strcmp(name,"SmbNumParams")==0) return SmbNumParamsEnum; else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum; else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum; else if (strcmp(name,"SmbARMAarlagcoefs")==0) return SmbARMAarlagcoefsEnum; + else if (strcmp(name,"SmbARMAdatebreaks")==0) return SmbARMAdatebreaksEnum; else if (strcmp(name,"SmbARMAmalagcoefs")==0) return SmbARMAmalagcoefsEnum; + else if (strcmp(name,"SmbARMApolyparams")==0) return SmbARMApolyparamsEnum; else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum; else if (strcmp(name,"SmbRefElevation")==0) return SmbRefElevationEnum; @@ -626,11 +629,11 @@ else if (strcmp(name,"TransientIsdebris")==0) return TransientIsdebrisEnum; else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum; - else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum; - else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; - else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; + if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum; + else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; + else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; + else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; @@ -749,11 +752,11 @@ else if (strcmp(name,"StrOld")==0) return StrOldEnum; else if (strcmp(name,"Str")==0) return StrEnum; - else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; - else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; - else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; + if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; + else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; + else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; + else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; @@ -872,11 +875,11 @@ else if (strcmp(name,"MaskOceanLevelset")==0) return MaskOceanLevelsetEnum; else if (strcmp(name,"MaskIceLevelset")==0) return MaskIceLevelsetEnum; - else if (strcmp(name,"MaskIceRefLevelset")==0) return MaskIceRefLevelsetEnum; - else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum; - else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum; + if (strcmp(name,"MaskIceRefLevelset")==0) return MaskIceRefLevelsetEnum; + else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum; + else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum; + else if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum; else if (strcmp(name,"MaterialsRheologyE")==0) return MaterialsRheologyEEnum; else if (strcmp(name,"MaterialsRheologyEbar")==0) return MaterialsRheologyEbarEnum; @@ -995,11 +998,11 @@ else if (strcmp(name,"SmbAccumulatedMelt")==0) return SmbAccumulatedMeltEnum; else if (strcmp(name,"SmbAccumulatedPrecipitation")==0) return SmbAccumulatedPrecipitationEnum; - else if (strcmp(name,"SmbAccumulatedRain")==0) return SmbAccumulatedRainEnum; - else if (strcmp(name,"SmbAccumulatedRefreeze")==0) return SmbAccumulatedRefreezeEnum; - else if (strcmp(name,"SmbAccumulatedRunoff")==0) return SmbAccumulatedRunoffEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SmbA")==0) return SmbAEnum; + if (strcmp(name,"SmbAccumulatedRain")==0) return SmbAccumulatedRainEnum; + else if (strcmp(name,"SmbAccumulatedRefreeze")==0) return SmbAccumulatedRefreezeEnum; + else if (strcmp(name,"SmbAccumulatedRunoff")==0) return SmbAccumulatedRunoffEnum; + else if (strcmp(name,"SmbA")==0) return SmbAEnum; else if (strcmp(name,"SmbAdiff")==0) return SmbAdiffEnum; else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum; @@ -1118,11 +1121,11 @@ else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; - else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; - else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum; - else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum; + if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; + else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum; + else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum; + else if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum; else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; @@ -1241,11 +1244,11 @@ else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum; - else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; - else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; - else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum; else stage=11; } if(stage==11){ - if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum; + if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; + else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; + else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum; + else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum; else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum; else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum; @@ -1364,11 +1367,11 @@ else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; - else if (strcmp(name,"DataSet")==0) return DataSetEnum; - else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; - else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"DebrisAnalysis")==0) return DebrisAnalysisEnum; + if (strcmp(name,"DataSet")==0) return DataSetEnum; + else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; + else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; + else if (strcmp(name,"DebrisAnalysis")==0) return DebrisAnalysisEnum; else if (strcmp(name,"DebrisSolution")==0) return DebrisSolutionEnum; else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; @@ -1487,11 +1490,11 @@ else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; else if (strcmp(name,"Loads")==0) return LoadsEnum; - else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; - else if (strcmp(name,"LoveHf")==0) return LoveHfEnum; - else if (strcmp(name,"LoveHt")==0) return LoveHtEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum; + if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; + else if (strcmp(name,"LoveHf")==0) return LoveHfEnum; + else if (strcmp(name,"LoveHt")==0) return LoveHtEnum; + else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum; else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum; else if (strcmp(name,"LoveKf")==0) return LoveKfEnum; @@ -1610,11 +1613,11 @@ else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; - else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; - else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; - else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; else stage=14; } if(stage==14){ - if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; + if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; + else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; + else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; + else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; else if (strcmp(name,"Scaled")==0) return ScaledEnum; Index: /issm/trunk-jpl/src/m/classes/SMBarma.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBarma.m (revision 27317) +++ /issm/trunk-jpl/src/m/classes/SMBarma.m (revision 27318) @@ -7,7 +7,6 @@ properties (SetAccess=public) num_basins = 0; - const = NaN; - trend = NaN; - arma_initialtime = 0; + num_breaks = 0; + num_params = 0; arma_timestep = 0; ar_order = 0; @@ -15,4 +14,6 @@ ma_order = 0; malag_coefs = NaN; + polynomialparams = NaN; + datebreaks = NaN; basin_id = NaN; lapserates = NaN; @@ -53,8 +54,4 @@ disp(' smb.ma_order (order of moving-average model) not specified: order of moving-average model set to 0'); end - if (self.arma_initialtime==0) - self.ar_initialtime = md.timestepping.start_time; %ARMA model has no prescribed initial time - disp(' smb.arma_initialtime (initial time in the ARMA model parameterization) not specified: set to md.timestepping.start_time'); - end if (self.arma_timestep==0) self.arma_timestep = md.timestepping.time_step; %ARMA model has no prescribed time step @@ -77,15 +74,33 @@ if ismember('MasstransportAnalysis',analyses), + nbas = md.smb.num_basins; + nprm = md.smb.num_params; + nbrk = md.smb.num_breaks; md = checkfield(md,'fieldname','smb.num_basins','numel',1,'NaN',1,'Inf',1,'>',0); + md = checkfield(md,'fieldname','smb.num_params','numel',1,'NaN',1,'Inf',1,'>',0); + md = checkfield(md,'fieldname','smb.num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','smb.basin_id','Inf',1,'>=',0,'<=',md.smb.num_basins,'size',[md.mesh.numberofelements,1]); - md = checkfield(md,'fieldname','smb.const','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); %scheme fails if passed as column vector - md = checkfield(md,'fieldname','smb.trend','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); %scheme fails if passed as column vector + if(nbas>1 && nbrk>=1 && nprm>1) + md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm); + elseif(nbas==1) + md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm); + elseif(nbrk==0) + md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm); + elseif(nprm==1) + md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1],'numel',nbas*(nbrk+1)*nprm); + end md = checkfield(md,'fieldname','smb.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','smb.ma_order','numel',1,'NaN',1,'Inf',1,'>=',0); - md = checkfield(md,'fieldname','smb.arma_initialtime','numel',1,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','smb.arma_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %arma time step cannot be finer than ISSM timestep md = checkfield(md,'fieldname','smb.arlag_coefs','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]); md = checkfield(md,'fieldname','smb.malag_coefs','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ma_order]); - + + if(nbrk>0) + md = checkfield(md,'fieldname','smb.datebreaks','NaN',1,'Inf',1,'size',[nbas,nbrk]); + elseif(numel(md.smb.datebreaks)==0 || all(isnan(md.smb.datebreaks))) + ; + else + error('md.smb.num_breaks is 0 but md.smb.datebreaks is not empty'); + end if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1) md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); @@ -114,9 +129,11 @@ fielddisplay(self,'num_basins','number of different basins [unitless]'); fielddisplay(self,'basin_id','basin number assigned to each element [unitless]'); - fielddisplay(self,'const','basin-specific constant values [m ice eq./yr]'); - fielddisplay(self,'trend','basin-specific trend values [m ice eq. yr^(-2)]'); + fielddisplay(self,'num_breaks','number of different breakpoints in the piecewise-polynomial (separating num_breaks+1 periods)'); + fielddisplay(self,'num_params','number of different parameters in the piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)'); + fielddisplay(self,'polynomialparams','coefficients for the polynomial (const,trend,quadratic,etc.),dim1 for basins,dim2 for periods,dim3 for orders'); + disp(sprintf('%51s ex: polyparams=cat(3,intercepts,trendlinearcoefs,trendquadraticcoefs)',' ')); + fielddisplay(self,'datebreaks','dates at which the breakpoints in the piecewise polynomial occur (1 row per basin) [yr]'); fielddisplay(self,'ar_order','order of the autoregressive model [unitless]'); fielddisplay(self,'ma_order','order of the moving-average model [unitless]'); - fielddisplay(self,'arma_initialtime','initial time assumed in the ARMA model parameterization [yr]'); fielddisplay(self,'arma_timestep','time resolution of the ARMA model [yr]'); fielddisplay(self,'arlag_coefs','basin-specific vectors of AR lag coefficients [unitless]'); @@ -136,4 +153,7 @@ yts=md.constants.yts; + nbas = md.smb.num_basins; + nprm = md.smb.num_params; + nper = md.smb.num_breaks+1; templapserates = md.smb.lapserates; @@ -163,15 +183,59 @@ [nbas,nbins] = size(templapserates); + % Scale the parameters % + polyparamsScaled = md.smb.polynomialparams; + polyparams2dScaled = zeros(nbas,nper*nprm); + if(nprm>1) + % Case 3D % + if(nbas>1 && nper>1) + for(ii=[1:nprm]) + polyparamsScaled(:,:,ii) = polyparamsScaled(:,:,ii)*((1/yts)^(ii)); + end + % Fit in 2D array % + for(ii=[1:nprm]) + jj = 1+(ii-1)*nper; + polyparams2dScaled(:,jj:jj+nper-1) = polyparamsScaled(:,:,ii); + end + % Case 2D and higher-order params at increasing row index % + elseif(nbas==1) + for(ii=[1:nprm]) + polyparamsScaled(ii,:) = polyparamsScaled(ii,:)*((1/yts)^(ii)); + end + % Fit in row array % + for(ii=[1:nprm]) + jj = 1+(ii-1)*nper; + polyparams2dScaled(1,jj:jj+nper-1) = polyparamsScaled(ii,:); + end + % Case 2D and higher-order params at incrasing column index % + elseif(nper==1) + for(ii=[1:nprm]) + polyparamsScaled(:,ii) = polyparamsScaled(:,ii)*((1/yts)^(ii)); + end + % 2D array is already in correct format % + polyparams2dScaled = polyparamsScaled; + end + else + polyparamsScaled = polyparamsScaled*(1/yts); + % 2D array is already in correct format % + polyparams2dScaled = polyparamsScaled; + end + if(nper==1) %a single period (no break date) + dbreaks = zeros(nbas,1); %dummy + else + dbreaks = md.smb.datebreaks; + end + WriteData(fid,prefix,'name','md.smb.model','data',13,'format','Integer'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','num_basins','format','Integer'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','num_params','format','Integer'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','num_breaks','format','Integer'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_order','format','Integer'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','ma_order','format','Integer'); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','arma_initialtime','format','Double','scale',yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','arma_timestep','format','Double','scale',yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','basin_id','data',self.basin_id-1,'name','md.smb.basin_id','format','IntMat','mattype',2); %0-indexed - WriteData(fid,prefix,'object',self,'class','smb','fieldname','const','format','DoubleMat','name','md.smb.const','scale',1./yts,'yts',yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','trend','format','DoubleMat','name','md.smb.trend','scale',1./(yts^2),'yts',yts); + WriteData(fid,prefix,'data',polyparams2dScaled,'name','md.smb.polynomialparams','format','DoubleMat'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','arlag_coefs','format','DoubleMat','name','md.smb.arlag_coefs','yts',yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','malag_coefs','format','DoubleMat','name','md.smb.malag_coefs','yts',yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','malag_coefs','format','DoubleMat','name','md.smb.malag_coefs','yts',yts); + WriteData(fid,prefix,'data',dbreaks,'name','md.smb.datebreaks','format','DoubleMat','scale',yts); WriteData(fid,prefix,'data',templapserates,'format','DoubleMat','name','md.smb.lapserates','scale',1./yts,'yts',yts); WriteData(fid,prefix,'data',tempelevationbins,'format','DoubleMat','name','md.smb.elevationbins'); Index: /issm/trunk-jpl/src/m/classes/SMBarma.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBarma.py (revision 27317) +++ /issm/trunk-jpl/src/m/classes/SMBarma.py (revision 27318) @@ -16,9 +16,8 @@ def __init__(self, *args): # {{{ self.num_basins = 0 - self.const = np.nan - self.trend = np.nan + self.num_params = 0 + self.num_breaks = 0 + self.polynomialparams = np.nan self.ar_order = 0 - self.arma_initialtime = 0 - self.arma_timestep = 0 self.arlag_coefs = np.nan self.malag_coefs = np.nan @@ -42,9 +41,10 @@ s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]')) s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'const', 'basin-specific constant values [m ice eq./yr]')) - s += '{}\n'.format(fielddisplay(self, 'trend', 'basin-specific trend values [m ice eq. yr^(-2)]')) + s += '{}\n'.format(fielddisplay(self, 'num_breaks', 'number of different breakpoints in the piecewise-polynomial (separating num_breaks+1 periods)')) + s += '{}\n'.format(fielddisplay(self, 'num_params', 'number of different parameters in the piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)')) + s += '{}\n'.format(fielddisplay(self, 'polynomialparams', 'coefficients for the polynomial (const,trend,quadratic,etc.),dim1 for basins,dim2 for periods,dim3 for orders, ex: polyparams=cat(num_params,intercepts,trendlinearcoefs,trendquadraticcoefs)')) + s += '{}\n'.format(fielddisplay(self, 'datebreaks', 'dates at which the breakpoints in the piecewise polynomial occur (1 row per basin) [yr]')) s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]')) s += '{}\n'.format(fielddisplay(self, 'ma_order', 'order of the moving-average model [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'arma_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]')) s += '{}\n'.format(fielddisplay(self, 'arma_timestep', 'time resolution of the ARMA model [yr]')) s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of AR lag coefficients [unitless]')) @@ -83,7 +83,4 @@ self.arlag_coefs = np.zeros((self.num_basins, self.ar_order)) # Autoregression coefficients all set to 0 print(' smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0') - if self.arma_initialtime == 0: - self.arma_initialtime = md.timestepping.start_time # ARMA model has no prescribed initial time - print(' smb.arma_initialtime (initial time in the ARMA model parameterization) not specified: set to md.timestepping.start_time') if self.arma_timestep == 0: self.arma_timestep = md.timestepping.time_step # ARMA model has no prescribed time step @@ -100,17 +97,32 @@ def checkconsistency(self, md, solution, analyses): # {{{ if 'MasstransportAnalysis' in analyses: + nbas = md.smb.num_basins; + nprm = md.smb.num_params; + nbrk = md.smb.num_breaks; md = checkfield(md, 'fieldname', 'smb.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) + md = checkfield(md, 'fieldname', 'smb.num_params', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) + md = checkfield(md, 'fieldname', 'smb.num_breaks', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) md = checkfield(md, 'fieldname', 'smb.basin_id', 'Inf', 1, '>=', 0, '<=', md.smb.num_basins, 'size', [md.mesh.numberofelements]) - if len(np.shape(self.const)) == 1: - self.const = np.array([self.const]) - self.trend = np.array([self.trend]) - md = checkfield(md, 'fieldname', 'smb.const', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) # Scheme fails if passed as column vector - md = checkfield(md, 'fieldname', 'smb.trend', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) # Scheme fails if passed as column vector; NOTE: As opposed to MATLAB implementation, pass list + if len(np.shape(self.polynomialparams)) == 1: + self.polynomialparams = np.array([[self.polynomialparams]]) + if(nbas>1 and nbrk>=1 and nprm>1): + md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm) + elif(nbas==1): + md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm) + elif(nbrk==0): + md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm) + elif(nprm==1): + md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk],'numel',nbas*(nbrk+1)*nprm) md = checkfield(md, 'fieldname', 'smb.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) - md = checkfield(md, 'fieldname', 'smb.arma_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'smb.arma_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep md = checkfield(md, 'fieldname', 'smb.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order]) md = checkfield(md, 'fieldname', 'smb.malag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ma_order]) - + if(nbrk>0): + md = checkfield(md, 'fieldname', 'smb.datebreaks', 'NaN', 1, 'Inf', 1, 'size', [nbas,nbrk]) + elif(np.size(md.smb.datebreaks)==0 or np.all(np.isnan(md.smb.datebreaks))): + pass + else: + raise RuntimeError('md.smb.num_breaks is 0 but md.smb.datebreaks is not empty') + if(np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1): if len(np.shape(self.refelevation)) == 1: @@ -149,8 +161,44 @@ def marshall(self, prefix, md, fid): # {{{ yts = md.constants.yts - + nbas = md.smb.num_basins; + nprm = md.smb.num_params; + nper = md.smb.num_breaks+1; templapserates = np.copy(md.smb.lapserates) tempelevationbins = np.copy(md.smb.elevationbins) temprefelevation = np.copy(md.smb.refelevation) + # Scale the parameters # + polyparamsScaled = np.copy(md.smb.polynomialparams) + polyparams2dScaled = np.zeros((nbas,nper*nprm)) + if(nprm>1): + # Case 3D # + if(nbas>1 and nper>1): + for ii in range(nprm): + polyparamsScaled[:,:,ii] = polyparamsScaled[:,:,ii]*(1/yts)**(ii+1) + # Fit in 2D array # + for ii in range(nprm): + polyparams2dScaled[:,ii*nper:(ii+1)*nper] = 1*polyparamsScaled[:,:,ii] + # Case 2D and higher-order params at increasing row index # + elif(nbas==1): + for ii in range(nprm): + polyparamsScaled[ii,:] = polyparamsScaled[ii,:]*(1/yts)**(ii+1) + # Fit in row array # + for ii in range(nprm): + polyparams2dScaled[0,ii*nper:(ii+1)*nper] = 1*polyparamsScaled[ii,:] + # Case 2D and higher-order params at incrasing column index # + elif(nper==1): + for ii in range(nprm): + polyparamsScaled[:,ii] = polyparamsScaled[:,ii]*(1/yts)**(ii+1) + # 2D array is already in correct format # + polyparams2dScaled = np.copy(polyparamsScaled) + else: + polyparamsScaled = polyparamsScaled*(1/yts) + # 2D array is already in correct format # + polyparams2dScaled = np.copy(polyparamsScaled) + + if(nper==1): + dbreaks = np.zeros((nbas,1)) + else: + dbreaks = np.copy(md.smb.datebreaks) + if(np.any(np.isnan(md.smb.lapserates))): templapserates = np.zeros((md.smb.num_basins,2)) @@ -172,13 +220,14 @@ WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 13, 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'num_basins', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'num_breaks', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'num_params', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ar_order', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ma_order', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'arma_initialtime', 'format', 'Double', 'scale', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'arma_timestep', 'format', 'Double', 'scale', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'basin_id', 'data', self.basin_id - 1, 'name', 'md.smb.basin_id', 'format', 'IntMat', 'mattype', 2) # 0-indexed - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'const', 'format', 'DoubleMat', 'name', 'md.smb.const', 'scale', 1 / yts, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'trend', 'format', 'DoubleMat', 'name', 'md.smb.trend', 'scale', 1 / (yts ** 2), 'yts', yts) + WriteData(fid, prefix, 'data', polyparams2dScaled, 'name', 'md.smb.polynomialparams', 'format', 'DoubleMat') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.smb.arlag_coefs', 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'malag_coefs', 'format', 'DoubleMat', 'name', 'md.smb.malag_coefs', 'yts', yts) + WriteData(fid, prefix, 'data', dbreaks, 'name', 'md.smb.datebreaks', 'format', 'DoubleMat','scale',yts) WriteData(fid, prefix, 'data', templapserates, 'name', 'md.smb.lapserates', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts) WriteData(fid, prefix, 'data', tempelevationbins, 'name', 'md.smb.elevationbins', 'format', 'DoubleMat') Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m (revision 27317) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m (revision 27318) @@ -7,9 +7,10 @@ properties (SetAccess=public) num_basins = 0; - const = NaN; - trend = NaN; + num_params = 0; + num_breaks = 0; + polynomialparams = NaN; + datebreaks = NaN; ar_order = 0; ma_order = 0; - arma_initialtime = 0; arma_timestep = 0; arlag_coefs = NaN; @@ -54,15 +55,33 @@ if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end + nbas = md.frontalforcings.num_basins; + nprm = md.frontalforcings.num_params; + nbrk = md.frontalforcings.num_breaks; md = checkfield(md,'fieldname','frontalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0); + md = checkfield(md,'fieldname','frontalforcings.num_params','numel',1,'NaN',1,'Inf',1,'>',0); + md = checkfield(md,'fieldname','frontalforcings.num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','frontalforcings.basin_id','Inf',1,'>=',0,'<=',md.frontalforcings.num_basins,'size',[md.mesh.numberofelements 1]); md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'NaN',1,'Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','frontalforcings.const','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins); - md = checkfield(md,'fieldname','frontalforcings.trend','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins); - md = checkfield(md,'fieldname','frontalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); + if(nbas>1 && nbrk>=1 && nprm>1) + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm); + elseif(nbas==1) + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm); + elseif(nbrk==0) + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm); + elseif(nprm==1) + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1],'numel',nbas*(nbrk+1)*nprm); + end + md = checkfield(md,'fieldname','frontalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','frontalforcings.ma_order','numel',1,'NaN',1,'Inf',1,'>=',0); - md = checkfield(md,'fieldname','frontalforcings.arma_initialtime','numel',1,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','frontalforcings.arma_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %ARMA time step cannot be finer than ISSM timestep md = checkfield(md,'fieldname','frontalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]); md = checkfield(md,'fieldname','frontalforcings.malag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ma_order]); + if(nbrk>0) + md = checkfield(md,'fieldname','frontalforcings.datebreaks','NaN',1,'Inf',1,'size',[nbas,nbrk]); + elseif(numel(md.frontalforcings.datebreaks)==0 || all(isnan(md.frontalforcings.datebreaks))) + ; + else + error('md.frontalforcings.num_breaks is 0 but md.frontalforcings.datebreaks is not empty'); + end if(any(~isnan(md.frontalforcings.monthly_effects))) md = checkfield(md,'fieldname','frontalforcings.monthly_effects','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,12]); @@ -75,12 +94,14 @@ function disp(self) % {{{ disp(sprintf(' Frontalforcings parameters:')); - fielddisplay(self,'num_basins','number of different basins [unitless]'); + fielddisplay(self,'num_basins','number of different basins'); fielddisplay(self,'basin_id','basin number assigned to each element [unitless]'); fielddisplay(self,'subglacial_discharge','sum of subglacial discharge for each basin [m/d]'); - fielddisplay(self,'const','basin-specific constant term [∘C]'); - fielddisplay(self,'trend','basin-specific trend values [∘C yr^(-1)]'); + fielddisplay(self,'num_breaks','number of different breakpoints in the piecewise-polynomial (separating num_breaks+1 periods)'); + fielddisplay(self,'num_params','number of different parameters in the piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)'); + fielddisplay(self,'polynomialparams','coefficients for the polynomial (const,trend,quadratic,etc.),dim1 for basins,dim2 for periods,dim3 for orders'); + disp(sprintf('%51s ex: polyparams=cat(3,intercepts,trendlinearcoefs,trendquadraticcoefs)',' ')); + fielddisplay(self,'datebreaks','dates at which the breakpoints in the piecewise polynomial occur (1 row per basin) [yr]'); fielddisplay(self,'ar_order','order of the autoregressive model [unitless]'); fielddisplay(self,'ma_order','order of the moving-average model [unitless]'); - fielddisplay(self,'arma_initialtime','initial time assumed in the ARMA model parameterization [yr]'); fielddisplay(self,'arma_timestep','time resolution of the autoregressive model [yr]'); fielddisplay(self,'arlag_coefs','basin-specific vectors of AR lag coefficients [unitless]'); @@ -90,4 +111,43 @@ function marshall(self,prefix,md,fid) % {{{ yts=md.constants.yts; + nbas = md.frontalforcings.num_basins; + nprm = md.frontalforcings.num_params; + nper = md.frontalforcings.num_breaks+1; + % Scale the parameters % + polyparamsScaled = md.frontalforcings.polynomialparams; + polyparams2dScaled = zeros(nbas,nper*nprm); + if(nprm>1) + % Case 3D % + if(nbas>1 && nper>1) + for(ii=[1:nprm]) + polyparamsScaled(:,:,ii) = polyparamsScaled(:,:,ii)*((1/yts)^(ii-1)); + end + % Fit in 2D array % + for(ii=[1:nprm]) + jj = 1+(ii-1)*nper; + polyparams2dScaled(:,jj:jj+nper-1) = polyparamsScaled(:,:,ii); + end + % Case 2D and higher-order params at increasing row index % + elseif(nbas==1) + for(ii=[1:nprm]) + polyparamsScaled(ii,:) = polyparamsScaled(ii,:)*((1/yts)^(ii-1)); + end + % Fit in row array % + for(ii=[1:nprm]) + jj = 1+(ii-1)*nper; + polyparams2dScaled(1,jj:jj+nper-1) = polyparamsScaled(ii,:); + end + % Case 2D and higher-order params at incrasing column index % + elseif(nper==1) + for(ii=[1:nprm]) + polyparamsScaled(:,ii) = polyparamsScaled(:,ii)*((1/yts)^(ii-1)); + end + % 2D array is already in correct format % + polyparams2dScaled = polyparamsScaled; + end + else + % 2D array is already in correct format and no need for scaling % + polyparams2dScaled = polyparamsScaled; + end if(any(isnan(md.frontalforcings.monthly_effects))) %monthly effects not provided, set to 0 meffects = zeros(md.frontalforcings.num_basins,12); @@ -95,17 +155,23 @@ meffects = md.frontalforcings.monthly_effects; end + if(nper==1) %a single period (no break date) + dbreaks = zeros(nbas,1); %dummy + else + dbreaks = md.frontalforcings.datebreaks; + end WriteData(fid,prefix,'name','md.frontalforcings.parameterization','data',3,'format','Integer'); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','num_basins','format','Integer'); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','num_breaks','format','Integer'); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','num_params','format','Integer'); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','subglacial_discharge','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_order','format','Integer'); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ma_order','format','Integer'); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','arma_initialtime','format','Double','scale',yts); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','arma_timestep','format','Double','scale',yts); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','basin_id','data',self.basin_id-1,'name','md.frontalforcings.basin_id','format','IntMat','mattype',2); %0-indexed - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','const','format','DoubleMat','name','md.frontalforcings.const'); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','trend','format','DoubleMat','name','md.frontalforcings.trend','scale',1./yts,'yts',yts); + WriteData(fid,prefix,'data',polyparams2dScaled,'name','md.frontalforcings.polynomialparams','format','DoubleMat'); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','arlag_coefs','format','DoubleMat','name','md.frontalforcings.arlag_coefs','yts',yts); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','malag_coefs','format','DoubleMat','name','md.frontalforcings.malag_coefs','yts',yts); + WriteData(fid,prefix,'data',dbreaks,'name','md.frontalforcings.datebreaks','format','DoubleMat','scale',yts); WriteData(fid,prefix,'data',meffects,'name','md.frontalforcings.monthly_effects','format','DoubleMat'); end % }}} Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py (revision 27317) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py (revision 27318) @@ -16,8 +16,10 @@ def __init__(self, *args): # {{{ self.num_basins = 0 - self.const = np.nan - self.trend = np.nan + self.num_params = 0 + self.num_breaks = 0 + self.polynomialparams = np.nan + self.datebreaks = np.nan self.ar_order = 0 - self.arma_initialtime = 0 + self.ma_order = 0 self.arma_timestep = 0 self.arlag_coefs = np.nan @@ -37,9 +39,10 @@ s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]')) s += '{}\n'.format(fielddisplay(self, 'subglacial_discharge', 'sum of subglacial discharge for each basin [m/d]')) - s += '{}\n'.format(fielddisplay(self, 'const', 'basin-specific constant values [°C]')) - s += '{}\n'.format(fielddisplay(self, 'trend', 'basin-specific trend values [°C yr^(-1)]')) + s += '{}\n'.format(fielddisplay(self, 'num_breaks', 'number of different breakpoints in the piecewise-polynomial (separating num_breaks+1 periods)')) + s += '{}\n'.format(fielddisplay(self, 'num_params', 'number of different parameters in the piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)')) + s += '{}\n'.format(fielddisplay(self, 'polynomialparams', 'coefficients for the polynomial (const,trend,quadratic,etc.),dim1 for basins,dim2 for periods,dim3 for orders, ex: polyparams=cat(num_params,intercepts,trendlinearcoefs,trendquadraticcoefs)')) + s += '{}\n'.format(fielddisplay(self, 'datebreaks', 'dates at which the breakpoints in the piecewise polynomial occur (1 row per basin) [yr]')) s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]')) s += '{}\n'.format(fielddisplay(self, 'ma_order', 'order of the moving-average model [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'arma_initialtime', 'initial time assumed in the ARMA model parameterization [yr]')) s += '{}\n'.format(fielddisplay(self, 'arma_timestep', 'time resolution of the ARMA model [yr]')) s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of AR lag coefficients [unitless]')) @@ -63,18 +66,33 @@ return md + nbas = md.frontalforcings.num_basins; + nprm = md.frontalforcings.num_params; + nbrk = md.frontalforcings.num_breaks; md = checkfield(md, 'fieldname', 'frontalforcings.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) + md = checkfield(md, 'fieldname', 'frontalforcings.num_params', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) + md = checkfield(md, 'fieldname', 'frontalforcings.num_breaks', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) md = checkfield(md, 'fieldname', 'frontalforcings.basin_id', 'Inf', 1, '>=', 0, '<=', md.frontalforcings.num_basins, 'size', [md.mesh.numberofelements]) md = checkfield(md, 'fieldname', 'frontalforcings.subglacial_discharge', '>=', 0, 'NaN', 1, 'Inf', 1, 'timeseries', 1) - if len(np.shape(self.const)) == 1: - self.const = np.array([self.const]) - self.trend = np.array([self.trend]) - md = checkfield(md, 'fieldname', 'frontalforcings.const', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins) - md = checkfield(md, 'fieldname', 'frontalforcings.trend', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins) + if len(np.shape(self.polynomialparams)) == 1: + self.polynomialparams = np.array([[self.polynomialparams]]) + if(nbas>1 and nbrk>=1 and nprm>1): + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm) + elif(nbas==1): + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm) + elif(nbrk==0): + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm) + elif(nprm==1): + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk],'numel',nbas*(nbrk+1)*nprm) md = checkfield(md, 'fieldname', 'frontalforcings.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) md = checkfield(md, 'fieldname', 'frontalforcings.ma_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) - md = checkfield(md, 'fieldname', 'frontalforcings.arma_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'frontalforcings.arma_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # ARMA time step cannot be finer than ISSM timestep md = checkfield(md, 'fieldname', 'frontalforcings.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ar_order]) md = checkfield(md, 'fieldname', 'frontalforcings.malag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ma_order]) + if(nbrk>0): + md = checkfield(md, 'fieldname', 'frontalforcings.datebreaks', 'NaN', 1, 'Inf', 1, 'size', [nbas,nbrk]) + elif(np.size(md.frontalforcings.datebreaks)==0 or np.all(np.isnan(md.frontalforcings.datebreaks))): + pass + else: + raise RuntimeError('md.frontalforcings.num_breaks is 0 but md.frontalforcings.datebreaks is not empty') if(np.any(np.isnan(md.frontalforcings.monthly_effects)==False)): md = checkfield(md, 'fieldname', 'frontalforcings.monthly_effects', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, 12]) @@ -91,20 +109,56 @@ def marshall(self, prefix, md, fid): # {{{ yts = md.constants.yts + nbas = md.frontalforcings.num_basins; + nprm = md.frontalforcings.num_params; + nper = md.frontalforcings.num_breaks+1; + # Scale the parameters # + polyparamsScaled = np.copy(md.frontalforcings.polynomialparams) + polyparams2dScaled = np.zeros((nbas,nper*nprm)) + if(nprm>1): + # Case 3D # + if(nbas>1 and nper>1): + for ii in range(nprm): + polyparamsScaled[:,:,ii] = polyparamsScaled[:,:,ii]*(1/yts)**ii + # Fit in 2D array # + for ii in range(nprm): + polyparams2dScaled[:,ii*nper:(ii+1)*nper] = 1*polyparamsScaled[:,:,ii] + # Case 2D and higher-order params at increasing row index # + elif(nbas==1): + for ii in range(nprm): + polyparamsScaled[ii,:] = polyparamsScaled[ii,:]*(1/yts)**ii + # Fit in row array # + for ii in range(nprm): + polyparams2dScaled[0,ii*nper:(ii+1)*nper] = 1*polyparamsScaled[ii,:] + # Case 2D and higher-order params at incrasing column index # + elif(nper==1): + for ii in range(nprm): + polyparamsScaled[:,ii] = polyparamsScaled[:,ii]*(1/yts)**ii + # 2D array is already in correct format # + polyparams2dScaled = np.copy(polyparamsScaled) + else: + # 2D array is already in correct format and no need for scaling # + polyparams2dScaled = np.copy(polyparamsScaled) if(np.any(np.isnan(md.frontalforcings.monthly_effects))): #monthly effects not provided, set to 0 meffects = np.zeros((md.frontalforcings.num_basins,12)) else: meffects = 1*md.frontalforcings.monthly_effects + if(nper==1): + dbreaks = np.zeros((nbas,1)) + else: + dbreaks = np.copy(md.frontalforcings.datebreaks) + WriteData(fid, prefix, 'name', 'md.frontalforcings.parameterization', 'data', 3, 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'num_basins', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'num_breaks', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'num_params', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'subglacial_discharge', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ar_order', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ma_order', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'arma_initialtime', 'format', 'Double', 'scale', yts) WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'arma_timestep', 'format', 'Double', 'scale', yts) WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'basin_id', 'data', self.basin_id - 1, 'name', 'md.frontalforcings.basin_id', 'format', 'IntMat', 'mattype', 2) # 0-indexed - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'const', 'format', 'DoubleMat', 'name', 'md.frontalforcings.const') - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'trend', 'format', 'DoubleMat', 'name', 'md.frontalforcings.trend', 'scale', 1 / yts, 'yts', yts) + WriteData(fid, prefix, 'data', polyparams2dScaled, 'name', 'md.frontalforcings.polynomialparams', 'format', 'DoubleMat') WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.frontalforcings.arlag_coefs', 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'malag_coefs', 'format', 'DoubleMat', 'name', 'md.frontalforcings.malag_coefs', 'yts', yts) + WriteData(fid, prefix, 'data', dbreaks, 'name', 'md.frontalforcings.datebreaks', 'format', 'DoubleMat','scale',yts) WriteData(fid, prefix, 'data', meffects, 'name', 'md.frontalforcings.monthly_effects', 'format', 'DoubleMat') # }}} Index: /issm/trunk-jpl/src/m/classes/linearbasalforcingsarma.m =================================================================== --- /issm/trunk-jpl/src/m/classes/linearbasalforcingsarma.m (revision 27317) +++ /issm/trunk-jpl/src/m/classes/linearbasalforcingsarma.m (revision 27318) @@ -8,9 +8,10 @@ properties (SetAccess=public) num_basins = 0; - const = NaN; - trend = NaN; + num_params = 0; + num_breaks = 0; + polynomialparams = NaN; + datebreaks = NaN; ar_order = 0; ma_order = 0; - arma_initialtime = 0; arma_timestep = 0; arlag_coefs = NaN; @@ -59,8 +60,4 @@ disp(' basalforcings.ma_order (order of moving-average model) not specified: order of moving-average model set to 0'); end - if (self.arma_initialtime==0) - self.ar_initialtime = md.timestepping.start_time; %ARMA model has no prescribed initial time - disp(' basalforcings.arma_initialtime (initial time in the ARMA model parameterization) not specified: set to md.timestepping.start_time'); - end if (self.arma_timestep==0) self.arma_timestep = md.timestepping.time_step; %ARMA model has no prescribed time step @@ -80,5 +77,10 @@ if ismember('MasstransportAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.ismasstransport==0), + nbas = md.basalforcings.num_basins; + nprm = md.basalforcings.num_params; + nbrk = md.basalforcings.num_breaks; md = checkfield(md,'fieldname','basalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0); + md = checkfield(md,'fieldname','basalforcings.num_params','numel',1,'NaN',1,'Inf',1,'>',0); + md = checkfield(md,'fieldname','basalforcings.num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1); md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); @@ -86,12 +88,25 @@ md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','NaN',1,'Inf',1,'<=',0,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); md = checkfield(md,'fieldname','basalforcings.basin_id','Inf',1,'>=',0,'<=',md.basalforcings.num_basins,'size',[md.mesh.numberofelements,1]); - md = checkfield(md,'fieldname','basalforcings.const','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); - md = checkfield(md,'fieldname','basalforcings.trend','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); + if(nbas>1 && nbrk>=1 && nprm>1) + md = checkfield(md,'fieldname','basalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm); + elseif(nbas==1) + md = checkfield(md,'fieldname','basalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm); + elseif(nbrk==0) + md = checkfield(md,'fieldname','basalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm); + elseif(nprm==1) + md = checkfield(md,'fieldname','basalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1],'numel',nbas*(nbrk+1)*nprm); + end md = checkfield(md,'fieldname','basalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','basalforcings.ma_order','numel',1,'NaN',1,'Inf',1,'>=',0); - md = checkfield(md,'fieldname','basalforcings.arma_initialtime','numel',1,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','basalforcings.arma_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %moving-average time step cannot be finer than ISSM timestep md = checkfield(md,'fieldname','basalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.basalforcings.num_basins,md.basalforcings.ar_order]); md = checkfield(md,'fieldname','basalforcings.malag_coefs','NaN',1,'Inf',1,'size',[md.basalforcings.num_basins,md.basalforcings.ma_order]); + if(nbrk>0) + md = checkfield(md,'fieldname','basalforcings.datebreaks','NaN',1,'Inf',1,'size',[nbas,nbrk]); + elseif(numel(md.basalforcings.datebreaks)==0 || all(isnan(md.basalforcings.datebreaks))) + ; + else + error('md.basalforcings.num_breaks is 0 but md.basalforcings.datebreaks is not empty'); + end end if ismember('BalancethicknessAnalysis',analyses), @@ -111,9 +126,11 @@ fielddisplay(self,'num_basins','number of different basins [unitless]'); fielddisplay(self,'basin_id','basin number assigned to each element [unitless]'); - fielddisplay(self,'const','basin-specific constant values [m/yr]'); - fielddisplay(self,'trend','basin-specific trend values [m yr^(-2)]'); - fielddisplay(self,'ar_order','order of the autoregressive model [unitless]'); + fielddisplay(self,'num_breaks','number of different breakpoints in the piecewise-polynomial (separating num_breaks+1 periods)'); + fielddisplay(self,'num_params','number of different parameters in the piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)'); + fielddisplay(self,'polynomialparams','coefficients for the polynomial (const,trend,quadratic,etc.),dim1 for basins,dim2 for periods,dim3 for orders'); + disp(sprintf('%51s ex: polyparams=cat(3,intercepts,trendlinearcoefs,trendquadraticcoefs)',' ')); + fielddisplay(self,'datebreaks','dates at which the breakpoints in the piecewise polynomial occur (1 row per basin) [yr]'); + fielddisplay(self,'ar_order','order of the autoregressive model [unitless]'); fielddisplay(self,'ar_order','order of the moving-average model [unitless]'); - fielddisplay(self,'arma_initialtime','initial time assumed in the ARMA model parameterization [yr]'); fielddisplay(self,'arma_timestep','time resolution of the ARMA model [yr]'); fielddisplay(self,'arlag_coefs','basin-specific vectors of AR lag coefficients [unitless]'); @@ -129,4 +146,49 @@ yts=md.constants.yts; + nbas = md.basalforcings.num_basins; + nprm = md.basalforcings.num_params; + nper = md.basalforcings.num_breaks+1; + % Scale the parameters % + polyparamsScaled = md.basalforcings.polynomialparams; + polyparams2dScaled = zeros(nbas,nper*nprm); + if(nprm>1) + % Case 3D % + if(nbas>1 && nper>1) + for(ii=[1:nprm]) + polyparamsScaled(:,:,ii) = polyparamsScaled(:,:,ii)*((1/yts)^(ii)); + end + % Fit in 2D array % + for(ii=[1:nprm]) + jj = 1+(ii-1)*nper; + polyparams2dScaled(:,jj:jj+nper-1) = polyparamsScaled(:,:,ii); + end + % Case 2D and higher-order params at increasing row index % + elseif(nbas==1) + for(ii=[1:nprm]) + polyparamsScaled(ii,:) = polyparamsScaled(ii,:)*((1/yts)^(ii)); + end + % Fit in row array % + for(ii=[1:nprm]) + jj = 1+(ii-1)*nper; + polyparams2dScaled(1,jj:jj+nper-1) = polyparamsScaled(ii,:); + end + % Case 2D and higher-order params at incrasing column index % + elseif(nper==1) + for(ii=[1:nprm]) + polyparamsScaled(:,ii) = polyparamsScaled(:,ii)*((1/yts)^(ii)); + end + % 2D array is already in correct format % + polyparams2dScaled = polyparamsScaled; + end + else + polyparamsScaled = polyparamsScaled*(1/yts); + % 2D array is already in correct format % + polyparams2dScaled = polyparamsScaled; + end + if(nper==1) %a single period (no break date) + dbreaks = zeros(nbas,1); %dummy + else + dbreaks = md.basalforcings.datebreaks; + end WriteData(fid,prefix,'name','md.basalforcings.model','data',9,'format','Integer'); @@ -134,13 +196,14 @@ WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','name','md.basalforcings.geothermalflux','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','num_basins','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','num_params','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','num_breaks','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','ar_order','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','ma_order','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','arma_initialtime','format','Double','scale',yts); WriteData(fid,prefix,'object',self,'fieldname','arma_timestep','format','Double','scale',yts); WriteData(fid,prefix,'object',self,'fieldname','basin_id','data',self.basin_id-1,'name','md.basalforcings.basin_id','format','IntMat','mattype',2); %0-indexed - WriteData(fid,prefix,'object',self,'fieldname','const','format','DoubleMat','name','md.basalforcings.const','scale',1./yts,'yts',yts); - WriteData(fid,prefix,'object',self,'fieldname','trend','format','DoubleMat','name','md.basalforcings.trend','scale',1./(yts^2),'yts',yts); WriteData(fid,prefix,'object',self,'fieldname','arlag_coefs','format','DoubleMat','name','md.basalforcings.arlag_coefs','yts',yts); WriteData(fid,prefix,'object',self,'fieldname','malag_coefs','format','DoubleMat','name','md.basalforcings.malag_coefs','yts',yts); + WriteData(fid,prefix,'data',polyparams2dScaled,'name','md.basalforcings.polynomialparams','format','DoubleMat'); + WriteData(fid,prefix,'data',dbreaks,'name','md.basalforcings.datebreaks','format','DoubleMat','scale',yts); WriteData(fid,prefix,'object',self,'fieldname','deepwater_elevation','format','DoubleMat','name','md.basalforcings.deepwater_elevation'); WriteData(fid,prefix,'object',self,'fieldname','upperwater_melting_rate','format','DoubleMat','name','md.basalforcings.upperwater_melting_rate','scale',1./yts); Index: /issm/trunk-jpl/src/m/classes/linearbasalforcingsarma.py =================================================================== --- /issm/trunk-jpl/src/m/classes/linearbasalforcingsarma.py (revision 27317) +++ /issm/trunk-jpl/src/m/classes/linearbasalforcingsarma.py (revision 27318) @@ -15,9 +15,10 @@ def __init__(self, *args): # {{{ self.num_basins = 0 - self.const = np.nan - self.trend = np.nan + self.num_params = 0 + self.num_breaks = 0 + self.polynomialparams = np.nan + self.datebreaks = np.nan self.ar_order = 0 self.ma_order = 0 - self.arma_initialtime = 0 self.arma_timestep = 0 self.arlag_coefs = np.nan @@ -42,9 +43,10 @@ s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]')) s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'const', 'basin-specific constant values [m ice eq./yr]')) - s += '{}\n'.format(fielddisplay(self, 'trend', 'basin-specific trend values [m ice eq. yr^(-2)]')) + s += '{}\n'.format(fielddisplay(self, 'num_breaks', 'number of different breakpoints in the piecewise-polynomial (separating num_breaks+1 periods)')) + s += '{}\n'.format(fielddisplay(self, 'num_params', 'number of different parameters in the piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)')) + s += '{}\n'.format(fielddisplay(self, 'polynomialparams', 'coefficients for the polynomial (const,trend,quadratic,etc.),dim1 for basins,dim2 for periods,dim3 for orders, ex: polyparams=cat(num_params,intercepts,trendlinearcoefs,trendquadraticcoefs)')) + s += '{}\n'.format(fielddisplay(self, 'datebreaks', 'dates at which the breakpoints in the piecewise polynomial occur (1 row per basin) [yr]')) s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]')) s += '{}\n'.format(fielddisplay(self, 'ma_order', 'order of the moving-average model [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'arma_initialtime', 'initial time assumed in the ARMA model parameterization [yr]')) s += '{}\n'.format(fielddisplay(self, 'arma_timestep', 'time resolution of the ARMA model [yr]')) s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of AR lag coefficients [unitless]')) @@ -79,7 +81,4 @@ self.arlag_coefs = np.zeros((self.num_basins, self.ar_order)) # Autoregression coefficients all set to 0 print(' basalforcings.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0') - if self.arma_initialtime == 0: - self.arma_initialtime = md.timestepping.start_time # ARMA model has no prescribed initial time - print(' basalforcings.arma_initialtime (initial time in the ARMA model parameterization) not specified: set to md.timestepping.start_time') if self.arma_timestep == 0: self.arma_timestep = md.timestepping.time_step # ARMA model has no prescribed time step @@ -96,6 +95,11 @@ def checkconsistency(self, md, solution, analyses): # {{{ if 'MasstransportAnalysis' in analyses: + nbas = md.basalforcings.num_basins; + nprm = md.basalforcings.num_params; + nbrk = md.basalforcings.num_breaks; md = checkfield(md, 'fieldname', 'basalforcings.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'basalforcings.num_params', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) + md = checkfield(md, 'fieldname', 'basalforcings.num_breaks', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) if len(np.shape(self.deepwater_elevation)) == 1: @@ -103,4 +107,14 @@ self.upperwater_elevation = np.array([self.upperwater_elevation]) self.upperwater_melting_rate = np.array([self.upperwater_melting_rate]) + if len(np.shape(self.polynomialparams)) == 1: + self.polynomialparams = np.array([[self.polynomialparams]]) + if(nbas>1 and nbrk>=1 and nprm>1): + md = checkfield(md,'fieldname','basalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm) + elif(nbas==1): + md = checkfield(md,'fieldname','basalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm) + elif(nbrk==0): + md = checkfield(md,'fieldname','basalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm) + elif(nprm==1): + md = checkfield(md,'fieldname','basalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk],'numel',nbas*(nbrk+1)*nprm) md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', 'NaN', 1, 'Inf', 1, '<=', 0, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) @@ -108,16 +122,17 @@ md = checkfield(md, 'fieldname', 'basalforcings.basin_id', 'Inf', 1, '>=', 0, '<=', md.basalforcings.num_basins, 'size', [md.mesh.numberofelements]) - if len(np.shape(self.const)) == 1: - self.const = np.array([self.const]) - self.trend = np.array([self.trend]) - - md = checkfield(md, 'fieldname', 'basalforcings.const', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) # Scheme fails if passed as column vector - md = checkfield(md, 'fieldname', 'basalforcings.trend', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) # Scheme fails if passed as column vector; NOTE: As opposed to MATLAB implementation, pass list + md = checkfield(md, 'fieldname', 'basalforcings.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) md = checkfield(md, 'fieldname', 'basalforcings.ma_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) - md = checkfield(md, 'fieldname', 'basalforcings.arma_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'basalforcings.arma_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep md = checkfield(md, 'fieldname', 'basalforcings.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.basalforcings.num_basins, md.basalforcings.ar_order]) md = checkfield(md, 'fieldname', 'basalforcings.malag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.basalforcings.num_basins, md.basalforcings.ma_order]) + if(nbrk>0): + md = checkfield(md, 'fieldname', 'basalforcings.datebreaks', 'NaN', 1, 'Inf', 1, 'size', [nbas,nbrk]) + elif(np.size(md.basalforcings.datebreaks)==0 or np.all(np.isnan(md.basalforcings.datebreaks))): + pass + else: + raise RuntimeError('md.basalforcings.num_breaks is 0 but md.basalforcings.datebreaks is not empty') + if 'BalancethicknessAnalysis' in analyses: raise Exception('not implemented yet!') @@ -130,4 +145,39 @@ def marshall(self, prefix, md, fid): # {{{ yts = md.constants.yts + nbas = md.basalforcings.num_basins; + nprm = md.basalforcings.num_params; + nper = md.basalforcings.num_breaks+1; + # Scale the parameters # + polyparamsScaled = np.copy(md.basalforcings.polynomialparams) + polyparams2dScaled = np.zeros((nbas,nper*nprm)) + if(nprm>1): + # Case 3D # + if(nbas>1 and nper>1): + for ii in range(nprm): + polyparamsScaled[:,:,ii] = polyparamsScaled[:,:,ii]*(1/yts)**(ii+1) + # Fit in 2D array # + for ii in range(nprm): + polyparams2dScaled[:,ii*nper:(ii+1)*nper] = 1*polyparamsScaled[:,:,ii] + # Case 2D and higher-order params at increasing row index # + elif(nbas==1): + for ii in range(nprm): + polyparamsScaled[ii,:] = polyparamsScaled[ii,:]*(1/yts)**(ii+1) + # Fit in row array # + for ii in range(nprm): + polyparams2dScaled[0,ii*nper:(ii+1)*nper] = 1*polyparamsScaled[ii,:] + # Case 2D and higher-order params at incrasing column index # + elif(nper==1): + for ii in range(nprm): + polyparamsScaled[:,ii] = polyparamsScaled[:,ii]*(1/yts)**(ii+1) + # 2D array is already in correct format # + polyparams2dScaled = np.copy(polyparamsScaled) + else: + polyparamsScaled = polyparamsScaled*(1/yts) + # 2D array is already in correct format # + polyparams2dScaled = np.copy(polyparamsScaled) + if(nper==1): + dbreaks = np.zeros((nbas,1)) + else: + dbreaks = np.copy(md.basalforcings.datebreaks) WriteData(fid, prefix, 'name', 'md.basalforcings.model', 'data', 9, 'format', 'Integer') @@ -135,13 +185,14 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'geothermalflux', 'name', 'md.basalforcings.geothermalflux', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'num_basins', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'num_params', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'num_breaks', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'ar_order', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'ma_order', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'arma_initialtime', 'format', 'Double', 'scale', yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'arma_timestep', 'format', 'Double', 'scale', yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'basin_id', 'data', self.basin_id - 1, 'name', 'md.basalforcings.basin_id', 'format', 'IntMat', 'mattype', 2) # 0-indexed - WriteData(fid, prefix, 'object', self, 'fieldname', 'const', 'format', 'DoubleMat', 'name', 'md.basalforcings.const', 'scale', 1 / yts, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'trend', 'format', 'DoubleMat', 'name', 'md.basalforcings.trend', 'scale', 1 / (yts ** 2), 'yts', yts) + WriteData(fid, prefix, 'data', polyparams2dScaled, 'name', 'md.basalforcings.polynomialparams', 'format', 'DoubleMat') WriteData(fid, prefix, 'object', self, 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.basalforcings.arlag_coefs', 'yts', yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'malag_coefs', 'format', 'DoubleMat', 'name', 'md.basalforcings.malag_coefs', 'yts', yts) + WriteData(fid, prefix, 'data', dbreaks, 'name', 'md.basalforcings.datebreaks', 'format', 'DoubleMat','scale',yts) WriteData(fid, prefix, 'object', self, 'fieldname', 'deepwater_elevation', 'format', 'DoubleMat', 'name', 'md.basalforcings.deepwater_elevation') WriteData(fid, prefix, 'object', self, 'fieldname', 'upperwater_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.upperwater_melting_rate', 'scale', 1 / yts, 'yts', yts) Index: /issm/trunk-jpl/test/NightlyRun/test257.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.m (revision 27317) +++ /issm/trunk-jpl/test/NightlyRun/test257.m (revision 27318) @@ -33,13 +33,21 @@ %SMB parameters +numparams = 2; +numbreaks = 1; +intercept = [0.5,1.0;1.0,0.6;,2.0,3.0]; %intercept values of SMB in basins [m ice eq./yr] +trendlin = [0.0,0.0;0.01,0.001;-0.01,0]; %trend values of SMB in basins [m ice eq./yr^2] +polynomialparams = cat(3,intercept,trendlin); +datebreaks = [3;3;3]; + md.timestepping.start_time = 0; md.timestepping.time_step = 1; -md.timestepping.final_time = 5; +md.timestepping.final_time = 7; md.smb = SMBarma(); md.smb.num_basins = 3; %number of basins md.smb.basin_id = idbasin; %prescribe basin ID number to elements -md.smb.const = [0.5,1.2,1.5]; %intercept values of SMB in basins [m ice eq./yr] -md.smb.trend = [0.0,0.01,-0.01]; %trend values of SMB in basins [m ice eq./yr^2] -md.smb.arma_initialtime = md.timestepping.start_time; +md.smb.num_params = numparams; %number of parameters in the polynomial +md.smb.num_breaks = numbreaks; %number of breakpoints +md.smb.polynomialparams = polynomialparams; +md.smb.datebreaks = datebreaks; md.smb.ar_order = 4; md.smb.ma_order = 1; @@ -63,5 +71,5 @@ field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,... 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13... - 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; + 1e-12,1e-12,1e-12,1e-12,1e-12,1e-12}; field_values={... (md.results.TransientSolution(1).Vx),... @@ -77,9 +85,9 @@ (md.results.TransientSolution(2).IceVolume),... (md.results.TransientSolution(2).SmbMassBalance),... - (md.results.TransientSolution(3).Vx),... - (md.results.TransientSolution(3).Vy),... - (md.results.TransientSolution(3).Vel),... - (md.results.TransientSolution(3).Thickness),... - (md.results.TransientSolution(3).IceVolume),... - (md.results.TransientSolution(3).SmbMassBalance),... + (md.results.TransientSolution(7).Vx),... + (md.results.TransientSolution(7).Vy),... + (md.results.TransientSolution(7).Vel),... + (md.results.TransientSolution(7).Thickness),... + (md.results.TransientSolution(7).IceVolume),... + (md.results.TransientSolution(7).SmbMassBalance),... }; Index: /issm/trunk-jpl/test/NightlyRun/test257.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.py (revision 27317) +++ /issm/trunk-jpl/test/NightlyRun/test257.py (revision 27318) @@ -38,13 +38,21 @@ # SMB parameters +numparams = 2 +numbreaks = 1 +intercept = np.array([[0.5,1.0],[1.0,0.6],[2.0,3.0]]) +trendlin = np.array([[0.0,0.0],[0.01,0.001],[-0.01,0]]) +polynomialparams = np.transpose(np.stack((intercept,trendlin)),(1,2,0)) +datebreaks = np.array([[3],[3],[3]]); + md.timestepping.start_time = 0 md.timestepping.time_step = 1 -md.timestepping.final_time = 5 +md.timestepping.final_time = 8 md.smb = SMBarma() md.smb.num_basins = 3 # number of basins md.smb.basin_id = idbasin # prescribe basin ID number to elements; -md.smb.const = np.array([[0.5, 1.2, 1.5]]) # intercept values of SMB in basins [m ice eq./yr] -md.smb.trend = np.array([[0.0, 0.01, -0.01]]) # trend values of SMB in basins [m ice eq./yr^2] -md.smb.arma_initialtime = md.timestepping.start_time +md.smb.num_params = numparams +md.smb.num_breaks = numbreaks +md.smb.polynomialparams = polynomialparams +md.smb.datebreaks = datebreaks md.smb.ar_order = 4 md.smb.ma_order = 1 @@ -73,5 +81,5 @@ 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, - 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13 + 1e-12, 1e-12, 1e-12, 1e-12, 1e-12, 1e-12 ] field_values = [ @@ -88,9 +96,9 @@ md.results.TransientSolution[1].IceVolume, md.results.TransientSolution[1].SmbMassBalance, - md.results.TransientSolution[2].Vx, - md.results.TransientSolution[2].Vy, - md.results.TransientSolution[2].Vel, - md.results.TransientSolution[2].Thickness, - md.results.TransientSolution[2].IceVolume, - md.results.TransientSolution[2].SmbMassBalance + md.results.TransientSolution[6].Vx, + md.results.TransientSolution[6].Vy, + md.results.TransientSolution[6].Vel, + md.results.TransientSolution[6].Thickness, + md.results.TransientSolution[6].IceVolume, + md.results.TransientSolution[6].SmbMassBalance ] Index: /issm/trunk-jpl/test/NightlyRun/test543.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.m (revision 27317) +++ /issm/trunk-jpl/test/NightlyRun/test543.m (revision 27318) @@ -44,12 +44,22 @@ md.levelset.spclevelset = NaN(md.mesh.numberofvertices,1); md.levelset.migration_max = 10.0; %avoid fast advance/retreat of the front -%Frontal forcing parameters +%%% Frontal forcing parameters %%% md.frontalforcings=frontalforcingsrignotarma(); +% Polynomial params % +numparams = 3; +numbreaks = 2; +intercept = [2.5,2.0,0.1;0.5,0.5,1.5]; +trendlin = [-0.5,-0.2,0.1;0,0,0]; +trendquad = [0,0.0,0;0.1,0.1,0.1]; +datebreaks = [4,7;4,7]; +polynomialparams = cat(numparams,intercept,trendlin,trendquad); + md.frontalforcings.num_basins = nb_tf; md.frontalforcings.basin_id = idb_tf; +md.frontalforcings.num_params = numparams; %number of parameters in the polynomial +md.frontalforcings.num_breaks = numbreaks; %number of breakpoints md.frontalforcings.subglacial_discharge = 0.01*ones(md.mesh.numberofvertices,1); -md.frontalforcings.const = [0.05,0.01]; %intercept values of TF in basins [C] -md.frontalforcings.trend = [0.0001,0.00001]; %trend values of TF in basins [C/yr] -md.frontalforcings.arma_initialtime = md.timestepping.start_time; %initial time in the AR model parameterization [yr] +md.frontalforcings.polynomialparams = polynomialparams; +md.frontalforcings.datebreaks = datebreaks; md.frontalforcings.ar_order = 4; md.frontalforcings.ma_order = 2; Index: /issm/trunk-jpl/test/NightlyRun/test543.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.py (revision 27317) +++ /issm/trunk-jpl/test/NightlyRun/test543.py (revision 27318) @@ -49,12 +49,22 @@ md.levelset.spclevelset = np.full((md.mesh.numberofvertices,), np.nan) md.levelset.migration_max = 10.0 -#Frontal forcing parameters +### Frontal forcing parameters ### md.frontalforcings = frontalforcingsrignotarma() +# Polynomial params # +numparams = 3; +numbreaks = 2; +intercept = np.array([[2.5,2.0,0.1],[0.5,0.5,1.5]]) +trendlin = np.array([[-0.5,-0.2,0.1],[0,0,0]]) +trendquad = np.array([[0,0.0,0],[0.1,0.1,0.1]]) +datebreaks = np.array([[4,7],[4,7]]) +polynomialparams = np.transpose(np.stack((intercept,trendlin,trendquad)),(1,2,0)) + md.frontalforcings.num_basins = nb_tf md.frontalforcings.basin_id = idb_tf md.frontalforcings.subglacial_discharge = 0.01 * np.ones((md.mesh.numberofvertices,)) -md.frontalforcings.const = np.array([[0.05, 0.01]]) # intercept values of TF in basins [C] -md.frontalforcings.trend = np.array([[0.0001, 0.00001]]) # trend values of TF in basins [C/yr] -md.frontalforcings.arma_initialtime = md.timestepping.start_time # initial time in the AR model parameterization [yr] +md.frontalforcings.num_params = numparams #number of parameters in the polynomial +md.frontalforcings.num_breaks = numbreaks #number of breakpoints +md.frontalforcings.polynomialparams = polynomialparams +md.frontalforcings.datebreaks = datebreaks md.frontalforcings.ar_order = 4 md.frontalforcings.ma_order = 2 Index: /issm/trunk-jpl/test/NightlyRun/test544.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test544.m (revision 27317) +++ /issm/trunk-jpl/test/NightlyRun/test544.m (revision 27318) @@ -24,10 +24,16 @@ %SMB +numparams = 1; +numbreaks = 0; +intercept = [0.5;0.01]; +polynomialparams = intercept; +datebreaks = NaN; md.smb = SMBarma(); md.smb.num_basins = nb_bas; %number of basins md.smb.basin_id = idb; %prescribe basin ID number to elements -md.smb.const = [0.5,1.2]; %intercept values of SMB in basins [m ice eq./yr] -md.smb.trend = [0.0,0.01]; %trend values of SMB in basins [m ice eq./yr^2] -md.smb.arma_initialtime = md.timestepping.start_time; +md.smb.num_params = numparams; %number of parameters in the polynomial +md.smb.num_breaks = numbreaks; %number of breakpoints +md.smb.polynomialparams = polynomialparams; +md.smb.datebreaks = datebreaks; md.smb.ar_order = 4; md.smb.ma_order = 4; @@ -44,10 +50,17 @@ % Basal forcing implementation -md.basalforcings = linearbasalforcingsarma(); +numparams = 2; +numbreaks = 1; +intercept = [3.0,4.0;1.0,0.5]; +trendlin = [0.0,0.1;0,0]; +polynomialparams = cat(3,intercept,trendlin); +datebreaks = [6;7]; +md.basalforcings = linearbasalforcingsarma(); md.basalforcings.num_basins = nb_bas; %number of basins md.basalforcings.basin_id = idb; %prescribe basin ID number to elements -md.basalforcings.const = [1.0,2.50]; %intercept values of DeepwaterMelt in basins [m/yr] -md.basalforcings.trend = [0.2,0.01]; %trend values of DeepwaterMelt in basins [m/yr^2] -md.basalforcings.arma_initialtime = md.timestepping.start_time; +md.basalforcings.polynomialparams = polynomialparams; +md.basalforcings.datebreaks = datebreaks; +md.basalforcings.num_params = numparams; %number of parameters in the polynomial +md.basalforcings.num_breaks = numbreaks; %number of breakpoints md.basalforcings.ar_order = 1; md.basalforcings.ma_order = 1; Index: /issm/trunk-jpl/test/NightlyRun/test544.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test544.py (revision 27317) +++ /issm/trunk-jpl/test/NightlyRun/test544.py (revision 27318) @@ -34,10 +34,16 @@ #SMB +numparams = 1; +numbreaks = 0; +intercept = np.array([[0.5],[0.01]]) +polynomialparams = np.copy(intercept) +datebreaks = np.nan md.smb = SMBarma() md.smb.num_basins = nb_bas # number of basins md.smb.basin_id = idb # prescribe basin ID number to elements; -md.smb.const = np.array([[0.5, 1.2]]) # intercept values of SMB in basins [m ice eq./yr] -md.smb.trend = np.array([[0.0, 0.01]]) # trend values of SMB in basins [m ice eq./yr^2] -md.smb.arma_initialtime = md.timestepping.start_time +md.smb.num_params = 1*numparams +md.smb.num_breaks = 1*numbreaks +md.smb.polynomialparams = 1*polynomialparams +md.smb.datebreaks = 1*datebreaks md.smb.ar_order = 4 md.smb.ma_order = 4 @@ -54,4 +60,11 @@ #Basal forcing implementation +numparams = 2 +numbreaks = 1 +intercept = np.array([[3.0,4.0],[1.0,0.5]]) +trendlin = np.array([[0.0,0.1],[0,0]]) +polynomialparams = np.transpose(np.stack((intercept,trendlin)),(1,2,0)) +datebreaks = np.array([[6],[7]]) + md.basalforcings = linearbasalforcingsarma() md.basalforcings.num_basins = nb_bas @@ -62,4 +75,8 @@ md.basalforcings.ar_order = 1 md.basalforcings.ma_order = 1 +md.basalforcings.polynomialparams = 1*polynomialparams; +md.basalforcings.datebreaks = 1*datebreaks; +md.basalforcings.num_params = 1*numparams +md.basalforcings.num_breaks = 1*numbreaks md.basalforcings.arma_timestep = 1.0 # timestep of the ARMA model [yr] md.basalforcings.arlag_coefs = np.array([[0.0], [0.1]]) # autoregressive parameters