Changeset 26495

Timestamp:

10/22/21 05:36:58 (3 years ago)

Author:

vverjans

Message:

CHG: added class frontalforcingsrignotautoregression (autoregressive calculation of thermal_forcing), autoregression routine in Element.cpp is now generic, corrected nightly runs 257 and 542, added nightly run 543

Location:

issm/trunk-jpl

Files:

• src/c/analyses/LevelsetAnalysis.cpp (modified) (2 diffs)
• src/c/classes/Elements/Element.cpp (modified) (2 diffs)
• src/c/classes/Elements/Element.h (modified) (2 diffs)
• src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (modified) (3 diffs)
• src/c/modules/FrontalForcingsx/FrontalForcingsx.h (modified) (1 diff)
• src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (6 diffs)
• src/c/shared/Enum/Enum.vim (modified) (12 diffs)
• src/c/shared/Enum/EnumDefinitions.h (modified) (3 diffs)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (3 diffs)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (25 diffs)
• src/c/shared/io/Marshalling/IoCodeConversions.cpp (modified) (1 diff)
• src/m/classes/frontalforcingsrignotautoregression.m (added)
• test/Archives/Archive257.arch (modified) ( previous)
• test/Archives/Archive542.arch (modified) ( previous)
• test/Archives/Archive543.arch (added)
• test/NightlyRun/test257.m (modified) (1 diff)
• test/NightlyRun/test542.m (modified) (1 diff)
• test/NightlyRun/test543.m (added)

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

@@ -93 +93 @@
                         break;
                 case FrontalForcingsRignotEnum:
-                        iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.basin_id",FrontalForcingsBasinIdEnum);
-                        iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.subglacial_discharge",FrontalForcingsSubglacialDischargeEnum);
-                        iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.thermalforcing",FrontalForcingsThermalForcingEnum);
-                        break;
+         /*Retrieve thermal forcing only in the case of non-autoregressive FrontalForcingsRignot*/
+         iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.thermalforcing",FrontalForcingsThermalForcingEnum);
+         /* Do not break here, still retrieve basin_ID,subglacial_discharge, etc.*/
+      case FrontalForcingsRignotAutoregressionEnum:
+         iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.basin_id",FrontalForcingsBasinIdEnum);
+         iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.subglacial_discharge",FrontalForcingsSubglacialDischargeEnum);
+         break; 
                 default:
                         _error_("Frontal forcings"<<EnumToStringx(melt_parameterization)<<" not supported yet");
@@ -137 +140 @@
         int melt_parameterization;
         iomodel->FindConstant(&melt_parameterization,"md.frontalforcings.parameterization");
+        int M,N;
+        IssmDouble* transparam = NULL;
         switch(melt_parameterization){
                 case FrontalForcingsDefaultEnum:
                         break;
+                case FrontalForcingsRignotAutoregressionEnum:
+                        /*Retrieve autoregressive parameters*/
+                        parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.randomflag",FrontalForcingsRandomflagEnum));
+         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_order",FrontalForcingsAutoregressiveOrderEnum));
+         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_initialtime",FrontalForcingsAutoregressionInitialTimeEnum));
+         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_timestep",FrontalForcingsAutoregressionTimestepEnum));
+                        iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.beta0");
+         parameters->AddObject(new DoubleVecParam(FrontalForcingsBeta0Enum,transparam,N));
+         xDelete<IssmDouble>(transparam);
+         iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.beta1");
+         parameters->AddObject(new DoubleVecParam(FrontalForcingsBeta1Enum,transparam,N));
+         xDelete<IssmDouble>(transparam);
+         iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.phi");
+         parameters->AddObject(new DoubleMatParam(FrontalForcingsPhiEnum,transparam,M,N));
+         xDelete<IssmDouble>(transparam);
+         iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.covmat");
+         parameters->AddObject(new DoubleMatParam(FrontalForcingsCovmatEnum,transparam,M,N));
+         xDelete<IssmDouble>(transparam);
+                        /*Do not break here, generic FrontalForcingsRignot parameters still to be retrieved*/
                 case FrontalForcingsRignotEnum:
                         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_basins",FrontalForcingsNumberofBasinsEnum));

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

@@ -59 +59 @@
         }
         return false;
+}/*}}}*/
+void       Element::AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noisespin,int enum_type){/*{{{*/
+
+        const int numvertices = this->GetNumberOfVertices();
+   int         basinid;
+   IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin;  
+   IssmDouble* phi_basin   = xNew<IssmDouble>(arorder);
+   IssmDouble* varspin     = xNewZeroInit<IssmDouble>(numvertices*arorder);
+
+   /*Get Basin ID and Basin coefficients*/
+        if(enum_type==SMBautoregressionEnum)                   this->GetInputValue(&basinid,SmbBasinsIdEnum);
+   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum);
+   for(int i=0;i<arorder;i++) phi_basin[i] = phi[basinid*arorder+i];
+   beta0_basin   = beta0[basinid];
+   beta1_basin   = beta1[basinid];
+
+   /*Loop over number of spin-up steps for all vertices*/   
+   for(int j=0;j<nspin;j++){  
+      tspin = starttime-((nspin-j)*tstep_ar); 
+      noisespin_basin = noisespin[j*numbasins+basinid];
+      IssmDouble* oldvarspin = xNew<IssmDouble>(numvertices*arorder);
+      for(int i=0;i<numvertices*arorder;i++) oldvarspin[i]=varspin[i];
+
+      for(int v=0;v<numvertices;v++){
+         IssmDouble autoregressionterm = 0.;
+         for(int i=0;i<arorder;i++) autoregressionterm += phi_basin[i]*varspin[v+i*numvertices];
+         varspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm+noisespin_basin;
+      }
+
+      /*Adjustt older values in varspin*/
+      for(int i=0;i<(arorder-1)*numvertices;i++) varspin[i+numvertices]=oldvarspin[i]; 
+
+      xDelete<IssmDouble>(oldvarspin); 
+   }
+        if(enum_type==SMBautoregressionEnum)                   this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder);
+   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder);
+
+ 
+   /*Cleanup and return*/
+   xDelete<IssmDouble>(varspin); 
+   xDelete<IssmDouble>(phi_basin);
+}/*}}}*/
+void       Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type){/*{{{*/
+   
+   const int numvertices = this->GetNumberOfVertices();
+   int         basinid,M,N;
+   IssmDouble  beta0_basin,beta1_basin,noise_basin;  
+   IssmDouble* phi_basin  = xNew<IssmDouble>(arorder);
+   IssmDouble* varlist     = xNew<IssmDouble>(numvertices);
+   IssmDouble* valuesautoregression = NULL;
+
+   /*Get Basin ID and Basin coefficients*/
+   if(enum_type==SMBautoregressionEnum)                   this->GetInputValue(&basinid,SmbBasinsIdEnum);
+   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum);
+   for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii]; 
+   beta0_basin   = beta0[basinid];
+   beta1_basin   = beta1[basinid];
+   noise_basin   = noiseterms[basinid];
+   if(enum_type==SMBautoregressionEnum)                   this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
+   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->GetArray(ThermalforcingValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
+
+   /*If not AR model timestep: take the old values of variable*/
+   if(isstepforar==false){
+      for(int i=0;i<numvertices;i++) varlist[i]=valuesautoregression[i];
+   }
+   /*If AR model timestep: compute variable values on vertices from AR*/
+   else{
+      for(int v=0;v<numvertices;v++){
+
+         /*Compute autoregressive term*/
+         IssmDouble autoregressionterm=0.;
+         for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*valuesautoregression[v+ii*numvertices];
+
+         /*Stochastic variable value*/
+         varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin;
+      }
+      /*Update autoregression TF values*/
+      IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
+      /*Assign newest values and shift older values*/
+      for(int i=0;i<numvertices;i++) temparray[i] = varlist[i]; 
+      for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = valuesautoregression[i];
+      if(enum_type==SMBautoregressionEnum)                   this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); 
+      if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); 
+      xDelete<IssmDouble>(temparray);
+   }
+   
+        //if(this->lid==55){_printf_("varlist: "<<varlist[0]<<"  "<<varlist[1]<<"  "<<varlist[2]<<"  "<<"noise_basin: "<<noise_basin<<"  "<<'\n');}
+
+        /*Add input to element*/
+   if(enum_type==SMBautoregressionEnum)                   this->AddInput(SmbMassBalanceEnum,varlist,P1Enum);
+   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum);
+
+   /*Cleanup*/
+   xDelete<IssmDouble>(phi_basin);
+   xDelete<IssmDouble>(varlist);
+   xDelete<IssmDouble>(valuesautoregression);
 }/*}}}*/
 void       Element::ComputeLambdaS(){/*{{{*/
@@ -3658 +3754 @@
 }
 /*}}}*/
-void       Element::SmbautoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noisespin){/*{{{*/
-
-        const int numvertices = this->GetNumberOfVertices();
-        int         basinid;
-        IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin;  
-        IssmDouble* phi_basin   = xNew<IssmDouble>(arorder);
-   IssmDouble* smbspin     = xNewZeroInit<IssmDouble>(numvertices*arorder);
-
-        /*Get Basin ID and Basin coefficients*/
-        this->GetInputValue(&basinid,SmbBasinsIdEnum);
-        for(int i=0;i<arorder;i++) phi_basin[i] = phi[basinid*arorder+i];
-        beta0_basin   = beta0[basinid];
-        beta1_basin   = beta1[basinid];
-
-        /*Loop over number of spin-up steps for all vertices*/  
-        for(int j=0;j<nspin;j++){       
-                tspin = starttime-((nspin-j)*tstep_ar); 
-      noisespin_basin = noisespin[j*numbasins+basinid];
-      IssmDouble* oldsmbspin = xNew<IssmDouble>(numvertices*arorder);
-      for(int i=0;i<numvertices*arorder;i++) oldsmbspin[i]=smbspin[i]; 
-
-      for(int v=0;v<numvertices;v++){
-         IssmDouble autoregressionterm = 0.;
-         for(int i=0;i<arorder;i++) autoregressionterm += phi_basin[i]*smbspin[v+i*numvertices];
-         smbspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm+noisespin_basin; //compute newest values in smbspin
-      }
-
-                /*Correct older values in smbspin*/
-      for(int i=0;i<(arorder-1)*numvertices;i++) smbspin[i+numvertices]=oldsmbspin[i]; 
-
-      xDelete<IssmDouble>(oldsmbspin); 
-        }
-        this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,smbspin,numvertices*arorder);
-
-        /*Cleanup and return*/
-   xDelete<IssmDouble>(smbspin); 
-   xDelete<IssmDouble>(phi_basin);
-}/*}}}*/
-void       Element::Smbautoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms){/*{{{*/
-        
-        const int numvertices = this->GetNumberOfVertices();
-        int         basinid,M,N;
-        IssmDouble  beta0_basin,beta1_basin,noise_basin;  
-        IssmDouble* phi_basin  = xNew<IssmDouble>(arorder);
-   IssmDouble* smblist    = xNew<IssmDouble>(numvertices);
-        IssmDouble* smbvaluesautoregression = NULL;
-
-        /*Get Basin ID and Basin coefficients*/
-        this->GetInputValue(&basinid,SmbBasinsIdEnum);
-        for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii]; 
-        beta0_basin   = beta0[basinid];
-        beta1_basin   = beta1[basinid];
-        noise_basin   = noiseterms[basinid]; 
-        this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&smbvaluesautoregression,&M);
-
-        /*If not AR model timestep: take the old SMB values*/
-        if(isstepforar==false){
-                /*VV do something with the lapse rate here if needed (12Oct2021)*/
-      for(int i=0;i<numvertices;i++) smblist[i]=smbvaluesautoregression[i]; 
-        }
-        /*If AR model timestep: compute SMB values on vertices from AR*/
-        else{
-                for(int v=0;v<numvertices;v++){
-
-                        /*Compute autoregressive term*/
-                        IssmDouble autoregressionterm=0.;
-                        for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];
-
-                        /*Stochastic SMB value*/
-                        smblist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin;
-                }
-                /*Update autoregression smb values*/
-                IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
-                /*Assign newest values and shift older values*/
-                for(int i=0;i<numvertices;i++) temparray[i] = smblist[i]; 
-                for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = smbvaluesautoregression[i];
-                this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); 
-                xDelete<IssmDouble>(temparray);
-        }
-        /*Add input to element*/
-        this->AddInput(SmbMassBalanceEnum,smblist,P1Enum);
-
-        /*Cleanup*/
-        xDelete<IssmDouble>(phi_basin);
-        xDelete<IssmDouble>(smblist);
-        xDelete<IssmDouble>(smbvaluesautoregression);
-}/*}}}*/
 void       Element::SmbGemb(IssmDouble timeinputs, int count){/*{{{*/
 

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

@@ -68 +68 @@
                 /*bool               AnyActive(void);*/
                 bool               AnyFSet(void);
+                void                                     AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noisespin,int enum_type);
+      void               Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type);
                 void               ComputeLambdaS(void);
                 void               ComputeNewDamage();
@@ -178 +180 @@
                 void               SmbSemic();
                 int                Sid();
-                void               SmbautoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noisespin);
-                void               Smbautoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms);
                 void               SmbGemb(IssmDouble timeinputs, int count);
                 void               StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);

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

@@ -6 +6 @@
 #include "../../shared/shared.h"
 #include "../../toolkits/toolkits.h"
+#include "../../shared/Random/random.h"
 
-void FrontalForcingsx(FemModel* femmodel){
+void FrontalForcingsx(FemModel* femmodel){/*{{{*/
 
         /*Recover melt_parameterization*/
@@ -17 +18 @@
                 case FrontalForcingsDefaultEnum:
                         break;
+                case FrontalForcingsRignotAutoregressionEnum:
+                        Thermalforcingautoregressionx(femmodel);
+                        /*Do not break here, call IcefrontAreax(),RignotMeltParameterizationx()*/
                 case FrontalForcingsRignotEnum:
                         femmodel->IcefrontAreax();
@@ -24 +28 @@
                         _error_("Frontal forcings "<<EnumToStringx(melt_parameterization)<<" not supported yet");
         }
-}
+}/*}}}*/
+void ThermalforcingautoregressionInitx(FemModel* femmodel){/*{{{*/
+
+   /*Initialization step of Thermalforcingautoregressionx*/
+   int M,N,Nphi,arorder,numbasins,my_rank;
+   IssmDouble starttime,tstep_ar,tinit_ar;
+   femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
+   femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
+   IssmDouble* beta0    = NULL;
+   IssmDouble* beta1    = NULL;
+   IssmDouble* phi      = NULL;
+   IssmDouble* covmat   = NULL;
+   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+   femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum);
+   femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
+   femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
+   femmodel->parameters->FindParam(&covmat,&M,&N,FrontalForcingsCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
+
+   /*AR model spin-up*/
+   int nspin{2*arorder+5};
+   IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin);
+   my_rank=IssmComm::GetRank();
+   if(my_rank==0){
+      bool randomflag{};
+      femmodel->parameters->FindParam(&randomflag,FrontalForcingsRandomflagEnum);
+      int fixedseed;
+      for(int i=0;i<nspin;i++){
+         IssmDouble* temparray = NULL;
+         /*Determine whether random seed is fixed to for loop step (randomflag==false) or random seed truly random (randomflag==true)*/
+         if(randomflag) fixedseed=-1;
+         else fixedseed = i;
+         multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed);
+         for(int j=0;j<numbasins;j++){noisespin[i*numbasins+j]=temparray[j];}
+         xDelete<IssmDouble>(temparray);
+      }
+   }
+   ISSM_MPI_Bcast(noisespin,numbasins*nspin,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+   /*Initialize ThermalforcingValuesAutoregressionEnum*/
+   for(Object* &object:femmodel->elements->objects){
+      Element* element      = xDynamicCast<Element*>(object); //generate element object
+      element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin,FrontalForcingsRignotAutoregressionEnum);
+   }
+
+   /*Cleanup*/
+   xDelete<IssmDouble>(beta0);
+   xDelete<IssmDouble>(beta1);
+   xDelete<IssmDouble>(phi);
+   xDelete<IssmDouble>(noisespin);
+   xDelete<IssmDouble>(covmat);
+}/*}}}*/
+void Thermalforcingautoregressionx(FemModel* femmodel){/*{{{*/
+
+   /*Get time parameters*/
+   IssmDouble time,dt,starttime,tstep_ar;
+   femmodel->parameters->FindParam(&time,TimeEnum);
+   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+   femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum);
+
+   /*Initialize module at first time step*/
+   if(time<=starttime+dt){ThermalforcingautoregressionInitx(femmodel);}
+   /*Determine if this is a time step for the AR model*/
+   bool isstepforar = false;
+
+   #ifndef _HAVE_AD_
+   if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt) isstepforar = true;
+   #else
+   _error_("not implemented yet");
+   #endif
+
+   /*Load parameters*/
+   int M,N,Nphi,arorder,numbasins,my_rank;
+   femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
+   femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
+   IssmDouble tinit_ar;
+   IssmDouble* beta0      = NULL;
+   IssmDouble* beta1      = NULL;
+   IssmDouble* phi        = NULL;
+   IssmDouble* covmat     = NULL;
+   IssmDouble* noiseterms = xNew<IssmDouble>(numbasins);
+   femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
+   femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum);    _assert_(M==numbasins);
+   femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
+   femmodel->parameters->FindParam(&covmat,&M,&N,FrontalForcingsCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
+
+   /*Time elapsed with respect to AR model initial time*/
+   IssmDouble telapsed_ar = time-tinit_ar;
+
+   /*Before looping through elements: compute noise term specific to each basin from covmat*/
+   my_rank=IssmComm::GetRank();
+   if(my_rank==0){
+      bool randomflag{};
+      femmodel->parameters->FindParam(&randomflag,FrontalForcingsRandomflagEnum);
+      /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
+      int fixedseed;
+      if(randomflag) fixedseed=-1;
+      else fixedseed = reCast<int,IssmDouble>((time-starttime)/dt);
+                /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
+      IssmDouble* temparray = NULL;
+      multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed);
+      for(int i=0;i<numbasins;i++) noiseterms[i]=temparray[i];
+      xDelete<IssmDouble>(temparray);
+   }
+   ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+   /*Loop over each element to compute SMB at vertices*/
+   for(Object* &object:femmodel->elements->objects){
+      Element* element = xDynamicCast<Element*>(object);
+      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,FrontalForcingsRignotAutoregressionEnum);
+   }
+
+   /*Cleanup*/
+   xDelete<IssmDouble>(beta0);
+   xDelete<IssmDouble>(beta1);
+   xDelete<IssmDouble>(phi);
+   xDelete<IssmDouble>(noiseterms);
+   xDelete<IssmDouble>(covmat);
+}/*}}}*/

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

@@ -7 +7 @@
 /* local prototypes: */
 void FrontalForcingsx(FemModel* femmodel);
+void ThermalforcingautoregressionInitx(FemModel* femmodel);
+void Thermalforcingautoregressionx(FemModel* femmodel);
 
 #endif

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

@@ -154 +154 @@
         femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
         femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
-        IssmDouble* beta0    = xNew<IssmDouble>(numbasins);
-        IssmDouble* beta1    = xNew<IssmDouble>(numbasins);
-        IssmDouble* phi      = xNew<IssmDouble>(numbasins*arorder);
-        IssmDouble* covmat   = xNew<IssmDouble>(numbasins*numbasins);
+        IssmDouble* beta0    = NULL; 
+        IssmDouble* beta1    = NULL;
+        IssmDouble* phi      = NULL;
+        IssmDouble* covmat   = NULL;
         femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
@@ -175 +175 @@
                 int fixedseed;
                 for(int i=0;i<nspin;i++){
-                        IssmDouble* temparray = xNew<IssmDouble>(numbasins);
+                        //IssmDouble* temparray = xNew<IssmDouble>(numbasins);
+                        IssmDouble* temparray = NULL;
                         /*Determine whether random seed is fixed to for loop step (randomflag==false) or random seed truly random (randomflag==true)*/
                         if(randomflag) fixedseed=-1;
@@ -189 +190 @@
         for(Object* &object:femmodel->elements->objects){
       Element* element      = xDynamicCast<Element*>(object); //generate element object
-                element->SmbautoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin);
+                element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin,SMBautoregressionEnum);
         }
         
@@ -224 +225 @@
         femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
         IssmDouble tinit_ar;
-        IssmDouble* beta0      = xNew<IssmDouble>(numbasins);
-        IssmDouble* beta1      = xNew<IssmDouble>(numbasins);
-        IssmDouble* phi        = xNew<IssmDouble>(numbasins*arorder);
+        IssmDouble* beta0      = NULL; 
+        IssmDouble* beta1      = NULL;
+        IssmDouble* phi        = NULL;
+        IssmDouble* covmat     = NULL;
         IssmDouble* noiseterms = xNew<IssmDouble>(numbasins);
-        IssmDouble* covmat     = xNew<IssmDouble>(numbasins*numbasins);
         femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
         femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
@@ -247 +248 @@
                 if(randomflag) fixedseed=-1;
                 else fixedseed =  reCast<int,IssmDouble>((time-starttime)/dt);
-                multivariateNormal(&noiseterms,numbasins,0.0,covmat,fixedseed);
+                /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/    
+                IssmDouble* temparray = NULL;
+                multivariateNormal(&temparray,numbasins,0.0,covmat,fixedseed);
+                for(int i=0;i<numbasins;i++) noiseterms[i]=temparray[i];
+                xDelete<IssmDouble>(temparray);
         }
         ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
@@ -254 +259 @@
         for(Object* &object:femmodel->elements->objects){
                 Element* element = xDynamicCast<Element*>(object);
-                element->Smbautoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms);
+                element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,SMBautoregressionEnum);
         }
 

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

@@ -177 +177 @@
 syn keyword cConstant FrictionVoidRatioEnum
 syn keyword cConstant FrontalForcingsBasinIcefrontAreaEnum
+syn keyword cConstant FrontalForcingsAutoregressionInitialTimeEnum
+syn keyword cConstant FrontalForcingsAutoregressionTimestepEnum
+syn keyword cConstant FrontalForcingsAutoregressiveOrderEnum
+syn keyword cConstant FrontalForcingsBeta0Enum
+syn keyword cConstant FrontalForcingsBeta1Enum
+syn keyword cConstant FrontalForcingsCovmatEnum
 syn keyword cConstant FrontalForcingsNumberofBasinsEnum
 syn keyword cConstant FrontalForcingsParamEnum
+syn keyword cConstant FrontalForcingsPhiEnum
+syn keyword cConstant FrontalForcingsRandomflagEnum
 syn keyword cConstant GrdModelEnum
 syn keyword cConstant GroundinglineFrictionInterpolationEnum
@@ -1016 +1024 @@
 syn keyword cConstant TemperaturePicardEnum
 syn keyword cConstant TemperatureSEMICEnum
+syn keyword cConstant ThermalforcingValuesAutoregressionEnum
 syn keyword cConstant ThermalSpctemperatureEnum
 syn keyword cConstant ThicknessAbsGradientEnum
@@ -1262 +1271 @@
 syn keyword cConstant FrontalForcingsDefaultEnum
 syn keyword cConstant FrontalForcingsRignotEnum
+syn keyword cConstant FrontalForcingsRignotAutoregressionEnum
 syn keyword cConstant FsetEnum
 syn keyword cConstant FullMeltOnPartiallyFloatingEnum
@@ -1553 +1563 @@
 syn keyword cType Cfsurfacesquare
 syn keyword cType Channel
+syn keyword cType classes
 syn keyword cType Constraint
 syn keyword cType Constraints
@@ -1559 +1570 @@
 syn keyword cType ControlInput
 syn keyword cType Covertree
+syn keyword cType DatasetInput
 syn keyword cType DataSetParam
-syn keyword cType DatasetInput
 syn keyword cType Definition
 syn keyword cType DependentObject
@@ -1573 +1584 @@
 syn keyword cType ElementInput
 syn keyword cType ElementMatrix
+syn keyword cType Elements
 syn keyword cType ElementVector
-syn keyword cType Elements
 syn keyword cType ExponentialVariogram
 syn keyword cType ExternalResult
@@ -1581 +1592 @@
 syn keyword cType Friction
 syn keyword cType Gauss
+syn keyword cType GaussianVariogram
+syn keyword cType gaussobjects
 syn keyword cType GaussPenta
 syn keyword cType GaussSeg
 syn keyword cType GaussTetra
 syn keyword cType GaussTria
-syn keyword cType GaussianVariogram
 syn keyword cType GenericExternalResult
 syn keyword cType GenericOption
@@ -1601 +1613 @@
 syn keyword cType IssmDirectApplicInterface
 syn keyword cType IssmParallelDirectApplicInterface
+syn keyword cType krigingobjects
 syn keyword cType Load
 syn keyword cType Loads
@@ -1611 +1624 @@
 syn keyword cType Matice
 syn keyword cType Matlitho
+syn keyword cType matrixobjects
 syn keyword cType MatrixParam
 syn keyword cType Misfit
@@ -1623 +1637 @@
 syn keyword cType Observations
 syn keyword cType Option
+syn keyword cType Options
 syn keyword cType OptionUtilities
-syn keyword cType Options
 syn keyword cType Param
 syn keyword cType Parameters
@@ -1638 +1652 @@
 syn keyword cType Regionaloutput
 syn keyword cType Results
+syn keyword cType Riftfront
 syn keyword cType RiftStruct
-syn keyword cType Riftfront
 syn keyword cType SealevelGeometry
 syn keyword cType Seg
 syn keyword cType SegInput
+syn keyword cType Segment
 syn keyword cType SegRef
-syn keyword cType Segment
 syn keyword cType SpcDynamic
 syn keyword cType SpcStatic
@@ -1663 +1677 @@
 syn keyword cType Vertex
 syn keyword cType Vertices
-syn keyword cType classes
-syn keyword cType gaussobjects
-syn keyword cType krigingobjects
-syn keyword cType matrixobjects
 syn keyword cType AdjointBalancethickness2Analysis
 syn keyword cType AdjointBalancethicknessAnalysis

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

@@ -171 +171 @@
         FrictionVoidRatioEnum,
         FrontalForcingsBasinIcefrontAreaEnum,
+        FrontalForcingsAutoregressionInitialTimeEnum,
+   FrontalForcingsAutoregressionTimestepEnum,
+   FrontalForcingsAutoregressiveOrderEnum,
+        FrontalForcingsBeta0Enum,
+   FrontalForcingsBeta1Enum,
+   FrontalForcingsCovmatEnum,
         FrontalForcingsNumberofBasinsEnum,
         FrontalForcingsParamEnum,
+   FrontalForcingsPhiEnum,
+        FrontalForcingsRandomflagEnum,
         GrdModelEnum,
         GroundinglineFrictionInterpolationEnum,
@@ -1013 +1021 @@
         TemperaturePicardEnum,
         TemperatureSEMICEnum,
+        ThermalforcingValuesAutoregressionEnum, 
         ThermalSpctemperatureEnum,
         ThicknessAbsGradientEnum,
@@ -1261 +1270 @@
         FrontalForcingsDefaultEnum,
         FrontalForcingsRignotEnum,
+        FrontalForcingsRignotAutoregressionEnum,
         FsetEnum,
         FullMeltOnPartiallyFloatingEnum,

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

@@ -179 +179 @@
                 case FrictionVoidRatioEnum : return "FrictionVoidRatio";
                 case FrontalForcingsBasinIcefrontAreaEnum : return "FrontalForcingsBasinIcefrontArea";
+                case FrontalForcingsAutoregressionInitialTimeEnum : return "FrontalForcingsAutoregressionInitialTime";
+                case FrontalForcingsAutoregressionTimestepEnum : return "FrontalForcingsAutoregressionTimestep";
+                case FrontalForcingsAutoregressiveOrderEnum : return "FrontalForcingsAutoregressiveOrder";
+                case FrontalForcingsBeta0Enum : return "FrontalForcingsBeta0";
+                case FrontalForcingsBeta1Enum : return "FrontalForcingsBeta1";
+                case FrontalForcingsCovmatEnum : return "FrontalForcingsCovmat";
                 case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins";
                 case FrontalForcingsParamEnum : return "FrontalForcingsParam";
+                case FrontalForcingsPhiEnum : return "FrontalForcingsPhi";
+                case FrontalForcingsRandomflagEnum : return "FrontalForcingsRandomflag";
                 case GrdModelEnum : return "GrdModel";
                 case GroundinglineFrictionInterpolationEnum : return "GroundinglineFrictionInterpolation";
@@ -1018 +1026 @@
                 case TemperaturePicardEnum : return "TemperaturePicard";
                 case TemperatureSEMICEnum : return "TemperatureSEMIC";
+                case ThermalforcingValuesAutoregressionEnum : return "ThermalforcingValuesAutoregression";
                 case ThermalSpctemperatureEnum : return "ThermalSpctemperature";
                 case ThicknessAbsGradientEnum : return "ThicknessAbsGradient";
@@ -1264 +1273 @@
                 case FrontalForcingsDefaultEnum : return "FrontalForcingsDefault";
                 case FrontalForcingsRignotEnum : return "FrontalForcingsRignot";
+                case FrontalForcingsRignotAutoregressionEnum : return "FrontalForcingsRignotAutoregression";
                 case FsetEnum : return "Fset";
                 case FullMeltOnPartiallyFloatingEnum : return "FullMeltOnPartiallyFloating";

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

@@ -182 +182 @@
               else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum;
               else if (strcmp(name,"FrontalForcingsBasinIcefrontArea")==0) return FrontalForcingsBasinIcefrontAreaEnum;
+              else if (strcmp(name,"FrontalForcingsAutoregressionInitialTime")==0) return FrontalForcingsAutoregressionInitialTimeEnum;
+              else if (strcmp(name,"FrontalForcingsAutoregressionTimestep")==0) return FrontalForcingsAutoregressionTimestepEnum;
+              else if (strcmp(name,"FrontalForcingsAutoregressiveOrder")==0) return FrontalForcingsAutoregressiveOrderEnum;
+              else if (strcmp(name,"FrontalForcingsBeta0")==0) return FrontalForcingsBeta0Enum;
+              else if (strcmp(name,"FrontalForcingsBeta1")==0) return FrontalForcingsBeta1Enum;
+              else if (strcmp(name,"FrontalForcingsCovmat")==0) return FrontalForcingsCovmatEnum;
               else if (strcmp(name,"FrontalForcingsNumberofBasins")==0) return FrontalForcingsNumberofBasinsEnum;
               else if (strcmp(name,"FrontalForcingsParam")==0) return FrontalForcingsParamEnum;
+              else if (strcmp(name,"FrontalForcingsPhi")==0) return FrontalForcingsPhiEnum;
+              else if (strcmp(name,"FrontalForcingsRandomflag")==0) return FrontalForcingsRandomflagEnum;
               else if (strcmp(name,"GrdModel")==0) return GrdModelEnum;
               else if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum;
@@ -252 +260 @@
               else if (strcmp(name,"InversionNsteps")==0) return InversionNstepsEnum;
               else if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum;
-              else if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
+         else stage=3;
+   }
+   if(stage==3){
+              if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
               else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
               else if (strcmp(name,"InversionType")==0) return InversionTypeEnum;
@@ -260 +271 @@
               else if (strcmp(name,"LevelsetReinitFrequency")==0) return LevelsetReinitFrequencyEnum;
               else if (strcmp(name,"LevelsetStabilization")==0) return LevelsetStabilizationEnum;
-         else stage=3;
-   }
-   if(stage==3){
-              if (strcmp(name,"LockFileName")==0) return LockFileNameEnum;
+              else if (strcmp(name,"LockFileName")==0) return LockFileNameEnum;
               else if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum;
               else if (strcmp(name,"LoveCoreMantleBoundary")==0) return LoveCoreMantleBoundaryEnum;
@@ -375 +383 @@
               else if (strcmp(name,"SolidearthSettingsViscous")==0) return SolidearthSettingsViscousEnum;
               else if (strcmp(name,"SealevelchangeGeometryDone")==0) return SealevelchangeGeometryDoneEnum;
-              else if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum;
+         else stage=4;
+   }
+   if(stage==4){
+              if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum;
               else if (strcmp(name,"SealevelchangeViscousTimes")==0) return SealevelchangeViscousTimesEnum;
               else if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum;
@@ -383 +394 @@
               else if (strcmp(name,"TidalLoveL")==0) return TidalLoveLEnum;
               else if (strcmp(name,"TidalLoveK2Secular")==0) return TidalLoveK2SecularEnum;
-         else stage=4;
-   }
-   if(stage==4){
-              if (strcmp(name,"LoadLoveH")==0) return LoadLoveHEnum;
+              else if (strcmp(name,"LoadLoveH")==0) return LoadLoveHEnum;
               else if (strcmp(name,"LoadLoveK")==0) return LoadLoveKEnum;
               else if (strcmp(name,"LoadLoveL")==0) return LoadLoveLEnum;
@@ -498 +506 @@
               else if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum;
               else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum;
-              else if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum;
+         else stage=5;
+   }
+   if(stage==5){
+              if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum;
               else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum;
               else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum;
@@ -506 +517 @@
               else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
               else if (strcmp(name,"ThermalNumRequestedOutputs")==0) return ThermalNumRequestedOutputsEnum;
-         else stage=5;
-   }
-   if(stage==5){
-              if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
+              else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
               else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
               else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
@@ -621 +629 @@
               else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
               else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
-              else if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum;
+         else stage=6;
+   }
+   if(stage==6){
+              if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum;
               else if (strcmp(name,"CalvingMeltingFluxLevelset")==0) return CalvingMeltingFluxLevelsetEnum;
               else if (strcmp(name,"Converged")==0) return ConvergedEnum;
@@ -629 +640 @@
               else if (strcmp(name,"DamageDbar")==0) return DamageDbarEnum;
               else if (strcmp(name,"DamageDbarOld")==0) return DamageDbarOldEnum;
-         else stage=6;
-   }
-   if(stage==6){
-              if (strcmp(name,"DamageF")==0) return DamageFEnum;
+              else if (strcmp(name,"DamageF")==0) return DamageFEnum;
               else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
               else if (strcmp(name,"DepthBelowSurface")==0) return DepthBelowSurfaceEnum;
@@ -744 +752 @@
               else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
               else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
-              else if (strcmp(name,"Ice")==0) return IceEnum;
+         else stage=7;
+   }
+   if(stage==7){
+              if (strcmp(name,"Ice")==0) return IceEnum;
               else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
               else if (strcmp(name,"Input")==0) return InputEnum;
@@ -752 +763 @@
               else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
               else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
-         else stage=7;
-   }
-   if(stage==7){
-              if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
+              else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
               else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
               else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
@@ -867 +875 @@
               else if (strcmp(name,"SealevelchangeGEsubelOcean")==0) return SealevelchangeGEsubelOceanEnum;
               else if (strcmp(name,"SealevelchangeGNsubelOcean")==0) return SealevelchangeGNsubelOceanEnum;
-              else if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum;
+         else stage=8;
+   }
+   if(stage==8){
+              if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum;
               else if (strcmp(name,"SealevelchangeGUsubelIce")==0) return SealevelchangeGUsubelIceEnum;
               else if (strcmp(name,"SealevelchangeGEsubelIce")==0) return SealevelchangeGEsubelIceEnum;
@@ -875 +886 @@
               else if (strcmp(name,"SealevelchangeGEsubelHydro")==0) return SealevelchangeGEsubelHydroEnum;
               else if (strcmp(name,"SealevelchangeGNsubelHydro")==0) return SealevelchangeGNsubelHydroEnum;
-         else stage=8;
-   }
-   if(stage==8){
-              if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum;
+              else if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum;
               else if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum;
               else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum;
@@ -990 +998 @@
               else if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
               else if (strcmp(name,"SmbValuesAutoregression")==0) return SmbValuesAutoregressionEnum;
-              else if (strcmp(name,"SmbV")==0) return SmbVEnum;
+         else stage=9;
+   }
+   if(stage==9){
+              if (strcmp(name,"SmbV")==0) return SmbVEnum;
               else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum;
               else if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
@@ -998 +1009 @@
               else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum;
               else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum;
-         else stage=9;
-   }
-   if(stage==9){
-              if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum;
+              else if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum;
               else if (strcmp(name,"SmbZY")==0) return SmbZYEnum;
               else if (strcmp(name,"SolidearthExternalDisplacementEastRate")==0) return SolidearthExternalDisplacementEastRateEnum;
@@ -1042 +1050 @@
               else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
               else if (strcmp(name,"TemperatureSEMIC")==0) return TemperatureSEMICEnum;
+              else if (strcmp(name,"ThermalforcingValuesAutoregression")==0) return ThermalforcingValuesAutoregressionEnum;
               else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum;
               else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
@@ -1112 +1121 @@
               else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
               else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
-              else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum;
+         else stage=10;
+   }
+   if(stage==10){
+              if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum;
               else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum;
               else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum;
@@ -1121 +1133 @@
               else if (strcmp(name,"Outputdefinition43")==0) return Outputdefinition43Enum;
               else if (strcmp(name,"Outputdefinition44")==0) return Outputdefinition44Enum;
-         else stage=10;
-   }
-   if(stage==10){
-              if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum;
+              else if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum;
               else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum;
               else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum;
@@ -1235 +1244 @@
               else if (strcmp(name,"Contact")==0) return ContactEnum;
               else if (strcmp(name,"Contour")==0) return ContourEnum;
-              else if (strcmp(name,"Contours")==0) return ContoursEnum;
+         else stage=11;
+   }
+   if(stage==11){
+              if (strcmp(name,"Contours")==0) return ContoursEnum;
               else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
               else if (strcmp(name,"ControlInputGrad")==0) return ControlInputGradEnum;
@@ -1244 +1256 @@
               else if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
               else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum;
-         else stage=11;
-   }
-   if(stage==11){
-              if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
+              else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
               else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
               else if (strcmp(name,"DataSet")==0) return DataSetEnum;
@@ -1294 +1303 @@
               else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum;
               else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum;
+              else if (strcmp(name,"FrontalForcingsRignotAutoregression")==0) return FrontalForcingsRignotAutoregressionEnum;
               else if (strcmp(name,"Fset")==0) return FsetEnum;
               else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum;
@@ -1357 +1367 @@
               else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
               else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
-              else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
+         else stage=12;
+   }
+   if(stage==12){
+              if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
               else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
               else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum;
@@ -1367 +1380 @@
               else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum;
               else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
-         else stage=12;
-   }
-   if(stage==12){
-              if (strcmp(name,"LoveKi")==0) return LoveKiEnum;
+              else if (strcmp(name,"LoveKi")==0) return LoveKiEnum;
               else if (strcmp(name,"LoveKr")==0) return LoveKrEnum;
               else if (strcmp(name,"LoveLi")==0) return LoveLiEnum;
@@ -1480 +1490 @@
               else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum;
               else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
-              else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
+         else stage=13;
+   }
+   if(stage==13){
+              if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
               else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
               else if (strcmp(name,"Scaled")==0) return ScaledEnum;
@@ -1490 +1503 @@
               else if (strcmp(name,"SealevelchangePolarMotion")==0) return SealevelchangePolarMotionEnum;
               else if (strcmp(name,"SealevelNmotion")==0) return SealevelNmotionEnum;
-         else stage=13;
-   }
-   if(stage==13){
-              if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum;
+              else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum;
               else if (strcmp(name,"SealevelchangeAnalysis")==0) return SealevelchangeAnalysisEnum;
               else if (strcmp(name,"Seg")==0) return SegEnum;

issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp

@@ -274 +274 @@
                 case 1: return FrontalForcingsDefaultEnum;
                 case 2: return FrontalForcingsRignotEnum;
+                case 55: return FrontalForcingsRignotAutoregressionEnum;
                 default: _error_("Marshalled Frontalforcings code \""<<enum_in<<"\" not supported yet");
         }

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

@@ -31 +31 @@
     end
 end
+
+%SMB parameters
+md.timestepping.start_time = 0;
+md.timestepping.time_step  = 1;
+md.timestepping.final_time = 5;
+md.smb                = SMBautoregression();
+md.smb.num_basins     = 3; %number of basins
+md.smb.basin_id       = idbasin; %prescribe basin ID number to elements
+md.smb.beta0          = [0.5,1.2,1.5]; %intercept values of SMB in basins [m ice eq./yr]
+md.smb.beta1          = [0.0,0.01,-0.01]; %trend values of SMB in basins [m ice eq./yr^2]
+md.smb.ar_initialtime = md.timestepping.start_time;
+md.smb.ar_order       = 4;
+md.smb.ar_timestep    = 2.0; %timestep of the autoregressive model [yr]
+md.smb.phi            = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1];[0.4,-0.4,0.1,-0.1]];
+md.smb.covmat         = [[0.15 0.08 -0.02];[0.08 0.12 -0.05];[-0.02 -0.05 0.1]];
+md.smb.randomflag     = 0; %fixed random seeds
 
 md=solve(md,'Transient');

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

@@ -9 +9 @@
 %separate domain in 2 basins
 idbasin = zeros(md.mesh.numberofelements,1);
-iid1    = find(md.mesh.y>=mean(md.mesh.y));
+iid1    = find(md.mesh.x<=-1.6e6);
 for ii=1:md.mesh.numberofelements
     for vertex=1:3