Changeset 27246


Ignore:
Timestamp:
08/29/22 08:52:51 (3 years ago)
Author:
vverjans
Message:

CHG: clean-up of autoregression method

Location:
issm/trunk-jpl
Files:
30 edited

Legend:

Unmodified
Added
Removed

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

    r27226 r27246  
    243243         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_initialtime",FrontalForcingsAutoregressionInitialTimeEnum));
    244244         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_timestep",FrontalForcingsAutoregressionTimestepEnum));
    245                         iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.beta0");
    246          parameters->AddObject(new DoubleVecParam(FrontalForcingsBeta0Enum,transparam,N));
     245                        iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.const");
     246         parameters->AddObject(new DoubleVecParam(FrontalForcingsautoregressionconstEnum,transparam,N));
    247247         xDelete<IssmDouble>(transparam);
    248          iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.beta1");
    249          parameters->AddObject(new DoubleVecParam(FrontalForcingsBeta1Enum,transparam,N));
     248         iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.trend");
     249         parameters->AddObject(new DoubleVecParam(FrontalForcingsautoregressiontrendEnum,transparam,N));
    250250         xDelete<IssmDouble>(transparam);
    251          iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.phi");
    252          parameters->AddObject(new DoubleMatParam(FrontalForcingsPhiEnum,transparam,M,N));
     251         iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.arlag_coefs");
     252         parameters->AddObject(new DoubleMatParam(FrontalForcingsarlagcoefsEnum,transparam,M,N));
    253253         xDelete<IssmDouble>(transparam);
    254254                        /*Do not break here, generic FrontalForcingsRignot parameters still to be retrieved*/

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

    r27229 r27246  
    6161        return false;
    6262}/*}}}*/
    63 void       Element::AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type){/*{{{*/
    64 
    65 const int numvertices = this->GetNumberOfVertices();
    66    int         basinid;
    67    IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin;
    68    IssmDouble* phi_basin   = xNew<IssmDouble>(arorder);
    69    IssmDouble* varspin     = xNewZeroInit<IssmDouble>(numvertices*arorder);
    70 
    71    /*Get Basin ID and Basin coefficients*/
    72    if(enum_type==SMBautoregressionEnum)                               this->GetInputValue(&basinid,SmbBasinsIdEnum);
    73    if(enum_type==FrontalForcingsRignotAutoregressionEnum)             this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum);
    74    if(enum_type==BasalforcingsDeepwaterMeltingRateAutoregressionEnum) this->GetInputValue(&basinid,BasalforcingsLinearBasinIdEnum);
    75    for(int i=0;i<arorder;i++) phi_basin[i] = phi[basinid*arorder+i];
    76    beta0_basin   = beta0[basinid];
    77    beta1_basin   = beta1[basinid];
    78 
    79    /*Loop over number of spin-up steps for all vertices*/
    80    for(int j=0;j<nspin;j++){
    81       tspin = starttime-((nspin-j)*tstep_ar);
    82       IssmDouble* oldvarspin = xNew<IssmDouble>(numvertices*arorder);
    83       for(int i=0;i<numvertices*arorder;i++) oldvarspin[i]=varspin[i];
    84 
    85       for(int v=0;v<numvertices;v++){
    86          IssmDouble autoregressionterm = 0.;
    87          for(int i=0;i<arorder;i++) autoregressionterm += phi_basin[i]*varspin[v+i*numvertices];
    88          varspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm;
    89       }
    90 
    91       /*Adjust older values in varspin*/
    92       for(int i=0;i<(arorder-1)*numvertices;i++) varspin[i+numvertices]=oldvarspin[i];
    93 
    94       xDelete<IssmDouble>(oldvarspin);
    95    }
    96    if(enum_type==SMBautoregressionEnum)                               this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder);
    97    if(enum_type==FrontalForcingsRignotAutoregressionEnum)             this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder);
    98    if(enum_type==BasalforcingsDeepwaterMeltingRateAutoregressionEnum) this->inputs->SetArrayInput(BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder);
    99 
    100    /*Cleanup and return*/
    101    xDelete<IssmDouble>(varspin);
    102    xDelete<IssmDouble>(phi_basin);     
    103 }/*}}}*/
    104 void       Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type){/*{{{*/
     63void       Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble tstep_ar,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* lagcoefs,bool isfieldstochastic,int enum_type){/*{{{*/
    10564
    10665   const int numvertices = this->GetNumberOfVertices();
    10766   int         basinid,M,N,arenum_type,basinenum_type,noiseenum_type,outenum_type;
    108    IssmDouble  beta0_basin,beta1_basin,noiseterm;
    109    IssmDouble* phi_basin  = xNew<IssmDouble>(arorder);
    110    IssmDouble* varlist     = xNew<IssmDouble>(numvertices);
     67   IssmDouble  time,dt,starttime,termconstant_basin,trend_basin,noiseterm;
     68   IssmDouble* lagcoefs_basin  = xNew<IssmDouble>(arorder);
     69   IssmDouble* varlist         = xNew<IssmDouble>(numvertices);
    11170   IssmDouble* valuesautoregression = NULL;
    11271   Input*      noiseterm_input      = NULL;
    … …  
    13493   }
    13594
     95        /*Get time parameters*/
     96   this->parameters->FindParam(&time,TimeEnum);
     97   this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     98   this->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
     99
     100   /*Get basin coefficients*/
     101   this->GetInputValue(&basinid,basinenum_type);
     102   for(int ii=0;ii<arorder;ii++) lagcoefs_basin[ii] = lagcoefs[basinid*arorder+ii];
     103   termconstant_basin   = termconstant[basinid];
     104   trend_basin   = trend[basinid];
     105
     106        /*Initialze autoregressive values at first time step*/
     107        if(time<=starttime+dt){
     108                IssmDouble* initvaluesautoregression = xNewZeroInit<IssmDouble>(numvertices*arorder);
     109                for(int i=0;i<numvertices*arorder;i++) initvaluesautoregression[i]=termconstant_basin;
     110      this->inputs->SetArrayInput(arenum_type,this->lid,initvaluesautoregression,numvertices*arorder);
     111      xDelete<IssmDouble>(initvaluesautoregression);
     112        }
     113
    136114   /*Get noise and autoregressive terms*/
    137    this->GetInputValue(&basinid,basinenum_type);
    138115   if(isfieldstochastic){
    139116      noiseterm_input = this->GetInput(noiseenum_type);
    … …  
    145122   this->inputs->GetArray(arenum_type,this->lid,&valuesautoregression,&M);
    146123
    147    /*Get basin coefficients*/
    148    for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii];
    149    beta0_basin   = beta0[basinid];
    150    beta1_basin   = beta1[basinid];
    151 
    152124        /*If not AR model timestep: take the old values of variable*/
    153125   if(isstepforar==false){
    … …  
    160132         /*Compute autoregressive term*/
    161133         IssmDouble autoregressionterm=0.;
    162          for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*valuesautoregression[v+ii*numvertices];
     134         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)));
    163135
    164136         /*Stochastic variable value*/
    165          varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noiseterm;
     137         varlist[v] = termconstant_basin+trend_basin*telapsed_ar+autoregressionterm+noiseterm;
    166138      }
    167139      /*Update autoregression values*/
    … …  
    178150
    179151   /*Cleanup*/
    180    xDelete<IssmDouble>(phi_basin);
     152   xDelete<IssmDouble>(lagcoefs_basin);
    181153   xDelete<IssmDouble>(varlist);
    182154   xDelete<IssmDouble>(valuesautoregression);

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

    r27131 r27246  
    6868                /*bool               AnyActive(void);*/
    6969                bool               AnyFSet(void);
    70                 void                                     AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type);
    71       void               Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type);
     70      void               Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble tstep_ar,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* lagcoefs,bool isfieldstochastic,int enum_type);
    7271                void               BasinLinearFloatingiceMeltingRate(IssmDouble* deepwaterel,IssmDouble* upperwatermelt,IssmDouble* upperwaterel,IssmDouble* perturbation);
    7372                void               CalvingSetZeroRate(void);

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

    r26990 r27246  
    210210}
    211211/*}}}*/
    212 void AutoregressionLinearFloatingiceMeltingRateInitx(FemModel* femmodel){/*{{{*/
    213 
    214         /*Initialization step of AutoregressionLinearFloatingiceMeltingRatex*/
    215    int M,N,Nphi,arorder,numbasins,my_rank;
    216    IssmDouble starttime,tstep_ar,tinit_ar;
    217    femmodel->parameters->FindParam(&numbasins,BasalforcingsLinearNumBasinsEnum);
    218    femmodel->parameters->FindParam(&arorder,BasalforcingsAutoregressiveOrderEnum);
    219    IssmDouble* beta0    = NULL;
    220    IssmDouble* beta1    = NULL;
    221    IssmDouble* phi      = NULL;
    222    femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    223    femmodel->parameters->FindParam(&tstep_ar,BasalforcingsAutoregressionTimestepEnum);
    224    femmodel->parameters->FindParam(&tinit_ar,BasalforcingsAutoregressionInitialTimeEnum);
    225    femmodel->parameters->FindParam(&beta0,&M,BasalforcingsBeta0Enum);    _assert_(M==numbasins);
    226    femmodel->parameters->FindParam(&beta1,&M,BasalforcingsBeta1Enum);    _assert_(M==numbasins);
    227    femmodel->parameters->FindParam(&phi,&M,&Nphi,BasalforcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
    228        
    229         /*AR model spin-up with 0 noise to initialize BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum (688 = log(0.001)/log(0.99): decaying time of inluence of phi[0]=0.99 to 0.001 of beta_0*/
    230    int nspin = 688;
    231    for(Object* &object:femmodel->elements->objects){
    232       Element* element = xDynamicCast<Element*>(object); //generate element object
    233       element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,BasalforcingsDeepwaterMeltingRateAutoregressionEnum);
    234         }
    235         /*Cleanup*/
    236    xDelete<IssmDouble>(beta0);
    237    xDelete<IssmDouble>(beta1);
    238    xDelete<IssmDouble>(phi);
    239 }/*}}}*/
    240212void AutoregressionLinearFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/
    241213
    … …  
    247219   femmodel->parameters->FindParam(&tstep_ar,BasalforcingsAutoregressionTimestepEnum);
    248220
    249    /*Initialize module at first time step*/
    250    if(time<=starttime+dt){AutoregressionLinearFloatingiceMeltingRateInitx(femmodel);}
    251221   /*Determine if this is a time step for the AR model*/
    252222   bool isstepforar = false;
    … …  
    261231   bool isstochastic;
    262232   bool isdeepmeltingstochastic = false;
    263    int M,N,Nphi,arorder,numbasins,my_rank;
     233   int M,N,Nlags,arorder,numbasins,my_rank;
    264234   femmodel->parameters->FindParam(&numbasins,BasalforcingsLinearNumBasinsEnum);
    265235        femmodel->parameters->FindParam(&arorder,BasalforcingsAutoregressiveOrderEnum);
    266236   IssmDouble tinit_ar;
    267    IssmDouble* beta0          = NULL;
    268    IssmDouble* beta1          = NULL;
    269    IssmDouble* phi            = NULL;
     237   IssmDouble* termconstant   = NULL;
     238   IssmDouble* trend          = NULL;
     239   IssmDouble* lagcoefs       = NULL;
    270240   IssmDouble* deepwaterel    = NULL;
    271241   IssmDouble* upperwaterel   = NULL;
    … …  
    275245        /*Get autoregressive parameters*/
    276246   femmodel->parameters->FindParam(&tinit_ar,BasalforcingsAutoregressionInitialTimeEnum);
    277    femmodel->parameters->FindParam(&beta0,&M,BasalforcingsBeta0Enum);               _assert_(M==numbasins);
    278    femmodel->parameters->FindParam(&beta1,&M,BasalforcingsBeta1Enum);               _assert_(M==numbasins);
    279    femmodel->parameters->FindParam(&phi,&M,&Nphi,BasalforcingsPhiEnum);             _assert_(M==numbasins); _assert_(Nphi==arorder);
     247   femmodel->parameters->FindParam(&termconstant,&M,BasalforcingsautoregressionconstEnum);               _assert_(M==numbasins);
     248   femmodel->parameters->FindParam(&trend,&M,BasalforcingsautoregressiontrendEnum);               _assert_(M==numbasins);
     249   femmodel->parameters->FindParam(&lagcoefs,&M,&Nlags,BasalforcingsarlagcoefsEnum);             _assert_(M==numbasins); _assert_(Nlags==arorder);
    280250
    281251        /*Get basin-specific parameters*/
    … …  
    303273      Element* element = xDynamicCast<Element*>(object);
    304274      /*Compute autoregression*/
    305       element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,isdeepmeltingstochastic,BasalforcingsDeepwaterMeltingRateAutoregressionEnum);
     275      element->Autoregression(isstepforar,arorder,telapsed_ar,tstep_ar,termconstant,trend,lagcoefs,isdeepmeltingstochastic,BasalforcingsDeepwaterMeltingRateAutoregressionEnum);
    306276                element->BasinLinearFloatingiceMeltingRate(deepwaterel,upperwatermelt,upperwaterel,perturbation);
    307277        }
    308278
    309279        /*Cleanup*/
    310         xDelete<IssmDouble>(beta0);
    311         xDelete<IssmDouble>(beta1);
    312         xDelete<IssmDouble>(phi);
     280        xDelete<IssmDouble>(termconstant);
     281        xDelete<IssmDouble>(trend);
     282        xDelete<IssmDouble>(lagcoefs);
    313283        xDelete<IssmDouble>(deepwaterel);
    314284        xDelete<IssmDouble>(upperwaterel);

  • issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h

    r26836 r27246  
    1616void FloatingiceMeltingRateIsmip6x(FemModel* femmodel);
    1717void BeckmannGoosseFloatingiceMeltingRatex(FemModel* femmodel);
    18 void AutoregressionLinearFloatingiceMeltingRateInitx(FemModel* femmodel);
    1918void AutoregressionLinearFloatingiceMeltingRatex(FemModel* femmodel);
    2019

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

    r26751 r27246  
    2929        }
    3030}/*}}}*/
    31 void ThermalforcingautoregressionInitx(FemModel* femmodel){/*{{{*/
    32 
    33    /*Initialization step of Thermalforcingautoregressionx*/
    34    int M,N,Nphi,arorder,numbasins,my_rank;
    35    IssmDouble starttime,tstep_ar,tinit_ar;
    36    femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
    37    femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
    38    IssmDouble* beta0    = NULL;
    39    IssmDouble* beta1    = NULL;
    40    IssmDouble* phi      = NULL;
    41    femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    42    femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum);
    43    femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
    44    femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum);    _assert_(M==numbasins);
    45    femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum);    _assert_(M==numbasins);
    46    femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
    47 
    48    /*AR model spin-up with 0 noise to initialize ThermalforcingValuesAutoregressionEnum (688 = log(0.001)/log(0.99): decaying time of inluence of phi[0]=0.99 to 0.001 of beta_0*/
    49         int nspin = 688;
    50    for(Object* &object:femmodel->elements->objects){
    51       Element* element      = xDynamicCast<Element*>(object); //generate element object
    52       element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,FrontalForcingsRignotAutoregressionEnum);
    53    }
    54 
    55    /*Cleanup*/
    56    xDelete<IssmDouble>(beta0);
    57    xDelete<IssmDouble>(beta1);
    58    xDelete<IssmDouble>(phi);
    59 }/*}}}*/
    6031void Thermalforcingautoregressionx(FemModel* femmodel){/*{{{*/
    6132
    … …  
    6738   femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum);
    6839
    69    /*Initialize module at first time step*/
    70    if(time<=starttime+dt){ThermalforcingautoregressionInitx(femmodel);}
    7140   /*Determine if this is a time step for the AR model*/
    7241   bool isstepforar = false;
    … …  
    8150        bool isstochastic;
    8251   bool istfstochastic = false;
    83         int M,N,Nphi,arorder,numbasins,my_rank;
     52        int M,N,Nlagcoefs,arorder,numbasins,my_rank;
    8453   femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
    8554   femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
    8655   IssmDouble tinit_ar;
    87    IssmDouble* beta0      = NULL;
    88    IssmDouble* beta1      = NULL;
    89    IssmDouble* phi        = NULL;
     56   IssmDouble* termconstant  = NULL;
     57   IssmDouble* trend         = NULL;
     58   IssmDouble* lagcoefs      = NULL;
    9059
    9160        femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
    92    femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum);    _assert_(M==numbasins);
    93    femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum);    _assert_(M==numbasins);
    94    femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
     61   femmodel->parameters->FindParam(&termconstant,&M,FrontalForcingsautoregressionconstEnum);    _assert_(M==numbasins);
     62   femmodel->parameters->FindParam(&trend,&M,FrontalForcingsautoregressiontrendEnum);    _assert_(M==numbasins);
     63   femmodel->parameters->FindParam(&lagcoefs,&M,&Nlagcoefs,FrontalForcingsarlagcoefsEnum);  _assert_(M==numbasins); _assert_(Nlagcoefs==arorder);
    9564
    9665        femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
    … …  
    11180   for(Object* &object:femmodel->elements->objects){
    11281      Element* element = xDynamicCast<Element*>(object);
    113       element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum);
     82      element->Autoregression(isstepforar,arorder,telapsed_ar,tstep_ar,termconstant,trend,lagcoefs,istfstochastic,FrontalForcingsRignotAutoregressionEnum);
    11483   }
    11584
    11685   /*Cleanup*/
    117    xDelete<IssmDouble>(beta0);
    118    xDelete<IssmDouble>(beta1);
    119    xDelete<IssmDouble>(phi);
     86   xDelete<IssmDouble>(termconstant);
     87   xDelete<IssmDouble>(trend);
     88   xDelete<IssmDouble>(lagcoefs);
    12089}/*}}}*/

  • issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h

    r26495 r27246  
    77/* local prototypes: */
    88void FrontalForcingsx(FemModel* femmodel);
    9 void ThermalforcingautoregressionInitx(FemModel* femmodel);
    109void Thermalforcingautoregressionx(FemModel* femmodel);
    1110

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

    r27181 r27246  
    260260         parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.ar_initialtime",BasalforcingsAutoregressionInitialTimeEnum));
    261261         parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.ar_timestep",BasalforcingsAutoregressionTimestepEnum));
    262                         iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.beta0");
    263          parameters->AddObject(new DoubleVecParam(BasalforcingsBeta0Enum,transparam,N));
    264          xDelete<IssmDouble>(transparam);
    265          iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.beta1");
    266          parameters->AddObject(new DoubleVecParam(BasalforcingsBeta1Enum,transparam,N));
    267          xDelete<IssmDouble>(transparam);
    268          iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.phi");
    269          parameters->AddObject(new DoubleMatParam(BasalforcingsPhiEnum,transparam,M,N));
     262                        iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.const");
     263         parameters->AddObject(new DoubleVecParam(BasalforcingsautoregressionconstEnum,transparam,N));
     264         xDelete<IssmDouble>(transparam);
     265         iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.trend");
     266         parameters->AddObject(new DoubleVecParam(BasalforcingsautoregressiontrendEnum,transparam,N));
     267         xDelete<IssmDouble>(transparam);
     268         iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.arlag_coefs");
     269         parameters->AddObject(new DoubleMatParam(BasalforcingsarlagcoefsEnum,transparam,M,N));
    270270         xDelete<IssmDouble>(transparam);
    271271                        iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.upperwater_melting_rate");
    … …  
    432432         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_timestep",SmbAutoregressionTimestepEnum));
    433433         parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_bins",SmbNumElevationBinsEnum));
    434          iomodel->FetchData(&transparam,&M,&N,"md.smb.beta0");
    435          parameters->AddObject(new DoubleVecParam(SmbBeta0Enum,transparam,N));
    436          xDelete<IssmDouble>(transparam);
    437          iomodel->FetchData(&transparam,&M,&N,"md.smb.beta1");
    438          parameters->AddObject(new DoubleVecParam(SmbBeta1Enum,transparam,N));
    439          xDelete<IssmDouble>(transparam);
    440          iomodel->FetchData(&transparam,&M,&N,"md.smb.phi");
    441          parameters->AddObject(new DoubleMatParam(SmbPhiEnum,transparam,M,N));
     434         iomodel->FetchData(&transparam,&M,&N,"md.smb.const");
     435         parameters->AddObject(new DoubleVecParam(SmbautoregressionconstEnum,transparam,N));
     436         xDelete<IssmDouble>(transparam);
     437         iomodel->FetchData(&transparam,&M,&N,"md.smb.trend");
     438         parameters->AddObject(new DoubleVecParam(SmbautoregressiontrendEnum,transparam,N));
     439         xDelete<IssmDouble>(transparam);
     440         iomodel->FetchData(&transparam,&M,&N,"md.smb.arlag_coefs");
     441         parameters->AddObject(new DoubleMatParam(SmbarlagcoefsEnum,transparam,M,N));
    442442         xDelete<IssmDouble>(transparam);
    443443         iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserates");

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

    r26947 r27246  
    147147
    148148}/*}}}*/
    149 void SmbautoregressionInitx(FemModel* femmodel){/*{{{*/
    150        
    151         /*Initialization step of Smbautoregressionx*/
    152         int M,N,Nphi,arorder,numbasins,my_rank;
    153         IssmDouble starttime,tstep_ar,tinit_ar;
    154         femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
    155         femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
    156         IssmDouble* beta0    = NULL;
    157         IssmDouble* beta1    = NULL;
    158         IssmDouble* phi      = NULL;
    159         femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    160    femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
    161         femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
    162         femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
    163         femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
    164         femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
    165        
    166         /*AR model spin-up with 0 noise to initialize SmbValuesAutoregressionEnum (688 = log(0.001)/log(0.99): decaying time of inluence of phi[0]=0.99 to 0.001 of beta_0*/
    167         int nspin = 688;
    168         for(Object* &object:femmodel->elements->objects){
    169       Element* element      = xDynamicCast<Element*>(object); //generate element object
    170                 element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,SMBautoregressionEnum);
    171         }
    172         /*Cleanup*/
    173         xDelete<IssmDouble>(beta0);
    174         xDelete<IssmDouble>(beta1);
    175         xDelete<IssmDouble>(phi);
    176 }/*}}}*/
    177149void Smbautoregressionx(FemModel* femmodel){/*{{{*/
    178150
    … …  
    184156   femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
    185157
    186    /*Initialize module at first time step*/
    187    if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
    188158   /*Determine if this is a time step for the AR model*/
    189159   bool isstepforar = false;
    … …  
    198168   bool isstochastic;
    199169   bool issmbstochastic = false;
    200    int M,N,Nphi,arorder,numbasins,numelevbins,my_rank;
     170   int M,N,Nlagcoefs,arorder,numbasins,numelevbins,my_rank;
    201171   femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
    202172   femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
    203173   femmodel->parameters->FindParam(&numelevbins,SmbNumElevationBinsEnum);
    204174   IssmDouble tinit_ar;
    205    IssmDouble* beta0         = NULL;
    206    IssmDouble* beta1         = NULL;
    207    IssmDouble* phi           = NULL;
     175   IssmDouble* termconstant  = NULL;
     176   IssmDouble* trend         = NULL;
     177   IssmDouble* lagcoefs      = NULL;
    208178   IssmDouble* lapserates    = NULL;
    209179   IssmDouble* elevbins      = NULL;
    … …  
    211181
    212182   femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
    213    femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);                  _assert_(M==numbasins);
    214    femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);                  _assert_(M==numbasins);
    215    femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);                _assert_(M==numbasins); _assert_(Nphi==arorder);
    216    femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum);     _assert_(M==numbasins); _assert_(N==numelevbins);
    217    femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum);    _assert_(M==numbasins); _assert_(N==numelevbins-1);
    218    femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum);    _assert_(M==numbasins);
     183   femmodel->parameters->FindParam(&termconstant,&M,SmbautoregressionconstEnum);  _assert_(M==numbasins);
     184   femmodel->parameters->FindParam(&trend,&M,SmbautoregressiontrendEnum);         _assert_(M==numbasins);
     185   femmodel->parameters->FindParam(&lagcoefs,&M,&Nlagcoefs,SmbarlagcoefsEnum);    _assert_(M==numbasins); _assert_(Nlagcoefs==arorder);
     186   femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum);          _assert_(M==numbasins); _assert_(N==numelevbins);
     187   femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum);         _assert_(M==numbasins); _assert_(N==numelevbins-1);
     188   femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum);         _assert_(M==numbasins);
    219189
    220190   femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
    … …  
    236206      Element* element = xDynamicCast<Element*>(object);
    237207      /*Compute autoregression*/
    238                 element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum);
     208                element->Autoregression(isstepforar,arorder,telapsed_ar,tstep_ar,termconstant,trend,lagcoefs,issmbstochastic,SMBautoregressionEnum);
    239209                /*Compute lapse rate adjustment*/
    240210                element->LapseRateBasinSMB(numelevbins,lapserates,elevbins,refelevation);
    … …  
    242212
    243213   /*Cleanup*/
    244    xDelete<IssmDouble>(beta0);
    245    xDelete<IssmDouble>(beta1);
    246    xDelete<IssmDouble>(phi);
     214   xDelete<IssmDouble>(termconstant);
     215   xDelete<IssmDouble>(trend);
     216   xDelete<IssmDouble>(lagcoefs);
    247217   xDelete<IssmDouble>(lapserates);
    248218   xDelete<IssmDouble>(elevbins);

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h

    r27240 r27246  
    1313void SmbGradientsx(FemModel* femmodel);
    1414void SmbGradientsElax(FemModel* femmodel);
    15 void SmbautoregressionInitx(FemModel* femmodel);
    1615void Smbautoregressionx(FemModel* femmodel);
    1716void Delta18oParameterizationx(FemModel* femmodel);

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

    r27229 r27246  
    7070syn keyword cConstant BasalforcingsAutoregressionTimestepEnum
    7171syn keyword cConstant BasalforcingsAutoregressiveOrderEnum
    72 syn keyword cConstant BasalforcingsBeta0Enum
    73 syn keyword cConstant BasalforcingsBeta1Enum
     72syn keyword cConstant BasalforcingsautoregressionconstEnum
     73syn keyword cConstant BasalforcingsautoregressiontrendEnum
    7474syn keyword cConstant BasalforcingsBottomplumedepthEnum
    7575syn keyword cConstant BasalforcingsCrustthicknessEnum
    … …  
    8989syn keyword cConstant BasalforcingsMantleconductivityEnum
    9090syn keyword cConstant BasalforcingsNusseltEnum
    91 syn keyword cConstant BasalforcingsPhiEnum
     91syn keyword cConstant BasalforcingsarlagcoefsEnum
    9292syn keyword cConstant BasalforcingsPicoAverageOverturningEnum
    9393syn keyword cConstant BasalforcingsPicoAverageSalinityEnum
    … …  
    199199syn keyword cConstant FrontalForcingsAutoregressionTimestepEnum
    200200syn keyword cConstant FrontalForcingsAutoregressiveOrderEnum
    201 syn keyword cConstant FrontalForcingsBeta0Enum
    202 syn keyword cConstant FrontalForcingsBeta1Enum
     201syn keyword cConstant FrontalForcingsautoregressionconstEnum
     202syn keyword cConstant FrontalForcingsautoregressiontrendEnum
    203203syn keyword cConstant FrontalForcingsNumberofBasinsEnum
    204204syn keyword cConstant FrontalForcingsParamEnum
    205 syn keyword cConstant FrontalForcingsPhiEnum
     205syn keyword cConstant FrontalForcingsarlagcoefsEnum
    206206syn keyword cConstant GrdModelEnum
    207207syn keyword cConstant GroundinglineFrictionInterpolationEnum
    … …  
    478478syn keyword cConstant SmbAutoregressiveOrderEnum
    479479syn keyword cConstant SmbAveragingEnum
    480 syn keyword cConstant SmbBeta0Enum
    481 syn keyword cConstant SmbBeta1Enum
     480syn keyword cConstant SmbautoregressionconstEnum
     481syn keyword cConstant SmbautoregressiontrendEnum
    482482syn keyword cConstant SmbDesfacEnum
    483483syn keyword cConstant SmbDpermilEnum
    … …  
    516516syn keyword cConstant SmbNumRequestedOutputsEnum
    517517syn keyword cConstant SmbPfacEnum
    518 syn keyword cConstant SmbPhiEnum
     518syn keyword cConstant SmbarlagcoefsEnum
    519519syn keyword cConstant SmbRdlEnum
    520520syn keyword cConstant SmbRefElevationEnum

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

    r27229 r27246  
    6464        BasalforcingsAutoregressionTimestepEnum,
    6565        BasalforcingsAutoregressiveOrderEnum,
    66         BasalforcingsBeta0Enum,
    67         BasalforcingsBeta1Enum,
     66        BasalforcingsautoregressionconstEnum,
     67        BasalforcingsautoregressiontrendEnum,
    6868        BasalforcingsBottomplumedepthEnum,
    6969        BasalforcingsCrustthicknessEnum,
    … …  
    8383        BasalforcingsMantleconductivityEnum,
    8484        BasalforcingsNusseltEnum,
    85         BasalforcingsPhiEnum,
     85        BasalforcingsarlagcoefsEnum,
    8686        BasalforcingsPicoAverageOverturningEnum,
    8787        BasalforcingsPicoAverageSalinityEnum,
    … …  
    193193   FrontalForcingsAutoregressionTimestepEnum,
    194194   FrontalForcingsAutoregressiveOrderEnum,
    195         FrontalForcingsBeta0Enum,
    196    FrontalForcingsBeta1Enum,
     195        FrontalForcingsautoregressionconstEnum,
     196   FrontalForcingsautoregressiontrendEnum,
    197197        FrontalForcingsNumberofBasinsEnum,
    198198        FrontalForcingsParamEnum,
    199    FrontalForcingsPhiEnum,
     199   FrontalForcingsarlagcoefsEnum,
    200200        GrdModelEnum,
    201201        GroundinglineFrictionInterpolationEnum,
    … …  
    472472   SmbAutoregressiveOrderEnum,
    473473        SmbAveragingEnum,
    474         SmbBeta0Enum,
    475    SmbBeta1Enum,
     474        SmbautoregressionconstEnum,
     475   SmbautoregressiontrendEnum,
    476476        SmbDesfacEnum,
    477477        SmbDpermilEnum,
    … …  
    510510        SmbNumRequestedOutputsEnum,
    511511        SmbPfacEnum,
    512         SmbPhiEnum,
     512        SmbarlagcoefsEnum,
    513513        SmbRdlEnum,
    514514        SmbRefElevationEnum,

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

    r27229 r27246  
    7272                case BasalforcingsAutoregressionTimestepEnum : return "BasalforcingsAutoregressionTimestep";
    7373                case BasalforcingsAutoregressiveOrderEnum : return "BasalforcingsAutoregressiveOrder";
    74                 case BasalforcingsBeta0Enum : return "BasalforcingsBeta0";
    75                 case BasalforcingsBeta1Enum : return "BasalforcingsBeta1";
     74                case BasalforcingsautoregressionconstEnum : return "Basalforcingsautoregressionconst";
     75                case BasalforcingsautoregressiontrendEnum : return "Basalforcingsautoregressiontrend";
    7676                case BasalforcingsBottomplumedepthEnum : return "BasalforcingsBottomplumedepth";
    7777                case BasalforcingsCrustthicknessEnum : return "BasalforcingsCrustthickness";
    … …  
    9191                case BasalforcingsMantleconductivityEnum : return "BasalforcingsMantleconductivity";
    9292                case BasalforcingsNusseltEnum : return "BasalforcingsNusselt";
    93                 case BasalforcingsPhiEnum : return "BasalforcingsPhi";
     93                case BasalforcingsarlagcoefsEnum : return "Basalforcingsarlagcoefs";
    9494                case BasalforcingsPicoAverageOverturningEnum : return "BasalforcingsPicoAverageOverturning";
    9595                case BasalforcingsPicoAverageSalinityEnum : return "BasalforcingsPicoAverageSalinity";
    … …  
    201201                case FrontalForcingsAutoregressionTimestepEnum : return "FrontalForcingsAutoregressionTimestep";
    202202                case FrontalForcingsAutoregressiveOrderEnum : return "FrontalForcingsAutoregressiveOrder";
    203                 case FrontalForcingsBeta0Enum : return "FrontalForcingsBeta0";
    204                 case FrontalForcingsBeta1Enum : return "FrontalForcingsBeta1";
     203                case FrontalForcingsautoregressionconstEnum : return "FrontalForcingsautoregressionconst";
     204                case FrontalForcingsautoregressiontrendEnum : return "FrontalForcingsautoregressiontrend";
    205205                case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins";
    206206                case FrontalForcingsParamEnum : return "FrontalForcingsParam";
    207                 case FrontalForcingsPhiEnum : return "FrontalForcingsPhi";
     207                case FrontalForcingsarlagcoefsEnum : return "FrontalForcingsarlagcoefs";
    208208                case GrdModelEnum : return "GrdModel";
    209209                case GroundinglineFrictionInterpolationEnum : return "GroundinglineFrictionInterpolation";
    … …  
    480480                case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder";
    481481                case SmbAveragingEnum : return "SmbAveraging";
    482                 case SmbBeta0Enum : return "SmbBeta0";
    483                 case SmbBeta1Enum : return "SmbBeta1";
     482                case SmbautoregressionconstEnum : return "Smbautoregressionconst";
     483                case SmbautoregressiontrendEnum : return "Smbautoregressiontrend";
    484484                case SmbDesfacEnum : return "SmbDesfac";
    485485                case SmbDpermilEnum : return "SmbDpermil";
    … …  
    518518                case SmbNumRequestedOutputsEnum : return "SmbNumRequestedOutputs";
    519519                case SmbPfacEnum : return "SmbPfac";
    520                 case SmbPhiEnum : return "SmbPhi";
     520                case SmbarlagcoefsEnum : return "Smbarlagcoefs";
    521521                case SmbRdlEnum : return "SmbRdl";
    522522                case SmbRefElevationEnum : return "SmbRefElevation";

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

    r27229 r27246  
    6363syn keyword juliaConstC BasalforcingsAutoregressionTimestepEnum
    6464syn keyword juliaConstC BasalforcingsAutoregressiveOrderEnum
    65 syn keyword juliaConstC BasalforcingsBeta0Enum
    66 syn keyword juliaConstC BasalforcingsBeta1Enum
     65syn keyword juliaConstC BasalforcingsautoregressionconstEnum
     66syn keyword juliaConstC BasalforcingsautoregressiontrendEnum
    6767syn keyword juliaConstC BasalforcingsBottomplumedepthEnum
    6868syn keyword juliaConstC BasalforcingsCrustthicknessEnum
    … …  
    8282syn keyword juliaConstC BasalforcingsMantleconductivityEnum
    8383syn keyword juliaConstC BasalforcingsNusseltEnum
    84 syn keyword juliaConstC BasalforcingsPhiEnum
     84syn keyword juliaConstC BasalforcingsarlagcoefsEnum
    8585syn keyword juliaConstC BasalforcingsPicoAverageOverturningEnum
    8686syn keyword juliaConstC BasalforcingsPicoAverageSalinityEnum
    … …  
    192192syn keyword juliaConstC FrontalForcingsAutoregressionTimestepEnum
    193193syn keyword juliaConstC FrontalForcingsAutoregressiveOrderEnum
    194 syn keyword juliaConstC FrontalForcingsBeta0Enum
    195 syn keyword juliaConstC FrontalForcingsBeta1Enum
     194syn keyword juliaConstC FrontalForcingsautoregressionconstEnum
     195syn keyword juliaConstC FrontalForcingsautoregressiontrendEnum
    196196syn keyword juliaConstC FrontalForcingsNumberofBasinsEnum
    197197syn keyword juliaConstC FrontalForcingsParamEnum
    198 syn keyword juliaConstC FrontalForcingsPhiEnum
     198syn keyword juliaConstC FrontalForcingsarlagcoefsEnum
    199199syn keyword juliaConstC GrdModelEnum
    200200syn keyword juliaConstC GroundinglineFrictionInterpolationEnum
    … …  
    471471syn keyword juliaConstC SmbAutoregressiveOrderEnum
    472472syn keyword juliaConstC SmbAveragingEnum
    473 syn keyword juliaConstC SmbBeta0Enum
    474 syn keyword juliaConstC SmbBeta1Enum
     473syn keyword juliaConstC SmbautoregressionconstEnum
     474syn keyword juliaConstC SmbautoregressiontrendEnum
    475475syn keyword juliaConstC SmbDesfacEnum
    476476syn keyword juliaConstC SmbDpermilEnum
    … …  
    509509syn keyword juliaConstC SmbNumRequestedOutputsEnum
    510510syn keyword juliaConstC SmbPfacEnum
    511 syn keyword juliaConstC SmbPhiEnum
     511syn keyword juliaConstC SmbarlagcoefsEnum
    512512syn keyword juliaConstC SmbRdlEnum
    513513syn keyword juliaConstC SmbRefElevationEnum
    … …  
    864864syn keyword juliaConstC SamplingBetaEnum
    865865syn keyword juliaConstC SamplingKappaEnum
    866 syn keyword juliaConstC SamplingPhiEnum
     866syn keyword juliaConstC SamplingarlagcoefsEnum
    867867syn keyword juliaConstC SamplingTauEnum
    868868syn keyword juliaConstC SealevelEnum

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

    r27229 r27246  
    7272              else if (strcmp(name,"BasalforcingsAutoregressionTimestep")==0) return BasalforcingsAutoregressionTimestepEnum;
    7373              else if (strcmp(name,"BasalforcingsAutoregressiveOrder")==0) return BasalforcingsAutoregressiveOrderEnum;
    74               else if (strcmp(name,"BasalforcingsBeta0")==0) return BasalforcingsBeta0Enum;
    75               else if (strcmp(name,"BasalforcingsBeta1")==0) return BasalforcingsBeta1Enum;
     74              else if (strcmp(name,"Basalforcingsautoregressionconst")==0) return BasalforcingsautoregressionconstEnum;
     75              else if (strcmp(name,"Basalforcingsautoregressiontrend")==0) return BasalforcingsautoregressiontrendEnum;
    7676              else if (strcmp(name,"BasalforcingsBottomplumedepth")==0) return BasalforcingsBottomplumedepthEnum;
    7777              else if (strcmp(name,"BasalforcingsCrustthickness")==0) return BasalforcingsCrustthicknessEnum;
    … …  
    9191              else if (strcmp(name,"BasalforcingsMantleconductivity")==0) return BasalforcingsMantleconductivityEnum;
    9292              else if (strcmp(name,"BasalforcingsNusselt")==0) return BasalforcingsNusseltEnum;
    93               else if (strcmp(name,"BasalforcingsPhi")==0) return BasalforcingsPhiEnum;
     93              else if (strcmp(name,"Basalforcingsarlagcoefs")==0) return BasalforcingsarlagcoefsEnum;
    9494              else if (strcmp(name,"BasalforcingsPicoAverageOverturning")==0) return BasalforcingsPicoAverageOverturningEnum;
    9595              else if (strcmp(name,"BasalforcingsPicoAverageSalinity")==0) return BasalforcingsPicoAverageSalinityEnum;
    … …  
    204204              else if (strcmp(name,"FrontalForcingsAutoregressionTimestep")==0) return FrontalForcingsAutoregressionTimestepEnum;
    205205              else if (strcmp(name,"FrontalForcingsAutoregressiveOrder")==0) return FrontalForcingsAutoregressiveOrderEnum;
    206               else if (strcmp(name,"FrontalForcingsBeta0")==0) return FrontalForcingsBeta0Enum;
    207               else if (strcmp(name,"FrontalForcingsBeta1")==0) return FrontalForcingsBeta1Enum;
     206              else if (strcmp(name,"FrontalForcingsautoregressionconst")==0) return FrontalForcingsautoregressionconstEnum;
     207              else if (strcmp(name,"FrontalForcingsautoregressiontrend")==0) return FrontalForcingsautoregressiontrendEnum;
    208208              else if (strcmp(name,"FrontalForcingsNumberofBasins")==0) return FrontalForcingsNumberofBasinsEnum;
    209209              else if (strcmp(name,"FrontalForcingsParam")==0) return FrontalForcingsParamEnum;
    210               else if (strcmp(name,"FrontalForcingsPhi")==0) return FrontalForcingsPhiEnum;
     210              else if (strcmp(name,"FrontalForcingsarlagcoefs")==0) return FrontalForcingsarlagcoefsEnum;
    211211              else if (strcmp(name,"GrdModel")==0) return GrdModelEnum;
    212212              else if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum;
    … …  
    489489              else if (strcmp(name,"SmbAutoregressiveOrder")==0) return SmbAutoregressiveOrderEnum;
    490490              else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum;
    491               else if (strcmp(name,"SmbBeta0")==0) return SmbBeta0Enum;
    492               else if (strcmp(name,"SmbBeta1")==0) return SmbBeta1Enum;
     491              else if (strcmp(name,"Smbautoregressionconst")==0) return SmbautoregressionconstEnum;
     492              else if (strcmp(name,"Smbautoregressiontrend")==0) return SmbautoregressiontrendEnum;
    493493              else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum;
    494494              else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum;
    … …  
    530530              else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
    531531              else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum;
    532               else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum;
     532              else if (strcmp(name,"Smbarlagcoefs")==0) return SmbarlagcoefsEnum;
    533533              else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum;
    534534              else if (strcmp(name,"SmbRefElevation")==0) return SmbRefElevationEnum;

  • issm/trunk-jpl/src/m/classes/SMBautoregression.m

    r26947 r27246  
    77        properties (SetAccess=public)
    88                num_basins        = 0;
    9                 beta0             = NaN;
    10                 beta1             = NaN;
     9                const             = NaN;
     10                trend             = NaN;
    1111                ar_order          = 0;
    1212                ar_initialtime    = 0;
    1313                ar_timestep       = 0;
    14                 phi               = NaN;
     14                arlag_coefs       = NaN;
    1515                basin_id          = NaN;
    1616                lapserates        = NaN;
    … …  
    3737                end % }}}
    3838                function self = initialize(self,md) % {{{
    39                         if isnan(self.beta1)
    40                                 self.beta1 = zeros(1,self.num_basins); %no trend in SMB
    41                                 disp('      smb.beta1 (trend) not specified: value set to 0');
     39                        if isnan(self.trend)
     40                                self.trend = zeros(1,self.num_basins); %no trend in SMB
     41                                disp('      smb.trend (trend) not specified: value set to 0');
    4242                        end
    4343                        if (self.ar_order==0)
    4444                                self.ar_order = 1; %dummy 1 value for autoregression
    45                                 self.phi      = zeros(self.num_basins,self.ar_order); %autoregression coefficients all set to 0
     45                                self.arlag_coefs      = zeros(self.num_basins,self.ar_order); %autoregression coefficients all set to 0
    4646                                disp('      smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0');
    4747                        end
    … …  
    5454                                disp('      smb.ar_timestep (timestep of autoregressive model) not specified: set to md.timestepping.time_step');
    5555                        end
    56                         if isnan(self.phi)
    57                                 self.phi = zeros(self.num_basins,self.ar_order); %autoregression model of order 0
    58                                 disp('      smb.phi (lag coefficients) not specified: order of autoregressive model set to 0');
     56                        if isnan(self.arlag_coefs)
     57                                self.arlag_coefs = zeros(self.num_basins,self.ar_order); %autoregression model of order 0
     58                                disp('      smb.arlag_coefs (lag coefficients) not specified: order of autoregressive model set to 0');
    5959                        end
    6060                end % }}}
    … …  
    6767                                md = checkfield(md,'fieldname','smb.num_basins','numel',1,'NaN',1,'Inf',1,'>',0);
    6868                                md = checkfield(md,'fieldname','smb.basin_id','Inf',1,'>=',0,'<=',md.smb.num_basins,'size',[md.mesh.numberofelements,1]);
    69                                 md = checkfield(md,'fieldname','smb.beta0','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); %scheme fails if passed as column vector
    70                                 md = checkfield(md,'fieldname','smb.beta1','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); %scheme fails if passed as column vector
     69                                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
     70                                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
    7171                                md = checkfield(md,'fieldname','smb.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0);
    7272                                md = checkfield(md,'fieldname','smb.ar_initialtime','numel',1,'NaN',1,'Inf',1);
    7373                                md = checkfield(md,'fieldname','smb.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep
    74                                 md = checkfield(md,'fieldname','smb.phi','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]);
     74                                md = checkfield(md,'fieldname','smb.arlag_coefs','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]);
    7575
    7676                                if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1)
    … …  
    100100                        fielddisplay(self,'num_basins','number of different basins [unitless]');
    101101                        fielddisplay(self,'basin_id','basin number assigned to each element [unitless]');
    102                         fielddisplay(self,'beta0','basin-specific intercept values [m ice eq./yr] (if beta_1==0 mean=beta_0/(1-sum(phi)))');
    103                         fielddisplay(self,'beta1','basin-specific trend values [m ice eq. yr^(-2)]');
     102                        fielddisplay(self,'const','basin-specific constant values [m ice eq./yr]');
     103                        fielddisplay(self,'trend','basin-specific trend values [m ice eq. yr^(-2)]');
    104104                        fielddisplay(self,'ar_order','order of the autoregressive model [unitless]');
    105105                        fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]');
    106106                        fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]');
    107                         fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]');
     107                        fielddisplay(self,'arlag_coefs','basin-specific vectors of lag coefficients [unitless]');
    108108                        fielddisplay(self,'lapserates','basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin [m ice eq yr^-1 m^-1] (default: no lapse rate)');
    109109                        fielddisplay(self,'elevationbins','basin-specific separations between elevation bins, 1 row per basin, 1 column per limit between bins [m] (default: no basin separation)');
    … …  
    153153                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_timestep','format','Double','scale',yts);
    154154                        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
    155                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta0','format','DoubleMat','name','md.smb.beta0','scale',1./yts,'yts',yts);
    156                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',yts);
    157                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',yts);
     155                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','const','format','DoubleMat','name','md.smb.const','scale',1./yts,'yts',yts);
     156                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','trend','format','DoubleMat','name','md.smb.trend','scale',1./(yts^2),'yts',yts);
     157                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','arlag_coefs','format','DoubleMat','name','md.smb.arlag_coefs','yts',yts);
    158158                        WriteData(fid,prefix,'data',templapserates,'format','DoubleMat','name','md.smb.lapserates','scale',1./yts,'yts',yts);
    159159                        WriteData(fid,prefix,'data',tempelevationbins,'format','DoubleMat','name','md.smb.elevationbins');

  • issm/trunk-jpl/src/m/classes/SMBautoregression.py

    r26949 r27246  
    1616    def __init__(self, *args):  # {{{
    1717        self.num_basins = 0
    18         self.beta0 = np.nan
    19         self.beta1 = np.nan
     18        self.const = np.nan
     19        self.trend = np.nan
    2020        self.ar_order = 0
    2121        self.ar_initialtime = 0
    2222        self.ar_timestep = 0
    23         self.phi = np.nan
     23        self.arlag_coefs = np.nan
    2424        self.basin_id = np.nan
    2525        self.lapserates = np.nan
    … …  
    4141        s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]'))
    4242        s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]'))
    43         s += '{}\n'.format(fielddisplay(self, 'beta0', 'basin-specific intercept values [m ice eq./yr] (if beta_1==0 mean=beta_0/(1-sum(phi)))'))
    44         s += '{}\n'.format(fielddisplay(self, 'beta1', 'basin-specific trend values [m ice eq. yr^(-2)]'))
     43        s += '{}\n'.format(fielddisplay(self, 'const', 'basin-specific constant values [m ice eq./yr]'))
     44        s += '{}\n'.format(fielddisplay(self, 'trend', 'basin-specific trend values [m ice eq. yr^(-2)]'))
    4545        s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]'))
    4646        s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]'))
    4747        s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]'))
    48         s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]'))
     48        s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of lag coefficients [unitless]'))
    4949        s += '{}\n'.format(fielddisplay(self, 'lapserates', 'basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin [m ice eq yr^-1 m^-1] (default: no lapse rate)'))
    5050        s += '{}\n'.format(fielddisplay(self, 'elevationbins', 'basin-specific SMB lapse rates applied in range of SMB<0 [m ice eq yr^-1 m^-1] (default: no lapse rate)'))
    … …  
    7272
    7373    def initialize(self, md):  # {{{
    74         if np.all(np.isnan(self.beta1)):
    75             self.beta1 = np.zeros((1, self.num_basins)) # No trend in SMB
    76             print('      smb.beta1 (trend) not specified: value set to 0')
     74        if np.all(np.isnan(self.trend)):
     75            self.trend = np.zeros((1, self.num_basins)) # No trend in SMB
     76            print('      smb.trend (trend) not specified: value set to 0')
    7777        if self.ar_order == 0:
    7878            self.ar_order = 1 # Dummy 1 value for autoregression
    79             self.phi = np.zeros((self.num_basins, self.ar_order)) # Autorgression coefficients all set to 0
     79            self.arlag_coefs = np.zeros((self.num_basins, self.ar_order)) # Autorgression coefficients all set to 0
    8080            print('      smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0')
    8181        if self.ar_initialtime == 0:
    … …  
    8585            self.ar_timestep = md.timestepping.time_step # Autoregression model has no prescribed time step
    8686            print('      smb.ar_timestep (timestep of autoregressive model) not specified: set to md.timestepping.time_step')
    87         if np.all(np.isnan(self.phi)):
    88             self.phi = np.zeros((self.num_basins, self.ar_order)) # Autoregression model of order 0
    89             print('      smb.phi (lag coefficients) not specified: order of autoregressive model set to 0')
     87        if np.all(np.isnan(self.arlag_coefs)):
     88            self.arlag_coefs = np.zeros((self.num_basins, self.ar_order)) # Autoregression model of order 0
     89            print('      smb.arlag_coefs (lag coefficients) not specified: order of autoregressive model set to 0')
    9090        return self
    9191    # }}}
    … …  
    9595            md = checkfield(md, 'fieldname', 'smb.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0)
    9696            md = checkfield(md, 'fieldname', 'smb.basin_id', 'Inf', 1, '>=', 0, '<=', md.smb.num_basins, 'size', [md.mesh.numberofelements])
    97             if len(np.shape(self.beta0)) == 1:
    98                 self.beta0 = np.array([self.beta0])
    99                 self.beta1 = np.array([self.beta1])
    100             md = checkfield(md, 'fieldname', 'smb.beta0', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) # Scheme fails if passed as column vector
    101             md = checkfield(md, 'fieldname', 'smb.beta1', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) # Scheme fails if passed as column vector; NOTE: As opposed to MATLAB implementation, pass list
     97            if len(np.shape(self.const)) == 1:
     98                self.const = np.array([self.const])
     99                self.trend = np.array([self.trend])
     100            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
     101            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
    102102            md = checkfield(md, 'fieldname', 'smb.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0)
    103103            md = checkfield(md, 'fieldname', 'smb.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1)
    104104            md = checkfield(md, 'fieldname', 'smb.ar_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep
    105             md = checkfield(md, 'fieldname', 'smb.phi', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order])
     105            md = checkfield(md, 'fieldname', 'smb.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order])
    106106           
    107107            if(np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1):
    … …  
    168168        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ar_timestep', 'format', 'Double', 'scale', yts)
    169169        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
    170         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'beta0', 'format', 'DoubleMat', 'name', 'md.smb.beta0', 'scale', 1 / yts, 'yts', yts)
    171         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.smb.beta1', 'scale', 1 / (yts ** 2), 'yts', yts)
    172         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.smb.phi', 'yts', yts)
     170        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'const', 'format', 'DoubleMat', 'name', 'md.smb.const', 'scale', 1 / yts, 'yts', yts)
     171        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'trend', 'format', 'DoubleMat', 'name', 'md.smb.trend', 'scale', 1 / (yts ** 2), 'yts', yts)
     172        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.smb.arlag_coefs', 'yts', yts)
    173173        WriteData(fid, prefix, 'data', templapserates, 'name', 'md.smb.lapserates', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts)
    174174        WriteData(fid, prefix, 'data', tempelevationbins, 'name', 'md.smb.elevationbins', 'format', 'DoubleMat')

  • issm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.m

    r26836 r27246  
    1010        properties (SetAccess=public)
    1111                num_basins                = 0;
    12                 beta0                     = NaN;
    13       beta1                     = NaN;
     12                const                     = NaN;
     13      trend                     = NaN;
    1414      ar_order                  = 0;
    1515      ar_initialtime            = 0;
    1616      ar_timestep               = 0;
    17       phi                       = NaN;
     17      arlag_coefs               = NaN;
    1818      basin_id                  = NaN;
    1919                groundedice_melting_rate  = NaN;
    … …  
    5959                                md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','NaN',1,'Inf',1,'<=',0,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins);
    6060            md = checkfield(md,'fieldname','basalforcings.basin_id','Inf',1,'>=',0,'<=',md.basalforcings.num_basins,'size',[md.mesh.numberofelements,1]);
    61             md = checkfield(md,'fieldname','basalforcings.beta0','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins);
    62                                 md = checkfield(md,'fieldname','basalforcings.beta1','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins);
     61            md = checkfield(md,'fieldname','basalforcings.const','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins);
     62                                md = checkfield(md,'fieldname','basalforcings.trend','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins);
    6363                                md = checkfield(md,'fieldname','basalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0);
    6464            md = checkfield(md,'fieldname','basalforcings.ar_initialtime','numel',1,'NaN',1,'Inf',1);
    6565            md = checkfield(md,'fieldname','basalforcings.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep
    66             md = checkfield(md,'fieldname','basalforcings.phi','NaN',1,'Inf',1,'size',[md.basalforcings.num_basins,md.basalforcings.ar_order]);
     66            md = checkfield(md,'fieldname','basalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.basalforcings.num_basins,md.basalforcings.ar_order]);
    6767                        end
    6868                        if ismember('BalancethicknessAnalysis',analyses),
    … …  
    8282                        fielddisplay(self,'num_basins','number of different basins [unitless]');
    8383         fielddisplay(self,'basin_id','basin number assigned to each element [unitless]');
    84          fielddisplay(self,'beta0','basin-specific intercept values [m/yr] (if beta_1==0 mean=beta_0/(1-sum(phi)))');
    85          fielddisplay(self,'beta1','basin-specific trend values [m  yr^(-2)]');
     84         fielddisplay(self,'const','basin-specific constant values [m/yr]');
     85         fielddisplay(self,'trend','basin-specific trend values [m  yr^(-2)]');
    8686         fielddisplay(self,'ar_order','order of the autoregressive model [unitless]');
    8787         fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]');
    8888         fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]');
    89          fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]');
     89         fielddisplay(self,'arlag_coefs','basin-specific vectors of lag coefficients [unitless]');
    9090                        fielddisplay(self,'deepwater_elevation','basin-specific elevation of ocean deepwater [m]');
    9191                        fielddisplay(self,'upperwater_melting_rate','basin-specific basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]');
    … …  
    107107         WriteData(fid,prefix,'object',self,'fieldname','ar_timestep','format','Double','scale',yts);
    108108         WriteData(fid,prefix,'object',self,'fieldname','basin_id','data',self.basin_id-1,'name','md.basalforcings.basin_id','format','IntMat','mattype',2); %0-indexed
    109          WriteData(fid,prefix,'object',self,'fieldname','beta0','format','DoubleMat','name','md.basalforcings.beta0','scale',1./yts,'yts',yts);
    110          WriteData(fid,prefix,'object',self,'fieldname','beta1','format','DoubleMat','name','md.basalforcings.beta1','scale',1./(yts^2),'yts',yts);
    111          WriteData(fid,prefix,'object',self,'fieldname','phi','format','DoubleMat','name','md.basalforcings.phi','yts',yts);   
     109         WriteData(fid,prefix,'object',self,'fieldname','const','format','DoubleMat','name','md.basalforcings.const','scale',1./yts,'yts',yts);
     110         WriteData(fid,prefix,'object',self,'fieldname','trend','format','DoubleMat','name','md.basalforcings.trend','scale',1./(yts^2),'yts',yts);
     111         WriteData(fid,prefix,'object',self,'fieldname','arlag_coefs','format','DoubleMat','name','md.basalforcings.arlag_coefs','yts',yts);   
    112112                        WriteData(fid,prefix,'object',self,'fieldname','deepwater_elevation','format','DoubleMat','name','md.basalforcings.deepwater_elevation');
    113113                        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/autoregressionlinearbasalforcings.py

    r26905 r27246  
    1515    def __init__(self, *args):  # {{{
    1616        self.num_basins = 0
    17         self.beta0 = np.nan
    18         self.beta1 = np.nan
     17        self.const = np.nan
     18        self.trend = np.nan
    1919        self.ar_order = 0
    2020        self.ar_initialtime = 0
    2121        self.ar_timestep = 0
    22         self.phi = np.nan
     22        self.arlag_coefs = np.nan
    2323        self.basin_id = np.nan
    2424        self.groundedice_melting_rate = np.nan
    … …  
    4040        s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]'))
    4141        s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]'))
    42         s += '{}\n'.format(fielddisplay(self, 'beta0', 'basin-specific intercept values [m ice eq./yr] (if beta_1==0 mean=beta_0/(1-sum(phi)))'))
    43         s += '{}\n'.format(fielddisplay(self, 'beta1', 'basin-specific trend values [m ice eq. yr^(-2)]'))
     42        s += '{}\n'.format(fielddisplay(self, 'const', 'basin-specific constant values [m ice eq./yr]'))
     43        s += '{}\n'.format(fielddisplay(self, 'trend', 'basin-specific trend values [m ice eq. yr^(-2)]'))
    4444        s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]'))
    4545        s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]'))
    4646        s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]'))
    47         s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]'))
     47        s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of lag coefficients [unitless]'))
    4848        s += '{}\n'.format(fielddisplay(self, 'deepwater_elevation', 'basin-specific elevation of ocean deepwater [m]'))
    4949        s += '{}\n'.format(fielddisplay(self, 'upperwater_melting_rate', 'basin-specic basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]'))
    … …  
    8484            md = checkfield(md, 'fieldname', 'basalforcings.basin_id', 'Inf', 1, '>=', 0, '<=', md.basalforcings.num_basins, 'size', [md.mesh.numberofelements])
    8585
    86             if len(np.shape(self.beta0)) == 1:
    87                 self.beta0 = np.array([self.beta0])
    88                 self.beta1 = np.array([self.beta1])
     86            if len(np.shape(self.const)) == 1:
     87                self.const = np.array([self.const])
     88                self.trend = np.array([self.trend])
    8989
    90             md = checkfield(md, 'fieldname', 'basalforcings.beta0', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) # Scheme fails if passed as column vector
    91             md = checkfield(md, 'fieldname', 'basalforcings.beta1', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) # Scheme fails if passed as column vector; NOTE: As opposed to MATLAB implementation, pass list
     90            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
     91            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
    9292            md = checkfield(md, 'fieldname', 'basalforcings.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0)
    9393            md = checkfield(md, 'fieldname', 'basalforcings.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1)
    9494            md = checkfield(md, 'fieldname', 'basalforcings.ar_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep
    95             md = checkfield(md, 'fieldname', 'basalforcings.phi', 'NaN', 1, 'Inf', 1, 'size', [md.basalforcings.num_basins, md.basalforcings.ar_order])
     95            md = checkfield(md, 'fieldname', 'basalforcings.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.basalforcings.num_basins, md.basalforcings.ar_order])
    9696        if 'BalancethicknessAnalysis' in analyses:
    9797            raise Exception('not implemented yet!')
    … …  
    113113        WriteData(fid, prefix, 'object', self, 'fieldname', 'ar_timestep', 'format', 'Double', 'scale', yts)
    114114        WriteData(fid, prefix, 'object', self, 'fieldname', 'basin_id', 'data', self.basin_id - 1, 'name', 'md.basalforcings.basin_id', 'format', 'IntMat', 'mattype', 2)  # 0-indexed
    115         WriteData(fid, prefix, 'object', self, 'fieldname', 'beta0', 'format', 'DoubleMat', 'name', 'md.basalforcings.beta0', 'scale', 1 / yts, 'yts', yts)
    116         WriteData(fid, prefix, 'object', self, 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.basalforcings.beta1', 'scale', 1 / (yts ** 2), 'yts', yts)
    117         WriteData(fid, prefix, 'object', self, 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.basalforcings.phi', 'yts', yts)
     115        WriteData(fid, prefix, 'object', self, 'fieldname', 'const', 'format', 'DoubleMat', 'name', 'md.basalforcings.const', 'scale', 1 / yts, 'yts', yts)
     116        WriteData(fid, prefix, 'object', self, 'fieldname', 'trend', 'format', 'DoubleMat', 'name', 'md.basalforcings.trend', 'scale', 1 / (yts ** 2), 'yts', yts)
     117        WriteData(fid, prefix, 'object', self, 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.basalforcings.arlag_coefs', 'yts', yts)
    118118        WriteData(fid, prefix, 'object', self, 'fieldname', 'deepwater_elevation', 'format', 'DoubleMat', 'name', 'md.basalforcings.deepwater_elevation')
    119119        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/src/m/classes/frontalforcingsrignotautoregression.m

    r26832 r27246  
    77        properties (SetAccess=public)
    88                num_basins           = 0;
    9                 beta0                = NaN;
    10                 beta1                = NaN;
     9                const                = NaN;
     10                trend                = NaN;
    1111                ar_order             = 0;
    1212                ar_initialtime       = 0;
    1313                ar_timestep          = 0;
    14                 phi                  = NaN;
     14                arlag_coef s         = NaN;
    1515                basin_id             = NaN;
    1616                subglacial_discharge = NaN;
    … …  
    5353         md = checkfield(md,'fieldname','frontalforcings.basin_id','Inf',1,'>=',0,'<=',md.frontalforcings.num_basins,'size',[md.mesh.numberofelements 1]);
    5454                        md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'NaN',1,'Inf',1,'timeseries',1);
    55                         md = checkfield(md,'fieldname','frontalforcings.beta0','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins);
    56          md = checkfield(md,'fieldname','frontalforcings.beta1','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins);
     55                        md = checkfield(md,'fieldname','frontalforcings.const','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins);
     56         md = checkfield(md,'fieldname','frontalforcings.trend','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins);
    5757         md = checkfield(md,'fieldname','frontalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0);
    5858         md = checkfield(md,'fieldname','frontalforcings.ar_initialtime','numel',1,'NaN',1,'Inf',1);
    5959         md = checkfield(md,'fieldname','frontalforcings.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep
    60                         md = checkfield(md,'fieldname','frontalforcings.phi','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]);
     60                        md = checkfield(md,'fieldname','frontalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]);
    6161
    6262                end % }}}
    … …  
    6666         fielddisplay(self,'basin_id','basin number assigned to each element [unitless]');
    6767         fielddisplay(self,'subglacial_discharge','sum of subglacial discharge for each basin [m/d]');
    68          fielddisplay(self,'beta0','basin-specific intercept values [∘C] (if beta_1==0 mean=beta_0/(1-sum(phi)))');
    69          fielddisplay(self,'beta1','basin-specific trend values [∘C yr^(-1)]');
     68         fielddisplay(self,'const','basin-specific constant term [∘C]');
     69         fielddisplay(self,'trend','basin-specific trend values [∘C yr^(-1)]');
    7070         fielddisplay(self,'ar_order','order of the autoregressive model [unitless]');
    7171         fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]');
    7272         fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]');
    73          fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]');
     73         fielddisplay(self,'arlag_coefs','basin-specific vectors of lag coefficients [unitless]');
    7474                end % }}}
    7575                function marshall(self,prefix,md,fid) % {{{
    … …  
    8282         WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_timestep','format','Double','scale',yts);
    8383         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
    84          WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','beta0','format','DoubleMat','name','md.frontalforcings.beta0');
    85          WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','beta1','format','DoubleMat','name','md.frontalforcings.beta1','scale',1./yts,'yts',yts);
    86          WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','phi','format','DoubleMat','name','md.frontalforcings.phi','yts',yts);
     84         WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','const','format','DoubleMat','name','md.frontalforcings.const');
     85         WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','trend','format','DoubleMat','name','md.frontalforcings.trend','scale',1./yts,'yts',yts);
     86         WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','arlag_coefs','format','DoubleMat','name','md.frontalforcings.arlag_coefs','yts',yts);
    8787                end % }}}
    8888        end

  • issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.py

    r26905 r27246  
    1616    def __init__(self, *args):  # {{{
    1717        self.num_basins = 0
    18         self.beta0 = np.nan
    19         self.beta1 = np.nan
     18        self.const = np.nan
     19        self.trend = np.nan
    2020        self.ar_order = 0
    2121        self.ar_initialtime = 0
    2222        self.ar_timestep = 0
    23         self.phi = np.nan
     23        self.arlag_coefs = np.nan
    2424        self.basin_id = np.nan
    2525        self.subglacial_discharge = np.nan
    … …  
    3535        s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]'))
    3636        s += '{}\n'.format(fielddisplay(self, 'subglacial_discharge', 'sum of subglacial discharge for each basin [m/d]'))
    37         s += '{}\n'.format(fielddisplay(self, 'beta0', 'basin-specific intercept values [°C] (if beta_1==0 mean=beta_0/(1-sum(phi)))'))
    38         s += '{}\n'.format(fielddisplay(self, 'beta1', 'basin-specific trend values [°C yr^(-1)]'))
     37        s += '{}\n'.format(fielddisplay(self, 'const', 'basin-specific constant values [°C]'))
     38        s += '{}\n'.format(fielddisplay(self, 'trend', 'basin-specific trend values [°C yr^(-1)]'))
    3939        s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]'))
    4040        s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]'))
    4141        s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]'))
    42         s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]'))
     42        s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of lag coefficients [unitless]'))
    4343        return s
    4444    #}}}
    … …  
    6060        md = checkfield(md, 'fieldname', 'frontalforcings.basin_id', 'Inf', 1, '>=', 0, '<=', md.frontalforcings.num_basins, 'size', [md.mesh.numberofelements])
    6161        md = checkfield(md, 'fieldname', 'frontalforcings.subglacial_discharge', '>=', 0, 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    62         if len(np.shape(self.beta0)) == 1:
    63             self.beta0 = np.array([self.beta0])
    64             self.beta1 = np.array([self.beta1])
    65         md = checkfield(md, 'fieldname', 'frontalforcings.beta0', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins)
    66         md = checkfield(md, 'fieldname', 'frontalforcings.beta1', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins)
     62        if len(np.shape(self.const)) == 1:
     63            self.const = np.array([self.const])
     64            self.trend = np.array([self.trend])
     65        md = checkfield(md, 'fieldname', 'frontalforcings.const', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins)
     66        md = checkfield(md, 'fieldname', 'frontalforcings.trend', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins)
    6767        md = checkfield(md, 'fieldname', 'frontalforcings.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0)
    6868        md = checkfield(md, 'fieldname', 'frontalforcings.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1)
    6969        md = checkfield(md, 'fieldname', 'frontalforcings.ar_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep
    70         md = checkfield(md, 'fieldname', 'frontalforcings.phi', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ar_order])
     70        md = checkfield(md, 'fieldname', 'frontalforcings.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ar_order])
    7171        return md
    7272    # }}}
    … …  
    8686        WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ar_timestep', 'format', 'Double', 'scale', yts)
    8787        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
    88         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'beta0', 'format', 'DoubleMat', 'name', 'md.frontalforcings.beta0')
    89         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.frontalforcings.beta1', 'scale', 1 / yts, 'yts', yts)
    90         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.frontalforcings.phi', 'yts', yts)
     88        WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'const', 'format', 'DoubleMat', 'name', 'md.frontalforcings.const')
     89        WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'trend', 'format', 'DoubleMat', 'name', 'md.frontalforcings.trend', 'scale', 1 / yts, 'yts', yts)
     90        WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.frontalforcings.arlag_coefs', 'yts', yts)
    9191    # }}}

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

    r26947 r27246  
    3939md.smb.num_basins     = 3; %number of basins
    4040md.smb.basin_id       = idbasin; %prescribe basin ID number to elements
    41 md.smb.beta0          = [0.5,1.2,1.5]; %intercept values of SMB in basins [m ice eq./yr]
    42 md.smb.beta1          = [0.0,0.01,-0.01]; %trend values of SMB in basins [m ice eq./yr^2]
     41md.smb.const          = [0.5,1.2,1.5]; %intercept values of SMB in basins [m ice eq./yr]
     42md.smb.trend          = [0.0,0.01,-0.01]; %trend values of SMB in basins [m ice eq./yr^2]
    4343md.smb.ar_initialtime = md.timestepping.start_time;
    4444md.smb.ar_order       = 4;
    4545md.smb.ar_timestep    = 2.0; %timestep of the autoregressive model [yr]
    46 md.smb.phi            = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1];[0.4,-0.4,0.1,-0.1]];
     46md.smb.arlag_coefs            = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1];[0.4,-0.4,0.1,-0.1]];
    4747md.smb.lapserates     = [0.01,0.0;0.01,-0.01;0.0,-0.01];
    4848md.smb.elevationbins  = [100;150;100];

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

    r26947 r27246  
    4444md.smb.num_basins = 3  # number of basins
    4545md.smb.basin_id = idbasin  # prescribe basin ID number to elements;
    46 md.smb.beta0 = np.array([[0.5, 1.2, 1.5]])  # intercept values of SMB in basins [m ice eq./yr]
    47 md.smb.beta1 = np.array([[0.0, 0.01, -0.01]])  # trend values of SMB in basins [m ice eq./yr^2]
     46md.smb.const = np.array([[0.5, 1.2, 1.5]])  # intercept values of SMB in basins [m ice eq./yr]
     47md.smb.trend = np.array([[0.0, 0.01, -0.01]])  # trend values of SMB in basins [m ice eq./yr^2]
    4848md.smb.ar_initialtime = md.timestepping.start_time
    4949md.smb.ar_order = 4
    5050md.smb.ar_timestep = 2.0  #timestep of the autoregressive model [yr]
    51 md.smb.phi = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1], [0.4, -0.4, 0.1, -0.1]])
     51md.smb.arlag_coefs = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1], [0.4, -0.4, 0.1, -0.1]])
    5252md.smb.lapserates        = np.array([[0.01,0.0],[0.01,-0.01],[0.0,-0.01]])
    5353md.smb.elevationbins  = np.array([100,150,100]).reshape(md.smb.num_basins,1)

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

    r26665 r27246  
    4949md.frontalforcings.basin_id             = idb_tf;
    5050md.frontalforcings.subglacial_discharge = 0.1*ones(md.mesh.numberofvertices,1);
    51 md.frontalforcings.beta0                = [0.05,0.01]; %intercept values of TF in basins [C]
    52 md.frontalforcings.beta1                = [0.001,0.0001]; %trend values of TF in basins [C/yr]
     51md.frontalforcings.const                = [0.05,0.01]; %intercept values of TF in basins [C]
     52md.frontalforcings.trend                = [0.001,0.0001]; %trend values of TF in basins [C/yr]
    5353md.frontalforcings.ar_initialtime       = md.timestepping.start_time; %initial time in the AR model parameterization [yr]
    5454md.frontalforcings.ar_order             = 4;
    5555md.frontalforcings.ar_timestep          = 2; %timestep of the autoregressive model [yr]
    56 md.frontalforcings.phi                  = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters
     56md.frontalforcings.arlag_coefs                  = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters
    5757
    5858% Floating Ice Melt parameters

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

    r26907 r27246  
    5454md.frontalforcings.basin_id = idb_tf
    5555md.frontalforcings.subglacial_discharge = 0.1 * np.ones((md.mesh.numberofvertices,))
    56 md.frontalforcings.beta0 = np.array([[0.05, 0.01]])  # intercept values of TF in basins [C]
    57 md.frontalforcings.beta1 = np.array([[0.001, 0.0001]])  # trend values of TF in basins [C/yr]
     56md.frontalforcings.const = np.array([[0.05, 0.01]])  # intercept values of TF in basins [C]
     57md.frontalforcings.trend = np.array([[0.001, 0.0001]])  # trend values of TF in basins [C/yr]
    5858md.frontalforcings.ar_initialtime = md.timestepping.start_time  # initial time in the AR model parameterization [yr]
    5959md.frontalforcings.ar_order = 4
    6060md.frontalforcings.ar_timestep = 2  # timestep of the autoregressive model [yr]
    61 md.frontalforcings.phi = np.array([[0.1, -0.1, 0.01, -0.01], [0.2, -0.2, 0.1, 0.0]])  # autoregressive parameters
     61md.frontalforcings.arlag_coefs = np.array([[0.1, -0.1, 0.01, -0.01], [0.2, -0.2, 0.1, 0.0]])  # autoregressive parameters
    6262#Floating Ice Melt parameters
    6363md.basalforcings.floatingice_melting_rate = 0.1 * np.ones((md.mesh.numberofvertices,))

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

    r26837 r27246  
    2727md.smb.num_basins     = nb_bas; %number of basins
    2828md.smb.basin_id       = idb; %prescribe basin ID number to elements
    29 md.smb.beta0          = [0.5,1.2]; %intercept values of SMB in basins [m ice eq./yr]
    30 md.smb.beta1          = [0.0,0.01]; %trend values of SMB in basins [m ice eq./yr^2]
     29md.smb.const          = [0.5,1.2]; %intercept values of SMB in basins [m ice eq./yr]
     30md.smb.trend          = [0.0,0.01]; %trend values of SMB in basins [m ice eq./yr^2]
    3131md.smb.ar_initialtime = md.timestepping.start_time;
    3232md.smb.ar_order       = 4;
    3333md.smb.ar_timestep    = 2.0; %timestep of the autoregressive model [yr]
    34 md.smb.phi            = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1]];
     34md.smb.arlag_coefs            = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1]];
    3535
    3636%Calving
    … …  
    4545md.basalforcings.num_basins     = nb_bas; %number of basins
    4646md.basalforcings.basin_id       = idb; %prescribe basin ID number to elements
    47 md.basalforcings.beta0          = [1.0,2.50]; %intercept values of DeepwaterMelt in basins [m/yr]
    48 md.basalforcings.beta1          = [0.2,0.01]; %trend values of DeepwaterMelt in basins [m/yr^2]
     47md.basalforcings.const          = [1.0,2.50]; %intercept values of DeepwaterMelt in basins [m/yr]
     48md.basalforcings.trend          = [0.2,0.01]; %trend values of DeepwaterMelt in basins [m/yr^2]
    4949md.basalforcings.ar_initialtime = md.timestepping.start_time;
    5050md.basalforcings.ar_order       = 1;
    5151md.basalforcings.ar_timestep    = 1.0; %timestep of the autoregressive model [yr]
    52 md.basalforcings.phi            = [0.0;0.1];
     52md.basalforcings.arlag_coefs            = [0.0;0.1];
    5353md.basalforcings.deepwater_elevation       = [-1000,-1520];
    5454md.basalforcings.upperwater_elevation      = [0,-50];

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

    r26907 r27246  
    3737md.smb.num_basins = nb_bas  # number of basins
    3838md.smb.basin_id = idb  # prescribe basin ID number to elements;
    39 md.smb.beta0 = np.array([[0.5, 1.2]])  # intercept values of SMB in basins [m ice eq./yr]
    40 md.smb.beta1 = np.array([[0.0, 0.01]])  # trend values of SMB in basins [m ice eq./yr^2]
     39md.smb.const = np.array([[0.5, 1.2]])  # intercept values of SMB in basins [m ice eq./yr]
     40md.smb.trend = np.array([[0.0, 0.01]])  # trend values of SMB in basins [m ice eq./yr^2]
    4141md.smb.ar_initialtime = md.timestepping.start_time
    4242md.smb.ar_order = 4
    4343md.smb.ar_timestep = 2.0  #timestep of the autoregressive model [yr]
    44 md.smb.phi = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1]])
     44md.smb.arlag_coefs = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1]])
    4545
    4646#Calving
    … …  
    5555md.basalforcings.num_basins = nb_bas
    5656md.basalforcings.basin_id  = idb
    57 md.basalforcings.beta0 = np.array([[1.0, 2.50]])  # intercept values of DeepwaterMelt in basins [m/yr]
    58 md.basalforcings.beta1  = np.array([[0.2, 0.01]])  # trend values of DeepwaterMelt in basins [m/yr^2]
     57md.basalforcings.const = np.array([[1.0, 2.50]])  # intercept values of DeepwaterMelt in basins [m/yr]
     58md.basalforcings.trend  = np.array([[0.2, 0.01]])  # trend values of DeepwaterMelt in basins [m/yr^2]
    5959md.basalforcings.ar_initialtime = md.timestepping.start_time  # initial time in the AR model parameterization [yr]
    6060md.basalforcings.ar_order  = 1
    6161md.basalforcings.ar_timestep  = 1.0  # timestep of the autoregressive model [yr]
    62 md.basalforcings.phi  = np.array([[0.0], [0.1]])  # autoregressive parameters
     62md.basalforcings.arlag_coefs  = np.array([[0.0], [0.1]])  # autoregressive parameters
    6363md.basalforcings.deepwater_elevation = np.array([[-1000, -1520]])
    6464md.basalforcings.upperwater_elevation = np.array([[0, -50]])
Note: See TracChangeset for help on using the changeset viewer.