Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27245) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27246) @@ -243,12 +243,12 @@ parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_initialtime",FrontalForcingsAutoregressionInitialTimeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_timestep",FrontalForcingsAutoregressionTimestepEnum)); - iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.beta0"); - parameters->AddObject(new DoubleVecParam(FrontalForcingsBeta0Enum,transparam,N)); + iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.const"); + parameters->AddObject(new DoubleVecParam(FrontalForcingsautoregressionconstEnum,transparam,N)); xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.beta1"); - parameters->AddObject(new DoubleVecParam(FrontalForcingsBeta1Enum,transparam,N)); + iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.trend"); + parameters->AddObject(new DoubleVecParam(FrontalForcingsautoregressiontrendEnum,transparam,N)); xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.phi"); - parameters->AddObject(new DoubleMatParam(FrontalForcingsPhiEnum,transparam,M,N)); + iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.arlag_coefs"); + parameters->AddObject(new DoubleMatParam(FrontalForcingsarlagcoefsEnum,transparam,M,N)); xDelete(transparam); /*Do not break here, generic FrontalForcingsRignot parameters still to be retrieved*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27245) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27246) @@ -61,52 +61,11 @@ return false; }/*}}}*/ -void Element::AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type){/*{{{*/ - -const int numvertices = this->GetNumberOfVertices(); - int basinid; - IssmDouble tspin,beta0_basin,beta1_basin,noisespin_basin; - IssmDouble* phi_basin = xNew(arorder); - IssmDouble* varspin = xNewZeroInit(numvertices*arorder); - - /*Get Basin ID and Basin coefficients*/ - if(enum_type==SMBautoregressionEnum) this->GetInputValue(&basinid,SmbBasinsIdEnum); - if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum); - if(enum_type==BasalforcingsDeepwaterMeltingRateAutoregressionEnum) this->GetInputValue(&basinid,BasalforcingsLinearBasinIdEnum); - for(int i=0;i(numvertices*arorder); - for(int i=0;i(oldvarspin); - } - if(enum_type==SMBautoregressionEnum) this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder); - if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder); - if(enum_type==BasalforcingsDeepwaterMeltingRateAutoregressionEnum) this->inputs->SetArrayInput(BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder); - - /*Cleanup and return*/ - xDelete(varspin); - xDelete(phi_basin); -}/*}}}*/ -void Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type){/*{{{*/ +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 beta0_basin,beta1_basin,noiseterm; - IssmDouble* phi_basin = xNew(arorder); - IssmDouble* varlist = xNew(numvertices); + 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; @@ -134,6 +93,24 @@ } + /*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*/ - this->GetInputValue(&basinid,basinenum_type); if(isfieldstochastic){ noiseterm_input = this->GetInput(noiseenum_type); @@ -145,9 +122,4 @@ this->inputs->GetArray(arenum_type,this->lid,&valuesautoregression,&M); - /*Get basin coefficients*/ - for(int ii=0;ii(phi_basin); + xDelete(lagcoefs_basin); xDelete(varlist); xDelete(valuesautoregression); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27245) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27246) @@ -68,6 +68,5 @@ /*bool AnyActive(void);*/ bool AnyFSet(void); - void AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type); - void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,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 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 27245) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 27246) @@ -210,32 +210,4 @@ } /*}}}*/ -void AutoregressionLinearFloatingiceMeltingRateInitx(FemModel* femmodel){/*{{{*/ - - /*Initialization step of AutoregressionLinearFloatingiceMeltingRatex*/ - int M,N,Nphi,arorder,numbasins,my_rank; - IssmDouble starttime,tstep_ar,tinit_ar; - femmodel->parameters->FindParam(&numbasins,BasalforcingsLinearNumBasinsEnum); - femmodel->parameters->FindParam(&arorder,BasalforcingsAutoregressiveOrderEnum); - IssmDouble* beta0 = NULL; - IssmDouble* beta1 = NULL; - IssmDouble* phi = NULL; - femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); - femmodel->parameters->FindParam(&tstep_ar,BasalforcingsAutoregressionTimestepEnum); - femmodel->parameters->FindParam(&tinit_ar,BasalforcingsAutoregressionInitialTimeEnum); - femmodel->parameters->FindParam(&beta0,&M,BasalforcingsBeta0Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&beta1,&M,BasalforcingsBeta1Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&phi,&M,&Nphi,BasalforcingsPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); - - /*AR model spin-up with 0 noise to initialize BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum (688 = log(0.001)/log(0.99): decaying time of inluence of phi[0]=0.99 to 0.001 of beta_0*/ - int nspin = 688; - for(Object* &object:femmodel->elements->objects){ - Element* element = xDynamicCast(object); //generate element object - element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,BasalforcingsDeepwaterMeltingRateAutoregressionEnum); - } - /*Cleanup*/ - xDelete(beta0); - xDelete(beta1); - xDelete(phi); -}/*}}}*/ void AutoregressionLinearFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/ @@ -247,6 +219,4 @@ femmodel->parameters->FindParam(&tstep_ar,BasalforcingsAutoregressionTimestepEnum); - /*Initialize module at first time step*/ - if(time<=starttime+dt){AutoregressionLinearFloatingiceMeltingRateInitx(femmodel);} /*Determine if this is a time step for the AR model*/ bool isstepforar = false; @@ -261,11 +231,11 @@ bool isstochastic; bool isdeepmeltingstochastic = false; - int M,N,Nphi,arorder,numbasins,my_rank; + int M,N,Nlags,arorder,numbasins,my_rank; femmodel->parameters->FindParam(&numbasins,BasalforcingsLinearNumBasinsEnum); femmodel->parameters->FindParam(&arorder,BasalforcingsAutoregressiveOrderEnum); IssmDouble tinit_ar; - IssmDouble* beta0 = NULL; - IssmDouble* beta1 = NULL; - IssmDouble* phi = NULL; + IssmDouble* termconstant = NULL; + IssmDouble* trend = NULL; + IssmDouble* lagcoefs = NULL; IssmDouble* deepwaterel = NULL; IssmDouble* upperwaterel = NULL; @@ -275,7 +245,7 @@ /*Get autoregressive parameters*/ femmodel->parameters->FindParam(&tinit_ar,BasalforcingsAutoregressionInitialTimeEnum); - femmodel->parameters->FindParam(&beta0,&M,BasalforcingsBeta0Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&beta1,&M,BasalforcingsBeta1Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&phi,&M,&Nphi,BasalforcingsPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); + femmodel->parameters->FindParam(&termconstant,&M,BasalforcingsautoregressionconstEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&trend,&M,BasalforcingsautoregressiontrendEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&lagcoefs,&M,&Nlags,BasalforcingsarlagcoefsEnum); _assert_(M==numbasins); _assert_(Nlags==arorder); /*Get basin-specific parameters*/ @@ -303,12 +273,12 @@ Element* element = xDynamicCast(object); /*Compute autoregression*/ - element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,isdeepmeltingstochastic,BasalforcingsDeepwaterMeltingRateAutoregressionEnum); + element->Autoregression(isstepforar,arorder,telapsed_ar,tstep_ar,termconstant,trend,lagcoefs,isdeepmeltingstochastic,BasalforcingsDeepwaterMeltingRateAutoregressionEnum); element->BasinLinearFloatingiceMeltingRate(deepwaterel,upperwatermelt,upperwaterel,perturbation); } /*Cleanup*/ - xDelete(beta0); - xDelete(beta1); - xDelete(phi); + xDelete(termconstant); + xDelete(trend); + xDelete(lagcoefs); xDelete(deepwaterel); xDelete(upperwaterel); Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h (revision 27245) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h (revision 27246) @@ -16,5 +16,4 @@ void FloatingiceMeltingRateIsmip6x(FemModel* femmodel); void BeckmannGoosseFloatingiceMeltingRatex(FemModel* femmodel); -void AutoregressionLinearFloatingiceMeltingRateInitx(FemModel* femmodel); void AutoregressionLinearFloatingiceMeltingRatex(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 27245) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 27246) @@ -29,33 +29,4 @@ } }/*}}}*/ -void ThermalforcingautoregressionInitx(FemModel* femmodel){/*{{{*/ - - /*Initialization step of Thermalforcingautoregressionx*/ - int M,N,Nphi,arorder,numbasins,my_rank; - IssmDouble starttime,tstep_ar,tinit_ar; - femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); - femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum); - IssmDouble* beta0 = NULL; - IssmDouble* beta1 = NULL; - IssmDouble* phi = NULL; - femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); - femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum); - femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum); - femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); - - /*AR model spin-up with 0 noise to initialize ThermalforcingValuesAutoregressionEnum (688 = log(0.001)/log(0.99): decaying time of inluence of phi[0]=0.99 to 0.001 of beta_0*/ - int nspin = 688; - for(Object* &object:femmodel->elements->objects){ - Element* element = xDynamicCast(object); //generate element object - element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,FrontalForcingsRignotAutoregressionEnum); - } - - /*Cleanup*/ - xDelete(beta0); - xDelete(beta1); - xDelete(phi); -}/*}}}*/ void Thermalforcingautoregressionx(FemModel* femmodel){/*{{{*/ @@ -67,6 +38,4 @@ femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum); - /*Initialize module at first time step*/ - if(time<=starttime+dt){ThermalforcingautoregressionInitx(femmodel);} /*Determine if this is a time step for the AR model*/ bool isstepforar = false; @@ -81,16 +50,16 @@ bool isstochastic; bool istfstochastic = false; - int M,N,Nphi,arorder,numbasins,my_rank; + int M,N,Nlagcoefs,arorder,numbasins,my_rank; femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum); IssmDouble tinit_ar; - IssmDouble* beta0 = NULL; - IssmDouble* beta1 = NULL; - IssmDouble* phi = NULL; + IssmDouble* termconstant = NULL; + IssmDouble* trend = NULL; + IssmDouble* lagcoefs = NULL; femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum); - femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); + femmodel->parameters->FindParam(&termconstant,&M,FrontalForcingsautoregressionconstEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&trend,&M,FrontalForcingsautoregressiontrendEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&lagcoefs,&M,&Nlagcoefs,FrontalForcingsarlagcoefsEnum); _assert_(M==numbasins); _assert_(Nlagcoefs==arorder); femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); @@ -111,10 +80,10 @@ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum); + element->Autoregression(isstepforar,arorder,telapsed_ar,tstep_ar,termconstant,trend,lagcoefs,istfstochastic,FrontalForcingsRignotAutoregressionEnum); } /*Cleanup*/ - xDelete(beta0); - xDelete(beta1); - xDelete(phi); + xDelete(termconstant); + xDelete(trend); + xDelete(lagcoefs); }/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h (revision 27245) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h (revision 27246) @@ -7,5 +7,4 @@ /* local prototypes: */ void FrontalForcingsx(FemModel* femmodel); -void ThermalforcingautoregressionInitx(FemModel* femmodel); void Thermalforcingautoregressionx(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27245) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27246) @@ -260,12 +260,12 @@ parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.ar_initialtime",BasalforcingsAutoregressionInitialTimeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.ar_timestep",BasalforcingsAutoregressionTimestepEnum)); - iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.beta0"); - parameters->AddObject(new DoubleVecParam(BasalforcingsBeta0Enum,transparam,N)); - xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.beta1"); - parameters->AddObject(new DoubleVecParam(BasalforcingsBeta1Enum,transparam,N)); - xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.phi"); - parameters->AddObject(new DoubleMatParam(BasalforcingsPhiEnum,transparam,M,N)); + iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.const"); + parameters->AddObject(new DoubleVecParam(BasalforcingsautoregressionconstEnum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.trend"); + parameters->AddObject(new DoubleVecParam(BasalforcingsautoregressiontrendEnum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.arlag_coefs"); + parameters->AddObject(new DoubleMatParam(BasalforcingsarlagcoefsEnum,transparam,M,N)); xDelete(transparam); iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.upperwater_melting_rate"); @@ -432,12 +432,12 @@ parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_timestep",SmbAutoregressionTimestepEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_bins",SmbNumElevationBinsEnum)); - iomodel->FetchData(&transparam,&M,&N,"md.smb.beta0"); - parameters->AddObject(new DoubleVecParam(SmbBeta0Enum,transparam,N)); - xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.smb.beta1"); - parameters->AddObject(new DoubleVecParam(SmbBeta1Enum,transparam,N)); - xDelete(transparam); - iomodel->FetchData(&transparam,&M,&N,"md.smb.phi"); - parameters->AddObject(new DoubleMatParam(SmbPhiEnum,transparam,M,N)); + iomodel->FetchData(&transparam,&M,&N,"md.smb.const"); + parameters->AddObject(new DoubleVecParam(SmbautoregressionconstEnum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.smb.trend"); + parameters->AddObject(new DoubleVecParam(SmbautoregressiontrendEnum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.smb.arlag_coefs"); + parameters->AddObject(new DoubleMatParam(SmbarlagcoefsEnum,transparam,M,N)); xDelete(transparam); iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserates"); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 27245) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 27246) @@ -147,32 +147,4 @@ }/*}}}*/ -void SmbautoregressionInitx(FemModel* femmodel){/*{{{*/ - - /*Initialization step of Smbautoregressionx*/ - int M,N,Nphi,arorder,numbasins,my_rank; - IssmDouble starttime,tstep_ar,tinit_ar; - femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum); - femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum); - IssmDouble* beta0 = NULL; - IssmDouble* beta1 = NULL; - IssmDouble* phi = NULL; - femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); - femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum); - femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum); - femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); - - /*AR model spin-up with 0 noise to initialize SmbValuesAutoregressionEnum (688 = log(0.001)/log(0.99): decaying time of inluence of phi[0]=0.99 to 0.001 of beta_0*/ - int nspin = 688; - for(Object* &object:femmodel->elements->objects){ - Element* element = xDynamicCast(object); //generate element object - element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,SMBautoregressionEnum); - } - /*Cleanup*/ - xDelete(beta0); - xDelete(beta1); - xDelete(phi); -}/*}}}*/ void Smbautoregressionx(FemModel* femmodel){/*{{{*/ @@ -184,6 +156,4 @@ femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum); - /*Initialize module at first time step*/ - if(time<=starttime+dt){SmbautoregressionInitx(femmodel);} /*Determine if this is a time step for the AR model*/ bool isstepforar = false; @@ -198,12 +168,12 @@ bool isstochastic; bool issmbstochastic = false; - int M,N,Nphi,arorder,numbasins,numelevbins,my_rank; + int M,N,Nlagcoefs,arorder,numbasins,numelevbins,my_rank; femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum); femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum); femmodel->parameters->FindParam(&numelevbins,SmbNumElevationBinsEnum); IssmDouble tinit_ar; - IssmDouble* beta0 = NULL; - IssmDouble* beta1 = NULL; - IssmDouble* phi = NULL; + IssmDouble* termconstant = NULL; + IssmDouble* trend = NULL; + IssmDouble* lagcoefs = NULL; IssmDouble* lapserates = NULL; IssmDouble* elevbins = NULL; @@ -211,10 +181,10 @@ femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum); - femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins); - femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); - 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); - femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&termconstant,&M,SmbautoregressionconstEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&trend,&M,SmbautoregressiontrendEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&lagcoefs,&M,&Nlagcoefs,SmbarlagcoefsEnum); _assert_(M==numbasins); _assert_(Nlagcoefs==arorder); + 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); + femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum); _assert_(M==numbasins); femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); @@ -236,5 +206,5 @@ Element* element = xDynamicCast(object); /*Compute autoregression*/ - element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum); + element->Autoregression(isstepforar,arorder,telapsed_ar,tstep_ar,termconstant,trend,lagcoefs,issmbstochastic,SMBautoregressionEnum); /*Compute lapse rate adjustment*/ element->LapseRateBasinSMB(numelevbins,lapserates,elevbins,refelevation); @@ -242,7 +212,7 @@ /*Cleanup*/ - xDelete(beta0); - xDelete(beta1); - xDelete(phi); + xDelete(termconstant); + xDelete(trend); + xDelete(lagcoefs); xDelete(lapserates); xDelete(elevbins); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 27245) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 27246) @@ -13,5 +13,4 @@ void SmbGradientsx(FemModel* femmodel); void SmbGradientsElax(FemModel* femmodel); -void SmbautoregressionInitx(FemModel* femmodel); void Smbautoregressionx(FemModel* femmodel); void Delta18oParameterizationx(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27245) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27246) @@ -70,6 +70,6 @@ syn keyword cConstant BasalforcingsAutoregressionTimestepEnum syn keyword cConstant BasalforcingsAutoregressiveOrderEnum -syn keyword cConstant BasalforcingsBeta0Enum -syn keyword cConstant BasalforcingsBeta1Enum +syn keyword cConstant BasalforcingsautoregressionconstEnum +syn keyword cConstant BasalforcingsautoregressiontrendEnum syn keyword cConstant BasalforcingsBottomplumedepthEnum syn keyword cConstant BasalforcingsCrustthicknessEnum @@ -89,5 +89,5 @@ syn keyword cConstant BasalforcingsMantleconductivityEnum syn keyword cConstant BasalforcingsNusseltEnum -syn keyword cConstant BasalforcingsPhiEnum +syn keyword cConstant BasalforcingsarlagcoefsEnum syn keyword cConstant BasalforcingsPicoAverageOverturningEnum syn keyword cConstant BasalforcingsPicoAverageSalinityEnum @@ -199,9 +199,9 @@ syn keyword cConstant FrontalForcingsAutoregressionTimestepEnum syn keyword cConstant FrontalForcingsAutoregressiveOrderEnum -syn keyword cConstant FrontalForcingsBeta0Enum -syn keyword cConstant FrontalForcingsBeta1Enum +syn keyword cConstant FrontalForcingsautoregressionconstEnum +syn keyword cConstant FrontalForcingsautoregressiontrendEnum syn keyword cConstant FrontalForcingsNumberofBasinsEnum syn keyword cConstant FrontalForcingsParamEnum -syn keyword cConstant FrontalForcingsPhiEnum +syn keyword cConstant FrontalForcingsarlagcoefsEnum syn keyword cConstant GrdModelEnum syn keyword cConstant GroundinglineFrictionInterpolationEnum @@ -478,6 +478,6 @@ syn keyword cConstant SmbAutoregressiveOrderEnum syn keyword cConstant SmbAveragingEnum -syn keyword cConstant SmbBeta0Enum -syn keyword cConstant SmbBeta1Enum +syn keyword cConstant SmbautoregressionconstEnum +syn keyword cConstant SmbautoregressiontrendEnum syn keyword cConstant SmbDesfacEnum syn keyword cConstant SmbDpermilEnum @@ -516,5 +516,5 @@ syn keyword cConstant SmbNumRequestedOutputsEnum syn keyword cConstant SmbPfacEnum -syn keyword cConstant SmbPhiEnum +syn keyword cConstant SmbarlagcoefsEnum syn keyword cConstant SmbRdlEnum syn keyword cConstant SmbRefElevationEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27245) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27246) @@ -64,6 +64,6 @@ BasalforcingsAutoregressionTimestepEnum, BasalforcingsAutoregressiveOrderEnum, - BasalforcingsBeta0Enum, - BasalforcingsBeta1Enum, + BasalforcingsautoregressionconstEnum, + BasalforcingsautoregressiontrendEnum, BasalforcingsBottomplumedepthEnum, BasalforcingsCrustthicknessEnum, @@ -83,5 +83,5 @@ BasalforcingsMantleconductivityEnum, BasalforcingsNusseltEnum, - BasalforcingsPhiEnum, + BasalforcingsarlagcoefsEnum, BasalforcingsPicoAverageOverturningEnum, BasalforcingsPicoAverageSalinityEnum, @@ -193,9 +193,9 @@ FrontalForcingsAutoregressionTimestepEnum, FrontalForcingsAutoregressiveOrderEnum, - FrontalForcingsBeta0Enum, - FrontalForcingsBeta1Enum, + FrontalForcingsautoregressionconstEnum, + FrontalForcingsautoregressiontrendEnum, FrontalForcingsNumberofBasinsEnum, FrontalForcingsParamEnum, - FrontalForcingsPhiEnum, + FrontalForcingsarlagcoefsEnum, GrdModelEnum, GroundinglineFrictionInterpolationEnum, @@ -472,6 +472,6 @@ SmbAutoregressiveOrderEnum, SmbAveragingEnum, - SmbBeta0Enum, - SmbBeta1Enum, + SmbautoregressionconstEnum, + SmbautoregressiontrendEnum, SmbDesfacEnum, SmbDpermilEnum, @@ -510,5 +510,5 @@ SmbNumRequestedOutputsEnum, SmbPfacEnum, - SmbPhiEnum, + SmbarlagcoefsEnum, SmbRdlEnum, SmbRefElevationEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27245) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27246) @@ -72,6 +72,6 @@ case BasalforcingsAutoregressionTimestepEnum : return "BasalforcingsAutoregressionTimestep"; case BasalforcingsAutoregressiveOrderEnum : return "BasalforcingsAutoregressiveOrder"; - case BasalforcingsBeta0Enum : return "BasalforcingsBeta0"; - case BasalforcingsBeta1Enum : return "BasalforcingsBeta1"; + case BasalforcingsautoregressionconstEnum : return "Basalforcingsautoregressionconst"; + case BasalforcingsautoregressiontrendEnum : return "Basalforcingsautoregressiontrend"; case BasalforcingsBottomplumedepthEnum : return "BasalforcingsBottomplumedepth"; case BasalforcingsCrustthicknessEnum : return "BasalforcingsCrustthickness"; @@ -91,5 +91,5 @@ case BasalforcingsMantleconductivityEnum : return "BasalforcingsMantleconductivity"; case BasalforcingsNusseltEnum : return "BasalforcingsNusselt"; - case BasalforcingsPhiEnum : return "BasalforcingsPhi"; + case BasalforcingsarlagcoefsEnum : return "Basalforcingsarlagcoefs"; case BasalforcingsPicoAverageOverturningEnum : return "BasalforcingsPicoAverageOverturning"; case BasalforcingsPicoAverageSalinityEnum : return "BasalforcingsPicoAverageSalinity"; @@ -201,9 +201,9 @@ case FrontalForcingsAutoregressionTimestepEnum : return "FrontalForcingsAutoregressionTimestep"; case FrontalForcingsAutoregressiveOrderEnum : return "FrontalForcingsAutoregressiveOrder"; - case FrontalForcingsBeta0Enum : return "FrontalForcingsBeta0"; - case FrontalForcingsBeta1Enum : return "FrontalForcingsBeta1"; + case FrontalForcingsautoregressionconstEnum : return "FrontalForcingsautoregressionconst"; + case FrontalForcingsautoregressiontrendEnum : return "FrontalForcingsautoregressiontrend"; case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins"; case FrontalForcingsParamEnum : return "FrontalForcingsParam"; - case FrontalForcingsPhiEnum : return "FrontalForcingsPhi"; + case FrontalForcingsarlagcoefsEnum : return "FrontalForcingsarlagcoefs"; case GrdModelEnum : return "GrdModel"; case GroundinglineFrictionInterpolationEnum : return "GroundinglineFrictionInterpolation"; @@ -480,6 +480,6 @@ case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder"; case SmbAveragingEnum : return "SmbAveraging"; - case SmbBeta0Enum : return "SmbBeta0"; - case SmbBeta1Enum : return "SmbBeta1"; + case SmbautoregressionconstEnum : return "Smbautoregressionconst"; + case SmbautoregressiontrendEnum : return "Smbautoregressiontrend"; case SmbDesfacEnum : return "SmbDesfac"; case SmbDpermilEnum : return "SmbDpermil"; @@ -518,5 +518,5 @@ case SmbNumRequestedOutputsEnum : return "SmbNumRequestedOutputs"; case SmbPfacEnum : return "SmbPfac"; - case SmbPhiEnum : return "SmbPhi"; + case SmbarlagcoefsEnum : return "Smbarlagcoefs"; 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 27245) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27246) @@ -63,6 +63,6 @@ syn keyword juliaConstC BasalforcingsAutoregressionTimestepEnum syn keyword juliaConstC BasalforcingsAutoregressiveOrderEnum -syn keyword juliaConstC BasalforcingsBeta0Enum -syn keyword juliaConstC BasalforcingsBeta1Enum +syn keyword juliaConstC BasalforcingsautoregressionconstEnum +syn keyword juliaConstC BasalforcingsautoregressiontrendEnum syn keyword juliaConstC BasalforcingsBottomplumedepthEnum syn keyword juliaConstC BasalforcingsCrustthicknessEnum @@ -82,5 +82,5 @@ syn keyword juliaConstC BasalforcingsMantleconductivityEnum syn keyword juliaConstC BasalforcingsNusseltEnum -syn keyword juliaConstC BasalforcingsPhiEnum +syn keyword juliaConstC BasalforcingsarlagcoefsEnum syn keyword juliaConstC BasalforcingsPicoAverageOverturningEnum syn keyword juliaConstC BasalforcingsPicoAverageSalinityEnum @@ -192,9 +192,9 @@ syn keyword juliaConstC FrontalForcingsAutoregressionTimestepEnum syn keyword juliaConstC FrontalForcingsAutoregressiveOrderEnum -syn keyword juliaConstC FrontalForcingsBeta0Enum -syn keyword juliaConstC FrontalForcingsBeta1Enum +syn keyword juliaConstC FrontalForcingsautoregressionconstEnum +syn keyword juliaConstC FrontalForcingsautoregressiontrendEnum syn keyword juliaConstC FrontalForcingsNumberofBasinsEnum syn keyword juliaConstC FrontalForcingsParamEnum -syn keyword juliaConstC FrontalForcingsPhiEnum +syn keyword juliaConstC FrontalForcingsarlagcoefsEnum syn keyword juliaConstC GrdModelEnum syn keyword juliaConstC GroundinglineFrictionInterpolationEnum @@ -471,6 +471,6 @@ syn keyword juliaConstC SmbAutoregressiveOrderEnum syn keyword juliaConstC SmbAveragingEnum -syn keyword juliaConstC SmbBeta0Enum -syn keyword juliaConstC SmbBeta1Enum +syn keyword juliaConstC SmbautoregressionconstEnum +syn keyword juliaConstC SmbautoregressiontrendEnum syn keyword juliaConstC SmbDesfacEnum syn keyword juliaConstC SmbDpermilEnum @@ -509,5 +509,5 @@ syn keyword juliaConstC SmbNumRequestedOutputsEnum syn keyword juliaConstC SmbPfacEnum -syn keyword juliaConstC SmbPhiEnum +syn keyword juliaConstC SmbarlagcoefsEnum syn keyword juliaConstC SmbRdlEnum syn keyword juliaConstC SmbRefElevationEnum @@ -864,5 +864,5 @@ syn keyword juliaConstC SamplingBetaEnum syn keyword juliaConstC SamplingKappaEnum -syn keyword juliaConstC SamplingPhiEnum +syn keyword juliaConstC SamplingarlagcoefsEnum syn keyword juliaConstC SamplingTauEnum syn keyword juliaConstC SealevelEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27245) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27246) @@ -72,6 +72,6 @@ else if (strcmp(name,"BasalforcingsAutoregressionTimestep")==0) return BasalforcingsAutoregressionTimestepEnum; else if (strcmp(name,"BasalforcingsAutoregressiveOrder")==0) return BasalforcingsAutoregressiveOrderEnum; - else if (strcmp(name,"BasalforcingsBeta0")==0) return BasalforcingsBeta0Enum; - else if (strcmp(name,"BasalforcingsBeta1")==0) return BasalforcingsBeta1Enum; + else if (strcmp(name,"Basalforcingsautoregressionconst")==0) return BasalforcingsautoregressionconstEnum; + else if (strcmp(name,"Basalforcingsautoregressiontrend")==0) return BasalforcingsautoregressiontrendEnum; else if (strcmp(name,"BasalforcingsBottomplumedepth")==0) return BasalforcingsBottomplumedepthEnum; else if (strcmp(name,"BasalforcingsCrustthickness")==0) return BasalforcingsCrustthicknessEnum; @@ -91,5 +91,5 @@ else if (strcmp(name,"BasalforcingsMantleconductivity")==0) return BasalforcingsMantleconductivityEnum; else if (strcmp(name,"BasalforcingsNusselt")==0) return BasalforcingsNusseltEnum; - else if (strcmp(name,"BasalforcingsPhi")==0) return BasalforcingsPhiEnum; + else if (strcmp(name,"Basalforcingsarlagcoefs")==0) return BasalforcingsarlagcoefsEnum; else if (strcmp(name,"BasalforcingsPicoAverageOverturning")==0) return BasalforcingsPicoAverageOverturningEnum; else if (strcmp(name,"BasalforcingsPicoAverageSalinity")==0) return BasalforcingsPicoAverageSalinityEnum; @@ -204,9 +204,9 @@ else if (strcmp(name,"FrontalForcingsAutoregressionTimestep")==0) return FrontalForcingsAutoregressionTimestepEnum; else if (strcmp(name,"FrontalForcingsAutoregressiveOrder")==0) return FrontalForcingsAutoregressiveOrderEnum; - else if (strcmp(name,"FrontalForcingsBeta0")==0) return FrontalForcingsBeta0Enum; - else if (strcmp(name,"FrontalForcingsBeta1")==0) return FrontalForcingsBeta1Enum; + else if (strcmp(name,"FrontalForcingsautoregressionconst")==0) return FrontalForcingsautoregressionconstEnum; + else if (strcmp(name,"FrontalForcingsautoregressiontrend")==0) return FrontalForcingsautoregressiontrendEnum; else if (strcmp(name,"FrontalForcingsNumberofBasins")==0) return FrontalForcingsNumberofBasinsEnum; else if (strcmp(name,"FrontalForcingsParam")==0) return FrontalForcingsParamEnum; - else if (strcmp(name,"FrontalForcingsPhi")==0) return FrontalForcingsPhiEnum; + else if (strcmp(name,"FrontalForcingsarlagcoefs")==0) return FrontalForcingsarlagcoefsEnum; else if (strcmp(name,"GrdModel")==0) return GrdModelEnum; else if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum; @@ -489,6 +489,6 @@ else if (strcmp(name,"SmbAutoregressiveOrder")==0) return SmbAutoregressiveOrderEnum; else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum; - else if (strcmp(name,"SmbBeta0")==0) return SmbBeta0Enum; - else if (strcmp(name,"SmbBeta1")==0) return SmbBeta1Enum; + else if (strcmp(name,"Smbautoregressionconst")==0) return SmbautoregressionconstEnum; + else if (strcmp(name,"Smbautoregressiontrend")==0) return SmbautoregressiontrendEnum; else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; @@ -530,5 +530,5 @@ else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum; else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum; - else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum; + else if (strcmp(name,"Smbarlagcoefs")==0) return SmbarlagcoefsEnum; else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum; else if (strcmp(name,"SmbRefElevation")==0) return SmbRefElevationEnum; Index: /issm/trunk-jpl/src/m/classes/SMBautoregression.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBautoregression.m (revision 27245) +++ /issm/trunk-jpl/src/m/classes/SMBautoregression.m (revision 27246) @@ -7,10 +7,10 @@ properties (SetAccess=public) num_basins = 0; - beta0 = NaN; - beta1 = NaN; + const = NaN; + trend = NaN; ar_order = 0; ar_initialtime = 0; ar_timestep = 0; - phi = NaN; + arlag_coefs = NaN; basin_id = NaN; lapserates = NaN; @@ -37,11 +37,11 @@ end % }}} function self = initialize(self,md) % {{{ - if isnan(self.beta1) - self.beta1 = zeros(1,self.num_basins); %no trend in SMB - disp(' smb.beta1 (trend) not specified: value set to 0'); + if isnan(self.trend) + self.trend = zeros(1,self.num_basins); %no trend in SMB + disp(' smb.trend (trend) not specified: value set to 0'); end if (self.ar_order==0) self.ar_order = 1; %dummy 1 value for autoregression - self.phi = zeros(self.num_basins,self.ar_order); %autoregression coefficients all set to 0 + self.arlag_coefs = zeros(self.num_basins,self.ar_order); %autoregression coefficients all set to 0 disp(' smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0'); end @@ -54,7 +54,7 @@ disp(' smb.ar_timestep (timestep of autoregressive model) not specified: set to md.timestepping.time_step'); end - if isnan(self.phi) - self.phi = zeros(self.num_basins,self.ar_order); %autoregression model of order 0 - disp(' smb.phi (lag coefficients) not specified: order of autoregressive model set to 0'); + if isnan(self.arlag_coefs) + self.arlag_coefs = zeros(self.num_basins,self.ar_order); %autoregression model of order 0 + disp(' smb.arlag_coefs (lag coefficients) not specified: order of autoregressive model set to 0'); end end % }}} @@ -67,10 +67,10 @@ md = checkfield(md,'fieldname','smb.num_basins','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.beta0','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.beta1','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.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 md = checkfield(md,'fieldname','smb.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','smb.ar_initialtime','numel',1,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','smb.ar_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.phi','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]); + md = checkfield(md,'fieldname','smb.arlag_coefs','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]); if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1) @@ -100,10 +100,10 @@ fielddisplay(self,'num_basins','number of different basins [unitless]'); fielddisplay(self,'basin_id','basin number assigned to each element [unitless]'); - fielddisplay(self,'beta0','basin-specific intercept values [m ice eq./yr] (if beta_1==0 mean=beta_0/(1-sum(phi)))'); - fielddisplay(self,'beta1','basin-specific trend values [m ice eq. yr^(-2)]'); + 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,'ar_order','order of the autoregressive model [unitless]'); fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]'); fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); - fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]'); + fielddisplay(self,'arlag_coefs','basin-specific vectors of lag coefficients [unitless]'); fielddisplay(self,'lapserates','basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin [m ice eq yr^-1 m^-1] (default: no lapse rate)'); fielddisplay(self,'elevationbins','basin-specific separations between elevation bins, 1 row per basin, 1 column per limit between bins [m] (default: no basin separation)'); @@ -153,7 +153,7 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_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','beta0','format','DoubleMat','name','md.smb.beta0','scale',1./yts,'yts',yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',yts); + 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,'object',self,'class','smb','fieldname','arlag_coefs','format','DoubleMat','name','md.smb.arlag_coefs','yts',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/SMBautoregression.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBautoregression.py (revision 27245) +++ /issm/trunk-jpl/src/m/classes/SMBautoregression.py (revision 27246) @@ -16,10 +16,10 @@ def __init__(self, *args): # {{{ self.num_basins = 0 - self.beta0 = np.nan - self.beta1 = np.nan + self.const = np.nan + self.trend = np.nan self.ar_order = 0 self.ar_initialtime = 0 self.ar_timestep = 0 - self.phi = np.nan + self.arlag_coefs = np.nan self.basin_id = np.nan self.lapserates = np.nan @@ -41,10 +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, 'beta0', 'basin-specific intercept values [m ice eq./yr] (if beta_1==0 mean=beta_0/(1-sum(phi)))')) - s += '{}\n'.format(fielddisplay(self, 'beta1', 'basin-specific trend values [m ice eq. yr^(-2)]')) + 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, 'ar_order', 'order of the autoregressive model [unitless]')) s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]')) s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) - s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of lag coefficients [unitless]')) s += '{}\n'.format(fielddisplay(self, 'lapserates', 'basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin [m ice eq yr^-1 m^-1] (default: no lapse rate)')) s += '{}\n'.format(fielddisplay(self, 'elevationbins', 'basin-specific SMB lapse rates applied in range of SMB<0 [m ice eq yr^-1 m^-1] (default: no lapse rate)')) @@ -72,10 +72,10 @@ def initialize(self, md): # {{{ - if np.all(np.isnan(self.beta1)): - self.beta1 = np.zeros((1, self.num_basins)) # No trend in SMB - print(' smb.beta1 (trend) not specified: value set to 0') + if np.all(np.isnan(self.trend)): + self.trend = np.zeros((1, self.num_basins)) # No trend in SMB + print(' smb.trend (trend) not specified: value set to 0') if self.ar_order == 0: self.ar_order = 1 # Dummy 1 value for autoregression - self.phi = np.zeros((self.num_basins, self.ar_order)) # Autorgression coefficients all set to 0 + self.arlag_coefs = np.zeros((self.num_basins, self.ar_order)) # Autorgression coefficients all set to 0 print(' smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0') if self.ar_initialtime == 0: @@ -85,7 +85,7 @@ self.ar_timestep = md.timestepping.time_step # Autoregression model has no prescribed time step print(' smb.ar_timestep (timestep of autoregressive model) not specified: set to md.timestepping.time_step') - if np.all(np.isnan(self.phi)): - self.phi = np.zeros((self.num_basins, self.ar_order)) # Autoregression model of order 0 - print(' smb.phi (lag coefficients) not specified: order of autoregressive model set to 0') + if np.all(np.isnan(self.arlag_coefs)): + self.arlag_coefs = np.zeros((self.num_basins, self.ar_order)) # Autoregression model of order 0 + print(' smb.arlag_coefs (lag coefficients) not specified: order of autoregressive model set to 0') return self # }}} @@ -95,13 +95,13 @@ md = checkfield(md, 'fieldname', 'smb.num_basins', '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.beta0)) == 1: - self.beta0 = np.array([self.beta0]) - self.beta1 = np.array([self.beta1]) - md = checkfield(md, 'fieldname', 'smb.beta0', '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.beta1', '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.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 md = checkfield(md, 'fieldname', 'smb.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) md = checkfield(md, 'fieldname', 'smb.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'smb.ar_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.phi', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order]) + md = checkfield(md, 'fieldname', 'smb.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order]) if(np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1): @@ -168,7 +168,7 @@ WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ar_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', 'beta0', 'format', 'DoubleMat', 'name', 'md.smb.beta0', 'scale', 1 / yts, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.smb.beta1', 'scale', 1 / (yts ** 2), 'yts', yts) - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.smb.phi', 'yts', yts) + 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, 'object', self, 'class', 'smb', 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.smb.arlag_coefs', 'yts', 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/autoregressionlinearbasalforcings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.m (revision 27245) +++ /issm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.m (revision 27246) @@ -10,10 +10,10 @@ properties (SetAccess=public) num_basins = 0; - beta0 = NaN; - beta1 = NaN; + const = NaN; + trend = NaN; ar_order = 0; ar_initialtime = 0; ar_timestep = 0; - phi = NaN; + arlag_coefs = NaN; basin_id = NaN; groundedice_melting_rate = NaN; @@ -59,10 +59,10 @@ 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.beta0','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); - md = checkfield(md,'fieldname','basalforcings.beta1','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); + 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); md = checkfield(md,'fieldname','basalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','basalforcings.ar_initialtime','numel',1,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','basalforcings.ar_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.phi','NaN',1,'Inf',1,'size',[md.basalforcings.num_basins,md.basalforcings.ar_order]); + md = checkfield(md,'fieldname','basalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.basalforcings.num_basins,md.basalforcings.ar_order]); end if ismember('BalancethicknessAnalysis',analyses), @@ -82,10 +82,10 @@ fielddisplay(self,'num_basins','number of different basins [unitless]'); fielddisplay(self,'basin_id','basin number assigned to each element [unitless]'); - fielddisplay(self,'beta0','basin-specific intercept values [m/yr] (if beta_1==0 mean=beta_0/(1-sum(phi)))'); - fielddisplay(self,'beta1','basin-specific trend values [m yr^(-2)]'); + 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,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]'); fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); - fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]'); + fielddisplay(self,'arlag_coefs','basin-specific vectors of lag coefficients [unitless]'); fielddisplay(self,'deepwater_elevation','basin-specific elevation of ocean deepwater [m]'); fielddisplay(self,'upperwater_melting_rate','basin-specific basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]'); @@ -107,7 +107,7 @@ WriteData(fid,prefix,'object',self,'fieldname','ar_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','beta0','format','DoubleMat','name','md.basalforcings.beta0','scale',1./yts,'yts',yts); - WriteData(fid,prefix,'object',self,'fieldname','beta1','format','DoubleMat','name','md.basalforcings.beta1','scale',1./(yts^2),'yts',yts); - WriteData(fid,prefix,'object',self,'fieldname','phi','format','DoubleMat','name','md.basalforcings.phi','yts',yts); + 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','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/autoregressionlinearbasalforcings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.py (revision 27245) +++ /issm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.py (revision 27246) @@ -15,10 +15,10 @@ def __init__(self, *args): # {{{ self.num_basins = 0 - self.beta0 = np.nan - self.beta1 = np.nan + self.const = np.nan + self.trend = np.nan self.ar_order = 0 self.ar_initialtime = 0 self.ar_timestep = 0 - self.phi = np.nan + self.arlag_coefs = np.nan self.basin_id = np.nan self.groundedice_melting_rate = np.nan @@ -40,10 +40,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, 'beta0', 'basin-specific intercept values [m ice eq./yr] (if beta_1==0 mean=beta_0/(1-sum(phi)))')) - s += '{}\n'.format(fielddisplay(self, 'beta1', 'basin-specific trend values [m ice eq. yr^(-2)]')) + 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, 'ar_order', 'order of the autoregressive model [unitless]')) s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]')) s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) - s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of lag coefficients [unitless]')) s += '{}\n'.format(fielddisplay(self, 'deepwater_elevation', 'basin-specific elevation of ocean deepwater [m]')) s += '{}\n'.format(fielddisplay(self, 'upperwater_melting_rate', 'basin-specic basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]')) @@ -84,14 +84,14 @@ md = checkfield(md, 'fieldname', 'basalforcings.basin_id', 'Inf', 1, '>=', 0, '<=', md.basalforcings.num_basins, 'size', [md.mesh.numberofelements]) - if len(np.shape(self.beta0)) == 1: - self.beta0 = np.array([self.beta0]) - self.beta1 = np.array([self.beta1]) + if len(np.shape(self.const)) == 1: + self.const = np.array([self.const]) + self.trend = np.array([self.trend]) - md = checkfield(md, 'fieldname', 'basalforcings.beta0', '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.beta1', '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.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.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'basalforcings.ar_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.phi', 'NaN', 1, 'Inf', 1, 'size', [md.basalforcings.num_basins, md.basalforcings.ar_order]) + md = checkfield(md, 'fieldname', 'basalforcings.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.basalforcings.num_basins, md.basalforcings.ar_order]) if 'BalancethicknessAnalysis' in analyses: raise Exception('not implemented yet!') @@ -113,7 +113,7 @@ WriteData(fid, prefix, 'object', self, 'fieldname', 'ar_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', 'beta0', 'format', 'DoubleMat', 'name', 'md.basalforcings.beta0', 'scale', 1 / yts, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.basalforcings.beta1', 'scale', 1 / (yts ** 2), 'yts', yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.basalforcings.phi', 'yts', yts) + 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', '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/src/m/classes/frontalforcingsrignotautoregression.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.m (revision 27245) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.m (revision 27246) @@ -7,10 +7,10 @@ properties (SetAccess=public) num_basins = 0; - beta0 = NaN; - beta1 = NaN; + const = NaN; + trend = NaN; ar_order = 0; ar_initialtime = 0; ar_timestep = 0; - phi = NaN; + arlag_coef s = NaN; basin_id = NaN; subglacial_discharge = NaN; @@ -53,10 +53,10 @@ 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.beta0','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins); - md = checkfield(md,'fieldname','frontalforcings.beta1','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins); + 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); md = checkfield(md,'fieldname','frontalforcings.ar_initialtime','numel',1,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','frontalforcings.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep - md = checkfield(md,'fieldname','frontalforcings.phi','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]); + md = checkfield(md,'fieldname','frontalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]); end % }}} @@ -66,10 +66,10 @@ 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,'beta0','basin-specific intercept values [∘C] (if beta_1==0 mean=beta_0/(1-sum(phi)))'); - fielddisplay(self,'beta1','basin-specific trend values [∘C yr^(-1)]'); + fielddisplay(self,'const','basin-specific constant term [∘C]'); + fielddisplay(self,'trend','basin-specific trend values [∘C yr^(-1)]'); fielddisplay(self,'ar_order','order of the autoregressive model [unitless]'); fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]'); fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); - fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]'); + fielddisplay(self,'arlag_coefs','basin-specific vectors of lag coefficients [unitless]'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -82,7 +82,7 @@ WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_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','beta0','format','DoubleMat','name','md.frontalforcings.beta0'); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','beta1','format','DoubleMat','name','md.frontalforcings.beta1','scale',1./yts,'yts',yts); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','phi','format','DoubleMat','name','md.frontalforcings.phi','yts',yts); + 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,'object',self,'class','frontalforcings','fieldname','arlag_coefs','format','DoubleMat','name','md.frontalforcings.arlag_coefs','yts',yts); end % }}} end Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.py (revision 27245) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.py (revision 27246) @@ -16,10 +16,10 @@ def __init__(self, *args): # {{{ self.num_basins = 0 - self.beta0 = np.nan - self.beta1 = np.nan + self.const = np.nan + self.trend = np.nan self.ar_order = 0 self.ar_initialtime = 0 self.ar_timestep = 0 - self.phi = np.nan + self.arlag_coefs = np.nan self.basin_id = np.nan self.subglacial_discharge = np.nan @@ -35,10 +35,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, 'beta0', 'basin-specific intercept values [°C] (if beta_1==0 mean=beta_0/(1-sum(phi)))')) - s += '{}\n'.format(fielddisplay(self, 'beta1', 'basin-specific trend values [°C yr^(-1)]')) + 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, 'ar_order', 'order of the autoregressive model [unitless]')) s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]')) s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) - s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of lag coefficients [unitless]')) return s #}}} @@ -60,13 +60,13 @@ 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.beta0)) == 1: - self.beta0 = np.array([self.beta0]) - self.beta1 = np.array([self.beta1]) - md = checkfield(md, 'fieldname', 'frontalforcings.beta0', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins) - md = checkfield(md, 'fieldname', 'frontalforcings.beta1', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins) + 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) md = checkfield(md, 'fieldname', 'frontalforcings.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) md = checkfield(md, 'fieldname', 'frontalforcings.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'frontalforcings.ar_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep - md = checkfield(md, 'fieldname', 'frontalforcings.phi', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ar_order]) + md = checkfield(md, 'fieldname', 'frontalforcings.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ar_order]) return md # }}} @@ -86,6 +86,6 @@ WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ar_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', 'beta0', 'format', 'DoubleMat', 'name', 'md.frontalforcings.beta0') - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.frontalforcings.beta1', 'scale', 1 / yts, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.frontalforcings.phi', 'yts', yts) + 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, 'object', self, 'class', 'frontalforcings', 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.frontalforcings.arlag_coefs', 'yts', yts) # }}} Index: /issm/trunk-jpl/test/NightlyRun/test257.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.m (revision 27245) +++ /issm/trunk-jpl/test/NightlyRun/test257.m (revision 27246) @@ -39,10 +39,10 @@ md.smb.num_basins = 3; %number of basins md.smb.basin_id = idbasin; %prescribe basin ID number to elements -md.smb.beta0 = [0.5,1.2,1.5]; %intercept values of SMB in basins [m ice eq./yr] -md.smb.beta1 = [0.0,0.01,-0.01]; %trend values of SMB in basins [m ice eq./yr^2] +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.ar_initialtime = md.timestepping.start_time; md.smb.ar_order = 4; md.smb.ar_timestep = 2.0; %timestep of the autoregressive model [yr] -md.smb.phi = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1];[0.4,-0.4,0.1,-0.1]]; +md.smb.arlag_coefs = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1];[0.4,-0.4,0.1,-0.1]]; md.smb.lapserates = [0.01,0.0;0.01,-0.01;0.0,-0.01]; md.smb.elevationbins = [100;150;100]; Index: /issm/trunk-jpl/test/NightlyRun/test257.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.py (revision 27245) +++ /issm/trunk-jpl/test/NightlyRun/test257.py (revision 27246) @@ -44,10 +44,10 @@ md.smb.num_basins = 3 # number of basins md.smb.basin_id = idbasin # prescribe basin ID number to elements; -md.smb.beta0 = np.array([[0.5, 1.2, 1.5]]) # intercept values of SMB in basins [m ice eq./yr] -md.smb.beta1 = np.array([[0.0, 0.01, -0.01]]) # trend values of SMB in basins [m ice eq./yr^2] +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.ar_initialtime = md.timestepping.start_time md.smb.ar_order = 4 md.smb.ar_timestep = 2.0 #timestep of the autoregressive model [yr] -md.smb.phi = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1], [0.4, -0.4, 0.1, -0.1]]) +md.smb.arlag_coefs = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1], [0.4, -0.4, 0.1, -0.1]]) md.smb.lapserates = np.array([[0.01,0.0],[0.01,-0.01],[0.0,-0.01]]) md.smb.elevationbins = np.array([100,150,100]).reshape(md.smb.num_basins,1) Index: /issm/trunk-jpl/test/NightlyRun/test543.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.m (revision 27245) +++ /issm/trunk-jpl/test/NightlyRun/test543.m (revision 27246) @@ -49,10 +49,10 @@ md.frontalforcings.basin_id = idb_tf; md.frontalforcings.subglacial_discharge = 0.1*ones(md.mesh.numberofvertices,1); -md.frontalforcings.beta0 = [0.05,0.01]; %intercept values of TF in basins [C] -md.frontalforcings.beta1 = [0.001,0.0001]; %trend values of TF in basins [C/yr] +md.frontalforcings.const = [0.05,0.01]; %intercept values of TF in basins [C] +md.frontalforcings.trend = [0.001,0.0001]; %trend values of TF in basins [C/yr] md.frontalforcings.ar_initialtime = md.timestepping.start_time; %initial time in the AR model parameterization [yr] md.frontalforcings.ar_order = 4; md.frontalforcings.ar_timestep = 2; %timestep of the autoregressive model [yr] -md.frontalforcings.phi = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters +md.frontalforcings.arlag_coefs = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters % Floating Ice Melt parameters Index: /issm/trunk-jpl/test/NightlyRun/test543.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.py (revision 27245) +++ /issm/trunk-jpl/test/NightlyRun/test543.py (revision 27246) @@ -54,10 +54,10 @@ md.frontalforcings.basin_id = idb_tf md.frontalforcings.subglacial_discharge = 0.1 * np.ones((md.mesh.numberofvertices,)) -md.frontalforcings.beta0 = np.array([[0.05, 0.01]]) # intercept values of TF in basins [C] -md.frontalforcings.beta1 = np.array([[0.001, 0.0001]]) # trend values of TF in basins [C/yr] +md.frontalforcings.const = np.array([[0.05, 0.01]]) # intercept values of TF in basins [C] +md.frontalforcings.trend = np.array([[0.001, 0.0001]]) # trend values of TF in basins [C/yr] md.frontalforcings.ar_initialtime = md.timestepping.start_time # initial time in the AR model parameterization [yr] md.frontalforcings.ar_order = 4 md.frontalforcings.ar_timestep = 2 # timestep of the autoregressive model [yr] -md.frontalforcings.phi = np.array([[0.1, -0.1, 0.01, -0.01], [0.2, -0.2, 0.1, 0.0]]) # autoregressive parameters +md.frontalforcings.arlag_coefs = np.array([[0.1, -0.1, 0.01, -0.01], [0.2, -0.2, 0.1, 0.0]]) # autoregressive parameters #Floating Ice Melt parameters md.basalforcings.floatingice_melting_rate = 0.1 * np.ones((md.mesh.numberofvertices,)) Index: /issm/trunk-jpl/test/NightlyRun/test544.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test544.m (revision 27245) +++ /issm/trunk-jpl/test/NightlyRun/test544.m (revision 27246) @@ -27,10 +27,10 @@ md.smb.num_basins = nb_bas; %number of basins md.smb.basin_id = idb; %prescribe basin ID number to elements -md.smb.beta0 = [0.5,1.2]; %intercept values of SMB in basins [m ice eq./yr] -md.smb.beta1 = [0.0,0.01]; %trend values of SMB in basins [m ice eq./yr^2] +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.ar_initialtime = md.timestepping.start_time; md.smb.ar_order = 4; md.smb.ar_timestep = 2.0; %timestep of the autoregressive model [yr] -md.smb.phi = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1]]; +md.smb.arlag_coefs = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1]]; %Calving @@ -45,10 +45,10 @@ md.basalforcings.num_basins = nb_bas; %number of basins md.basalforcings.basin_id = idb; %prescribe basin ID number to elements -md.basalforcings.beta0 = [1.0,2.50]; %intercept values of DeepwaterMelt in basins [m/yr] -md.basalforcings.beta1 = [0.2,0.01]; %trend values of DeepwaterMelt in basins [m/yr^2] +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.ar_initialtime = md.timestepping.start_time; md.basalforcings.ar_order = 1; md.basalforcings.ar_timestep = 1.0; %timestep of the autoregressive model [yr] -md.basalforcings.phi = [0.0;0.1]; +md.basalforcings.arlag_coefs = [0.0;0.1]; md.basalforcings.deepwater_elevation = [-1000,-1520]; md.basalforcings.upperwater_elevation = [0,-50]; Index: /issm/trunk-jpl/test/NightlyRun/test544.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test544.py (revision 27245) +++ /issm/trunk-jpl/test/NightlyRun/test544.py (revision 27246) @@ -37,10 +37,10 @@ md.smb.num_basins = nb_bas # number of basins md.smb.basin_id = idb # prescribe basin ID number to elements; -md.smb.beta0 = np.array([[0.5, 1.2]]) # intercept values of SMB in basins [m ice eq./yr] -md.smb.beta1 = np.array([[0.0, 0.01]]) # trend values of SMB in basins [m ice eq./yr^2] +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.ar_initialtime = md.timestepping.start_time md.smb.ar_order = 4 md.smb.ar_timestep = 2.0 #timestep of the autoregressive model [yr] -md.smb.phi = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1]]) +md.smb.arlag_coefs = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1]]) #Calving @@ -55,10 +55,10 @@ md.basalforcings.num_basins = nb_bas md.basalforcings.basin_id = idb -md.basalforcings.beta0 = np.array([[1.0, 2.50]]) # intercept values of DeepwaterMelt in basins [m/yr] -md.basalforcings.beta1 = np.array([[0.2, 0.01]]) # trend values of DeepwaterMelt in basins [m/yr^2] +md.basalforcings.const = np.array([[1.0, 2.50]]) # intercept values of DeepwaterMelt in basins [m/yr] +md.basalforcings.trend = np.array([[0.2, 0.01]]) # trend values of DeepwaterMelt in basins [m/yr^2] md.basalforcings.ar_initialtime = md.timestepping.start_time # initial time in the AR model parameterization [yr] md.basalforcings.ar_order = 1 md.basalforcings.ar_timestep = 1.0 # timestep of the autoregressive model [yr] -md.basalforcings.phi = np.array([[0.0], [0.1]]) # autoregressive parameters +md.basalforcings.arlag_coefs = np.array([[0.0], [0.1]]) # autoregressive parameters md.basalforcings.deepwater_elevation = np.array([[-1000, -1520]]) md.basalforcings.upperwater_elevation = np.array([[0, -50]])