Changeset 27318

Timestamp:

10/21/22 11:51:59 (3 years ago)

Author:

vverjans

Message:

NEW piecewise polynomials implemented in ARMA model schemes

Location:

issm/trunk-jpl

Files:

• src/c/analyses/LevelsetAnalysis.cpp (modified) (1 diff)
• src/c/classes/Elements/Element.cpp (modified) (5 diffs)
• src/c/classes/Elements/Element.h (modified) (1 diff)
• src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (modified) (4 diffs)
• src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (modified) (4 diffs)
• src/c/modules/ModelProcessorx/CreateParameters.cpp (modified) (2 diffs)
• src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (4 diffs)
• src/c/shared/Enum/Enum.vim (modified) (18 diffs)
• src/c/shared/Enum/EnumDefinitions.h (modified) (5 diffs)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (5 diffs)
• src/c/shared/Enum/Enumjl.vim (modified) (8 diffs)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (17 diffs)
• src/m/classes/SMBarma.m (modified) (7 diffs)
• src/m/classes/SMBarma.py (modified) (6 diffs)
• src/m/classes/frontalforcingsrignotarma.m (modified) (5 diffs)
• src/m/classes/frontalforcingsrignotarma.py (modified) (4 diffs)
• src/m/classes/linearbasalforcingsarma.m (modified) (7 diffs)
• src/m/classes/linearbasalforcingsarma.py (modified) (8 diffs)
• test/Archives/Archive257.arch (modified) ( previous)
• test/Archives/Archive543.arch (modified) ( previous)
• test/Archives/Archive544.arch (modified) ( previous)
• test/NightlyRun/test257.m (modified) (3 diffs)
• test/NightlyRun/test257.py (modified) (3 diffs)
• test/NightlyRun/test543.m (modified) (1 diff)
• test/NightlyRun/test543.py (modified) (1 diff)
• test/NightlyRun/test544.m (modified) (2 diffs)
• test/NightlyRun/test544.py (modified) (3 diffs)

issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

@@ -246 +246 @@
                 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<IssmDouble>(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<IssmDouble>(transparam);
          iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.arlag_coefs");

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -61 +61 @@
         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<IssmDouble>(arorder);
    IssmDouble* malagcoefs_basin     = xNew<IssmDouble>(maorder);
+   IssmDouble* datebreaks_basin     = xNew<IssmDouble>(numbreaks);
+   IssmDouble* polyparams_basin     = xNew<IssmDouble>(numperiods*numparams);
    IssmDouble* varlist              = xNew<IssmDouble>(numvertices);
+   IssmDouble* sumpoly              = xNewZeroInit<IssmDouble>(arorder+1);
    IssmDouble* valuesautoregression = NULL;
    IssmDouble* valuesmovingaverage  = NULL;
@@ -104 +108 @@
    /*Get basin coefficients*/
    this->GetInputValue(&basinid,basinenum_type);
-        for(int ii=0;ii<arorder;ii++) arlagcoefs_basin[ii] = arlagcoefs[basinid*arorder+ii];
-        for(int ii=0;ii<maorder;ii++) malagcoefs_basin[ii] = malagcoefs[basinid*maorder+ii];
-   termconstant_basin   = termconstant[basinid];
-   trend_basin          = trend[basinid];
-
+        for(int i=0;i<arorder;i++) arlagcoefs_basin[i]   = arlagcoefs[basinid*arorder+i];
+        for(int i=0;i<maorder;i++) malagcoefs_basin[i]   = malagcoefs[basinid*maorder+i];
+        for(int i=0;i<numparams;i++){
+                for(int j=0;j<numperiods;j++) polyparams_basin[i*numperiods+j] = polyparams[basinid*numparams*numperiods+i*numperiods+j];
+        }
+        if(numbreaks>0){
+                for(int i=0;i<numbreaks;i++) datebreaks_basin[i] = datebreaks[basinid*numbreaks+i];
+        }
+
+        /*Compute terms from polynomial parameters from arorder steps back to present*/
+        IssmDouble telapsed_break;
+        IssmDouble tatstep;
+        for(int s=0;s<arorder+1;s++){
+                tatstep = time-s*tstep_arma;
+                if(numbreaks>0){
+                        /*Find index of tatstep compared to the breakpoints*/
+                        indperiod = 0;
+                        for(int i=0;i<numbreaks;i++){
+                                if(tatstep>datebreaks_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<numparams;j++)   sumpoly[s] = sumpoly[s]+polyparams_basin[indperiod+j*numperiods]*pow(telapsed_break,j);
+                }
+                else for(int j=0;j<numparams;j++) sumpoly[s] = sumpoly[s]+polyparams_basin[j*numperiods]*pow(tatstep,j);
+        }
         /*Initialze autoregressive and moving-average values at first time step*/
         if(time<=starttime+dt){
                 IssmDouble* initvaluesautoregression = xNewZeroInit<IssmDouble>(numvertices*arorder);
                 IssmDouble* initvaluesmovingaverage  = xNewZeroInit<IssmDouble>(numvertices*maorder);
-                for(int i=0;i<numvertices*arorder;i++) initvaluesautoregression[i]=termconstant_basin;
+                for(int i=0;i<numvertices*arorder;i++) initvaluesautoregression[i]=polyparams_basin[0];
       this->inputs->SetArrayInput(arenum_type,this->lid,initvaluesautoregression,numvertices*arorder);
       this->inputs->SetArrayInput(maenum_type,this->lid,initvaluesmovingaverage,numvertices*maorder);
@@ -141 +167 @@
          /*Compute autoregressive term*/
          IssmDouble autoregressionterm=0.;
-         for(int ii=0;ii<arorder;ii++) autoregressionterm += arlagcoefs_basin[ii]*(valuesautoregression[v+ii*numvertices]-(termconstant_basin+trend_basin*(telapsed_arma-(ii+1)*tstep_arma)));
-         /*Compute moving-average term*/
+         for(int ii=0;ii<arorder;ii++) autoregressionterm += arlagcoefs_basin[ii]*(valuesautoregression[v+ii*numvertices]-sumpoly[ii+1]);
+                        /*Compute moving-average term*/
          IssmDouble movingaverageterm=0.;
          for(int ii=0;ii<maorder;ii++) movingaverageterm  += malagcoefs_basin[ii]*valuesmovingaverage[v+ii*numvertices];
 
                         /*Stochastic variable value*/
-         varlist[v] = termconstant_basin+trend_basin*telapsed_arma+autoregressionterm+movingaverageterm+noiseterm;
+         varlist[v] = sumpoly[0]+autoregressionterm+movingaverageterm+noiseterm;
       }
       /*Update autoregression and moving-average values*/
@@ -169 +195 @@
    xDelete<IssmDouble>(arlagcoefs_basin);
    xDelete<IssmDouble>(malagcoefs_basin);
+   xDelete<IssmDouble>(datebreaks_basin);
+   xDelete<IssmDouble>(polyparams_basin);
+   xDelete<IssmDouble>(sumpoly);
    xDelete<IssmDouble>(varlist);
    xDelete<IssmDouble>(valuesautoregression);
    xDelete<IssmDouble>(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<IssmDouble>(arorder);
-   IssmDouble* varlist         = xNew<IssmDouble>(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<arorder;ii++) lagcoefs_basin[ii] = lagcoefs[basinid*arorder+ii];
-   termconstant_basin   = termconstant[basinid];
-   trend_basin   = trend[basinid];
-
-        /*Initialze autoregressive values at first time step*/
-        if(time<=starttime+dt){
-                IssmDouble* initvaluesautoregression = xNewZeroInit<IssmDouble>(numvertices*arorder);
-                for(int i=0;i<numvertices*arorder;i++) initvaluesautoregression[i]=termconstant_basin;
-      this->inputs->SetArrayInput(arenum_type,this->lid,initvaluesautoregression,numvertices*arorder);
-      xDelete<IssmDouble>(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;i++) varlist[i]=valuesautoregression[i];
-   }
-   /*If AR model timestep: compute variable values on vertices from AR*/
-   else{
-      for(int v=0;v<numvertices;v++){
-
-         /*Compute autoregressive term*/
-         IssmDouble autoregressionterm=0.;
-         for(int ii=0;ii<arorder;ii++) autoregressionterm += lagcoefs_basin[ii]*(valuesautoregression[v+ii*numvertices]-(termconstant_basin+trend_basin*(telapsed_ar-(ii+1)*tstep_ar)));
-
-         /*Stochastic variable value*/
-         varlist[v] = termconstant_basin+trend_basin*telapsed_ar+autoregressionterm+noiseterm;
-      }
-      /*Update autoregression values*/
-      IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
-      /*Assign newest values and shift older values*/
-      for(int i=0;i<numvertices;i++) temparray[i] = varlist[i];
-      for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = valuesautoregression[i];
-      this->inputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder);
-      xDelete<IssmDouble>(temparray);
-   }
-
-   /*Add input to element*/
-   this->AddInput(outenum_type,varlist,P1Enum);
-
-   /*Cleanup*/
-   xDelete<IssmDouble>(lagcoefs_basin);
-   xDelete<IssmDouble>(varlist);
-   xDelete<IssmDouble>(valuesautoregression);
-}/*}}}*/
+}
+
+/*}}}*/
 void       Element::BasinLinearFloatingiceMeltingRate(IssmDouble* deepwaterel,IssmDouble* upperwatermelt,IssmDouble* upperwaterel,IssmDouble* perturbation){/*{{{*/
 
@@ -3889 +3827 @@
          /* 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<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -68 +68 @@
                 /*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);

issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp

@@ -231 +231 @@
    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 +247 @@
 
         /*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 +269 @@
       xDelete<int>(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 +274 @@
       Element* element = xDynamicCast<Element*>(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<IssmDouble>(termconstant);
-        xDelete<IssmDouble>(trend);
         xDelete<IssmDouble>(arlagcoefs);
         xDelete<IssmDouble>(malagcoefs);
+        xDelete<IssmDouble>(polyparams);
+        xDelete<IssmDouble>(datebreaks);
         xDelete<IssmDouble>(deepwaterel);
         xDelete<IssmDouble>(upperwaterel);

issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp

@@ -50 +50 @@
         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 +79 @@
                 xDelete<int>(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 +84 @@
       Element* element = xDynamicCast<Element*>(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 +90 @@
 
    /*Cleanup*/
-   xDelete<IssmDouble>(termconstant);
-   xDelete<IssmDouble>(trend);
    xDelete<IssmDouble>(arlagcoefs);
    xDelete<IssmDouble>(malagcoefs);
    xDelete<IssmDouble>(monthlyeff);
+   xDelete<IssmDouble>(polyparams);
+   xDelete<IssmDouble>(datebreaks);
 }/*}}}*/

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -259 +259 @@
                         /*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<IssmDouble>(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<IssmDouble>(transparam);
+         iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.polynomialparams");
+         parameters->AddObject(new DoubleMatParam(BasalforcingsARMApolyparamsEnum,transparam,M,N));
          xDelete<IssmDouble>(transparam);
          iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.arlag_coefs");
@@ -435 +436 @@
          /*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<IssmDouble>(transparam);
-         iomodel->FetchData(&transparam,&M,&N,"md.smb.trend");
-         parameters->AddObject(new DoubleVecParam(SmbARMAtrendEnum,transparam,N));
-         xDelete<IssmDouble>(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<IssmDouble>(transparam);
+         iomodel->FetchData(&transparam,&M,&N,"md.smb.polynomialparams");
+         parameters->AddObject(new DoubleMatParam(SmbARMApolyparamsEnum,transparam,M,N));
+         xDelete<IssmDouble>(transparam);
+                        iomodel->FetchData(&transparam,&M,&N,"md.smb.arlag_coefs");
          parameters->AddObject(new DoubleMatParam(SmbARMAarlagcoefsEnum,transparam,M,N));
          xDelete<IssmDouble>(transparam);

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

@@ -168 +168 @@
    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 +202 @@
       xDelete<int>(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 +207 @@
       Element* element = xDynamicCast<Element*>(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 +213 @@
 
    /*Cleanup*/
-   xDelete<IssmDouble>(termconstant);
-   xDelete<IssmDouble>(trend);
    xDelete<IssmDouble>(arlagcoefs);
    xDelete<IssmDouble>(malagcoefs);
+        xDelete<IssmDouble>(polyparams);
+   xDelete<IssmDouble>(datebreaks);
    xDelete<IssmDouble>(lapserates);
    xDelete<IssmDouble>(elevbins);

issm/trunk-jpl/src/c/shared/Enum/Enum.vim

@@ -67 +67 @@
 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 +84 @@
 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 +214 @@
 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 +320 @@
 syn keyword cConstant LoveMinIntegrationStepsEnum
 syn keyword cConstant LoveMaxIntegrationdrEnum
+syn keyword cConstant LoveIntegrationSchemeEnum
 syn keyword cConstant LoveKernelsEnum
 syn keyword cConstant LoveMu0Enum
@@ -493 +496 @@
 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 +533 @@
 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 +943 @@
 syn keyword cConstant SealevelUNorthEsaEnum
 syn keyword cConstant SealevelchangeIndicesEnum
+syn keyword cConstant SealevelchangeConvolutionVerticesEnum
 syn keyword cConstant SealevelchangeAlphaIndexEnum
 syn keyword cConstant SealevelchangeAzimuthIndexEnum
@@ -957 +962 @@
 syn keyword cConstant SealevelchangeViscousNEnum
 syn keyword cConstant SealevelchangeViscousEEnum
+syn keyword cConstant CouplingTransferCountEnum
 syn keyword cConstant SedimentHeadEnum
 syn keyword cConstant SedimentHeadOldEnum
@@ -1680 +1686 @@
 syn keyword cType Cfsurfacesquare
 syn keyword cType Channel
+syn keyword cType classes
 syn keyword cType Constraint
 syn keyword cType Constraints
@@ -1686 +1693 @@
 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 +1707 @@
 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 +1715 @@
 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 +1737 @@
 syn keyword cType IssmDirectApplicInterface
 syn keyword cType IssmParallelDirectApplicInterface
+syn keyword cType krigingobjects
 syn keyword cType Load
 syn keyword cType Loads
@@ -1739 +1748 @@
 syn keyword cType Matice
 syn keyword cType Matlitho
+syn keyword cType matrixobjects
 syn keyword cType MatrixParam
 syn keyword cType Misfit
@@ -1751 +1761 @@
 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 +1776 @@
 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 +1801 @@
 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 +1811 @@
 syn keyword cType BalancevelocityAnalysis
 syn keyword cType DamageEvolutionAnalysis
+syn keyword cType DebrisAnalysis
 syn keyword cType DepthAverageAnalysis
 syn keyword cType EnthalpyAnalysis

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -61 +61 @@
         BalancethicknessStabilizationEnum,
         BarystaticContributionsEnum,
-        BasalforcingsARMAInitialTimeEnum,
         BasalforcingsARMATimestepEnum,
         BasalforcingsARMAarOrderEnum,
         BasalforcingsARMAmaOrderEnum,
-        BasalforcingsARMAconstEnum,
-        BasalforcingsARMAtrendEnum,
         BasalforcingsBottomplumedepthEnum,
         BasalforcingsCrustthicknessEnum,
@@ -81 +78 @@
         BasalforcingsIsmip6TfDepthsEnum,
         BasalforcingsLinearNumBasinsEnum,
-        BasalforcingsLowercrustheatEnum,
+        BasalforcingsLinearNumBreaksEnum,
+        BasalforcingsLinearNumParamsEnum,
         BasalforcingsMantleconductivityEnum,
         BasalforcingsNusseltEnum,
         BasalforcingsARMAarlagcoefsEnum,
+        BasalforcingsARMAdatebreaksEnum,
         BasalforcingsARMAmalagcoefsEnum,
+        BasalforcingsARMApolyparamsEnum,
         BasalforcingsIsThermalForcingEnum,
+        BasalforcingsLowercrustheatEnum,
         BasalforcingsPicoAverageOverturningEnum,
         BasalforcingsPicoAverageSalinityEnum,
@@ -207 +208 @@
         FrictionVoidRatioEnum,
         FrontalForcingsBasinIcefrontAreaEnum,
-        FrontalForcingsARMAInitialTimeEnum,
    FrontalForcingsARMATimestepEnum,
    FrontalForcingsARMAarOrderEnum,
    FrontalForcingsARMAmaOrderEnum,
-        FrontalForcingsARMAconstEnum,
-   FrontalForcingsARMAtrendEnum,
+   FrontalForcingsARMAdatebreaksEnum,
+   FrontalForcingsARMApolyparamsEnum,
         FrontalForcingsNumberofBasinsEnum,
+        FrontalForcingsNumberofBreaksEnum,
+        FrontalForcingsNumberofParamsEnum,
         FrontalForcingsParamEnum,
    FrontalForcingsARMAarlagcoefsEnum,
@@ -488 +490 @@
         SmbAccurefEnum,
         SmbAdThreshEnum,
-        SmbARMAInitialTimeEnum,
         SmbARMATimestepEnum,
    SmbARMAarOrderEnum,
    SmbARMAmaOrderEnum,
         SmbAveragingEnum,
-        SmbARMAconstEnum,
-   SmbARMAtrendEnum,
         SmbDesfacEnum,
         SmbDpermilEnum,
@@ -528 +527 @@
    SmbLapseRatesEnum,
         SmbNumBasinsEnum,
+        SmbNumBreaksEnum,
         SmbNumElevationBinsEnum,
+        SmbNumParamsEnum,
         SmbNumRequestedOutputsEnum,
         SmbPfacEnum,
         SmbARMAarlagcoefsEnum,
+        SmbARMAdatebreaksEnum,
         SmbARMAmalagcoefsEnum,
+        SmbARMApolyparamsEnum,
         SmbRdlEnum,
         SmbRefElevationEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -69 +69 @@
                 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 +86 @@
                 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 +216 @@
                 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 +498 @@
                 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 +535 @@
                 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";

issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim

@@ -60 +60 @@
 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 +77 @@
 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 +207 @@
 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 +313 @@
 syn keyword juliaConstC LoveMinIntegrationStepsEnum
 syn keyword juliaConstC LoveMaxIntegrationdrEnum
+syn keyword juliaConstC LoveIntegrationSchemeEnum
 syn keyword juliaConstC LoveKernelsEnum
 syn keyword juliaConstC LoveMu0Enum
@@ -486 +489 @@
 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 +526 @@
 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 +936 @@
 syn keyword juliaConstC SealevelUNorthEsaEnum
 syn keyword juliaConstC SealevelchangeIndicesEnum
+syn keyword juliaConstC SealevelchangeConvolutionVerticesEnum
 syn keyword juliaConstC SealevelchangeAlphaIndexEnum
 syn keyword juliaConstC SealevelchangeAzimuthIndexEnum
@@ -950 +955 @@
 syn keyword juliaConstC SealevelchangeViscousNEnum
 syn keyword juliaConstC SealevelchangeViscousEEnum
+syn keyword juliaConstC CouplingTransferCountEnum
 syn keyword juliaConstC SedimentHeadEnum
 syn keyword juliaConstC SedimentHeadOldEnum

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -69 +69 @@
               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 +86 @@
               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 +137 @@
               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 +219 @@
               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 +260 @@
               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 +383 @@
               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 +506 @@
               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 +547 @@
               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 +629 @@
               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 +752 @@
               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 +875 @@
               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 +998 @@
               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 +1121 @@
               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 +1244 @@
               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 +1367 @@
               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 +1490 @@
               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 +1613 @@
               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;

issm/trunk-jpl/src/m/classes/SMBarma.m

@@ -7 +7 @@
         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 +14 @@
                 ma_order          = 0;
                 malag_coefs       = NaN;
+                polynomialparams  = NaN;
+                datebreaks        = NaN;
                 basin_id          = NaN;
                 lapserates        = NaN;
@@ -53 +54 @@
                                 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 +74 @@
 
                         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 +129 @@
                         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 +153 @@
 
                         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 +183 @@
                         [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');

issm/trunk-jpl/src/m/classes/SMBarma.py

@@ -16 +16 @@
     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 +41 @@
         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 +83 @@
             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 +97 @@
     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 +161 @@
     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 +220 @@
         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')

issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m

@@ -7 +7 @@
         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 +55 @@
                         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 +94 @@
                 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 +111 @@
                 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 +155 @@
                                 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 % }}}

issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py

@@ -16 +16 @@
     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 +39 @@
         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 +66 @@
             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 +109 @@
     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')
     # }}}

issm/trunk-jpl/src/m/classes/linearbasalforcingsarma.m

@@ -8 +8 @@
         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 +60 @@
             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 +77 @@
 
                         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 +88 @@
                                 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 +126 @@
                         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 +146 @@
 
                         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 +196 @@
                         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);

issm/trunk-jpl/src/m/classes/linearbasalforcingsarma.py

@@ -15 +15 @@
     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 +43 @@
         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 +81 @@
             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 +95 @@
     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 +107 @@
                 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 +122 @@
             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 +145 @@
     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 +185 @@
         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)

issm/trunk-jpl/test/NightlyRun/test257.m

@@ -33 +33 @@
 
 %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 +71 @@
 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 +85 @@
         (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),...
         };

issm/trunk-jpl/test/NightlyRun/test257.py

@@ -38 +38 @@
 
 # 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 +81 @@
     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 +96 @@
     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
 ]

issm/trunk-jpl/test/NightlyRun/test543.m

@@ -44 +44 @@
 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;

issm/trunk-jpl/test/NightlyRun/test543.py

@@ -49 +49 @@
 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

issm/trunk-jpl/test/NightlyRun/test544.m

@@ -24 +24 @@
 
 %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 +50 @@
 
 % 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;

issm/trunk-jpl/test/NightlyRun/test544.py

@@ -34 +34 @@
 
 #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 +60 @@
 
 #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 +75 @@
 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