Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 26476)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 26477)
@@ -187,4 +187,5 @@
 	./shared/Elements/DrainageFunctionWaterfraction.cpp \
 	./shared/Elements/EstarComponents.cpp \
+	./shared/Random/random.cpp \
 	./shared/String/DescriptorIndex.cpp \
 	./toolkits/issm/IssmToolkitUtils.cpp \
Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp	(revision 26476)
+++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp	(revision 26477)
@@ -179,4 +179,7 @@
 			iomodel->FetchDataToInput(inputs,elements,"md.smb.b_min",SmbBMinEnum);
 			break;
+		case SMBautoregressionEnum:
+         iomodel->FetchDataToInput(inputs,elements,"md.smb.basin_id",SmbBasinsIdEnum);
+         break;
 		case SMBhenningEnum:
 			iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum,0.);
@@ -375,4 +378,7 @@
 			/*Nothing to add to parameters*/
 			break;
+		case SMBautoregressionEnum:
+         /*Nothing to add to parameters*/
+         break;
 		case SMBhenningEnum:
 			/*Nothing to add to parameters*/
@@ -467,4 +473,8 @@
 			SmbGradientsElax(femmodel);
 			break;
+		case SMBautoregressionEnum:
+         if(VerboseSolution())_printf0_("    call smb autoregression module\n");
+         Smbautoregressionx(femmodel);
+         break;
 		case SMBhenningEnum:
 			if(VerboseSolution())_printf0_("  call smb Henning module\n");
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26476)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26477)
@@ -3574,4 +3574,70 @@
 }
 /*}}}*/
+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; //initialize scalars 
+	IssmDouble* phi_basin               = xNew<IssmDouble>(arorder);
+   IssmDouble* smbspin                 = xNew<IssmDouble>(numvertices*arorder);
+	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];
+	for(int jj{0};jj<nspin;jj++){	
+		tspin = starttime-((nspin-jj)*tstep_ar); //current time in spin-up loop
+      noisespin_basin = noisespin[jj*numbasins+basinid];
+      IssmDouble* oldsmbspin = xNewZeroInit<IssmDouble>(numvertices*arorder);
+      for(int ii{0};ii<numvertices*arorder;ii++){oldsmbspin[ii]=smbspin[ii];} //copy smbspin in oldsmbspin
+      for(int v=0;v<numvertices;v++){
+         double autoregressionterm{0.0};
+         for(int ii{0};ii<arorder;ii++){autoregressionterm += phi_basin[ii]*smbspin[v+ii*numvertices];}
+         smbspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm+noisespin_basin; //compute newest values in smbspin
+      }
+      for(int ii{0};ii<(arorder-1)*numvertices;ii++){smbspin[ii+numvertices]=oldsmbspin[ii];} //correct older values in smbspin
+      xDelete<IssmDouble>(oldsmbspin); //cleanup
+	}
+	this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,smbspin,numvertices*arorder);//save spin-up autoregression values
+   xDelete<IssmDouble>(smbspin); //cleanup
+   xDelete<IssmDouble>(phi_basin); //cleanup
+}/*}}}*/
+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; //initialize scalars 
+	IssmDouble* phi_basin               = xNew<IssmDouble>(arorder);
+   IssmDouble* smblist                 = xNew<IssmDouble>(numvertices);
+	IssmDouble* smbvaluesautoregression = NULL; //array for past SMB values that we are about to retrieve
+	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]; //note that noiseterms is computed at every timestep
+	this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&smbvaluesautoregression,&M); //get array of past SMB values to compute AR model
+	/*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 ii{0};ii<numvertices;ii++){smblist[ii]=smbvaluesautoregression[ii];}
+	}
+	/*If AR model timestep: compute SMB values on vertices from AR*/
+	else{
+		for(int v=0;v<numvertices;v++){
+			double autoregressionterm{0.0};
+			for(int ii{0};ii<arorder;ii++){autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];} //compute autoregressive term
+			smblist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin; //stochastic SMB value
+		}
+		/*Update autoregression smb values*/
+		IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
+		for(int ii{0};ii<numvertices;ii++){temparray[ii] = smblist[ii];} //first store newly computed smb values
+		for(int ii{0};ii<(arorder-1)*numvertices;ii++){temparray[ii+numvertices] = smbvaluesautoregression[ii];} //second shift older smb values
+		this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); //save updated autoregression values
+		xDelete<IssmDouble>(temparray); //cleanup
+	}
+	/*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){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26476)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26477)
@@ -176,4 +176,6 @@
 		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);
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 26476)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 26477)
@@ -396,4 +396,23 @@
 			/*Nothing to add*/
 			break;
+		case SMBautoregressionEnum:
+         /*Add parameters that are not in standard nbvertices format*/
+         parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_basins",SmbNumBasinsEnum));
+         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_order",SmbAutoregressiveOrderEnum));
+         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_initialtime",SmbAutoregressionInitialTimeEnum));
+         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_timestep",SmbAutoregressionTimestepEnum));
+         iomodel->FetchData(&transparam,&M,&N,"md.smb.beta0");
+         parameters->AddObject(new DoubleVecParam(SmbBeta0Enum,transparam,N));
+         xDelete<IssmDouble>(transparam);
+         iomodel->FetchData(&transparam,&M,&N,"md.smb.beta1");
+         parameters->AddObject(new DoubleVecParam(SmbBeta1Enum,transparam,N));
+         xDelete<IssmDouble>(transparam);
+         iomodel->FetchData(&transparam,&M,&N,"md.smb.phi");
+         parameters->AddObject(new DoubleMatParam(SmbPhiEnum,transparam,M,N));
+         xDelete<IssmDouble>(transparam);
+         iomodel->FetchData(&transparam,&M,&N,"md.smb.covmat");
+         parameters->AddObject(new DoubleMatParam(SmbCovmatEnum,transparam,M,N));
+         xDelete<IssmDouble>(transparam);
+         break;
 		case SMBgembEnum:
 			parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 26476)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 26477)
@@ -8,4 +8,5 @@
 #include "../modules.h"
 #include "../../classes/Inputs/TransientInput.h"
+#include "../../shared/Random/random.h"
 
 void SmbForcingx(FemModel* femmodel){/*{{{*/
@@ -144,4 +145,84 @@
 	}
 
+}/*}}}*/
+void SmbautoregressionInitx(FemModel* femmodel){/*{{{*/
+	/*Initialization step of Smbautoregressionx*/
+	int M,N,Nphi,arorder,numbasins;
+	IssmDouble starttime,tstep_ar,tinit_ar;
+	IssmDouble* beta0    = xNew<IssmDouble>(numbasins);
+	IssmDouble* beta1    = xNew<IssmDouble>(numbasins);
+	IssmDouble* phi      = xNew<IssmDouble>(numbasins*arorder);
+	IssmDouble* covmat   = xNew<IssmDouble>(numbasins*numbasins);
+	femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+   femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
+	femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
+	femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
+	femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
+	femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
+	femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
+	femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
+	femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
+	/*AR model spin-up*/
+	int nspin{2*arorder+5}; //number of spin-up steps to be executed
+	IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin); //sample of basin-specific noise at each spinup step
+	for(int ii{0};ii<nspin;ii++){
+		IssmDouble* temparray = xNew<IssmDouble>(numbasins);
+		multivariateNormal(&temparray,numbasins,0.0,covmat);
+		for(int jj{0};jj<numbasins;jj++){noisespin[ii*numbasins+jj]=temparray[jj];}
+		xDelete<IssmDouble>(temparray);
+	}
+	/*Initialize SmbValuesAutoregressionEnum*/
+	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);
+	}
+	/*Cleanup*/
+	xDelete<IssmDouble>(beta0);
+	xDelete<IssmDouble>(beta1);
+	xDelete<IssmDouble>(phi);
+	xDelete<IssmDouble>(noisespin);
+	xDelete<IssmDouble>(covmat);
+}/*}}}*/
+void Smbautoregressionx(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,SmbAutoregressionTimestepEnum);
+	/*Initialize module at first time step*/
+	if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
+	/*Determine if this is a time step for the AR model*/
+	bool isstepforar{false};
+   if((std::fmod(time,tstep_ar)<std::fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt){isstepforar = true;}
+	/*Load parameters*/
+	int M,N,Nphi,arorder,numbasins;
+	femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
+	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* 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);
+	femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
+	femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
+	femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
+	IssmDouble telapsed_ar = time-tinit_ar; //time elapsed with respect to AR model initial time
+	/*Before looping through elements: compute noise term specific to each basin from covmat*/
+	multivariateNormal(&noiseterms,numbasins,0.0,covmat);
+	/*Loop over each element to compute SMB at vertices*/
+	for(Object* &object:femmodel->elements->objects){
+		Element* element      = xDynamicCast<Element*>(object); //generate element object
+		element->Smbautoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms);
+	}
+	/*Cleanup*/
+	xDelete<IssmDouble>(beta0);
+	xDelete<IssmDouble>(beta1);
+	xDelete<IssmDouble>(phi);
+	xDelete<IssmDouble>(noiseterms);
+	xDelete<IssmDouble>(covmat);
 }/*}}}*/
 void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 26476)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 26477)
@@ -13,4 +13,6 @@
 void SmbGradientsx(FemModel* femmodel);
 void SmbGradientsElax(FemModel* femmodel);
+void SmbautoregressionInitx(FemModel* femmodel);
+void Smbautoregressionx(FemModel* femmodel);
 void Delta18oParameterizationx(FemModel* femmodel);
 void MungsmtpParameterizationx(FemModel* femmodel);
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26476)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26477)
@@ -409,5 +409,11 @@
 	SmbAccurefEnum,
 	SmbAdThreshEnum,
+	SmbAutoregressionInitialTimeEnum,
+   SmbAutoregressionTimestepEnum,
+   SmbAutoregressiveOrderEnum,
 	SmbAveragingEnum,
+	SmbBeta0Enum,
+   SmbBeta1Enum,
+   SmbCovmatEnum,
 	SmbDesfacEnum,
 	SmbDpermilEnum,
@@ -438,6 +444,8 @@
 	SmbIsturbulentfluxEnum,
 	SmbKEnum,
+	SmbNumBasinsEnum,
 	SmbNumRequestedOutputsEnum,
 	SmbPfacEnum,
+	SmbPhiEnum,
 	SmbRdlEnum,
 	SmbRequestedOutputsEnum,
@@ -863,4 +871,5 @@
 	SmbAdiffiniEnum,
 	SmbAiniEnum,
+	SmbBasinsIdEnum,
 	SmbBMaxEnum,
 	SmbBMinEnum,
@@ -949,4 +958,5 @@
 	SmbTmeanEnum,
 	SmbTzEnum,
+	SmbValuesAutoregressionEnum,
 	SmbVEnum,
 	SmbVmeanEnum,
@@ -1414,4 +1424,5 @@
 	SamplingSolutionEnum,
 	SIAApproximationEnum,
+	SMBautoregressionEnum,
 	SMBcomponentsEnum,
 	SMBd18opddEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26476)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26477)
@@ -417,5 +417,11 @@
 		case SmbAccurefEnum : return "SmbAccuref";
 		case SmbAdThreshEnum : return "SmbAdThresh";
+		case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime";
+		case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep";
+		case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder";
 		case SmbAveragingEnum : return "SmbAveraging";
+		case SmbBeta0Enum : return "SmbBeta0";
+		case SmbBeta1Enum : return "SmbBeta1";
+		case SmbCovmatEnum : return "SmbCovmat";
 		case SmbDesfacEnum : return "SmbDesfac";
 		case SmbDpermilEnum : return "SmbDpermil";
@@ -446,6 +452,8 @@
 		case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux";
 		case SmbKEnum : return "SmbK";
+		case SmbNumBasinsEnum : return "SmbNumBasins";
 		case SmbNumRequestedOutputsEnum : return "SmbNumRequestedOutputs";
 		case SmbPfacEnum : return "SmbPfac";
+		case SmbPhiEnum : return "SmbPhi";
 		case SmbRdlEnum : return "SmbRdl";
 		case SmbRequestedOutputsEnum : return "SmbRequestedOutputs";
@@ -542,5 +550,9 @@
 		case AdjointpEnum : return "Adjointp";
 		case AdjointxEnum : return "Adjointx";
+		case AdjointxBaseEnum : return "AdjointxBase";
+		case AdjointxShearEnum : return "AdjointxShear";
 		case AdjointyEnum : return "Adjointy";
+		case AdjointyBaseEnum : return "AdjointyBase";
+		case AdjointyShearEnum : return "AdjointyShear";
 		case AdjointzEnum : return "Adjointz";
 		case AirEnum : return "Air";
@@ -865,4 +877,5 @@
 		case SmbAdiffiniEnum : return "SmbAdiffini";
 		case SmbAiniEnum : return "SmbAini";
+		case SmbBasinsIdEnum : return "SmbBasinsId";
 		case SmbBMaxEnum : return "SmbBMax";
 		case SmbBMinEnum : return "SmbBMin";
@@ -950,4 +963,5 @@
 		case SmbTmeanEnum : return "SmbTmean";
 		case SmbTzEnum : return "SmbTz";
+		case SmbValuesAutoregressionEnum : return "SmbValuesAutoregression";
 		case SmbVEnum : return "SmbV";
 		case SmbVmeanEnum : return "SmbVmean";
@@ -1413,4 +1427,5 @@
 		case SamplingSolutionEnum : return "SamplingSolution";
 		case SIAApproximationEnum : return "SIAApproximation";
+		case SMBautoregressionEnum : return "SMBautoregression";
 		case SMBcomponentsEnum : return "SMBcomponents";
 		case SMBd18opddEnum : return "SMBd18opdd";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26476)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26477)
@@ -426,5 +426,11 @@
 	      else if (strcmp(name,"SmbAccuref")==0) return SmbAccurefEnum;
 	      else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
+	      else if (strcmp(name,"SmbAutoregressionInitialTime")==0) return SmbAutoregressionInitialTimeEnum;
+	      else if (strcmp(name,"SmbAutoregressionTimestep")==0) return SmbAutoregressionTimestepEnum;
+	      else if (strcmp(name,"SmbAutoregressiveOrder")==0) return SmbAutoregressiveOrderEnum;
 	      else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum;
+	      else if (strcmp(name,"SmbBeta0")==0) return SmbBeta0Enum;
+	      else if (strcmp(name,"SmbBeta1")==0) return SmbBeta1Enum;
+	      else if (strcmp(name,"SmbCovmat")==0) return SmbCovmatEnum;
 	      else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum;
 	      else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum;
@@ -455,6 +461,8 @@
 	      else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum;
 	      else if (strcmp(name,"SmbK")==0) return SmbKEnum;
+	      else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum;
 	      else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
 	      else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum;
+	      else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum;
 	      else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum;
 	      else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
@@ -498,5 +506,8 @@
 	      else if (strcmp(name,"ThermalNumRequestedOutputs")==0) return ThermalNumRequestedOutputsEnum;
 	      else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
-	      else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
+         else stage=5;
+   }
+   if(stage==5){
+	      if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
 	      else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
 	      else if (strcmp(name,"ThermalReltol")==0) return ThermalReltolEnum;
@@ -506,8 +517,5 @@
 	      else if (strcmp(name,"Time")==0) return TimeEnum;
 	      else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
-         else stage=5;
-   }
-   if(stage==5){
-	      if (strcmp(name,"TimesteppingCouplingTime")==0) return TimesteppingCouplingTimeEnum;
+	      else if (strcmp(name,"TimesteppingCouplingTime")==0) return TimesteppingCouplingTimeEnum;
 	      else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
 	      else if (strcmp(name,"TimesteppingInterpForcing")==0) return TimesteppingInterpForcingEnum;
@@ -554,5 +562,9 @@
 	      else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
 	      else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
+	      else if (strcmp(name,"AdjointxBase")==0) return AdjointxBaseEnum;
+	      else if (strcmp(name,"AdjointxShear")==0) return AdjointxShearEnum;
 	      else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
+	      else if (strcmp(name,"AdjointyBase")==0) return AdjointyBaseEnum;
+	      else if (strcmp(name,"AdjointyShear")==0) return AdjointyShearEnum;
 	      else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
 	      else if (strcmp(name,"Air")==0) return AirEnum;
@@ -617,5 +629,8 @@
 	      else if (strcmp(name,"DamageDbarOld")==0) return DamageDbarOldEnum;
 	      else if (strcmp(name,"DamageF")==0) return DamageFEnum;
-	      else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
+         else stage=6;
+   }
+   if(stage==6){
+	      if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
 	      else if (strcmp(name,"DepthBelowSurface")==0) return DepthBelowSurfaceEnum;
 	      else if (strcmp(name,"DeltaIceThickness")==0) return DeltaIceThicknessEnum;
@@ -629,8 +644,5 @@
 	      else if (strcmp(name,"Str")==0) return StrEnum;
 	      else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
-         else stage=6;
-   }
-   if(stage==6){
-	      if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
+	      else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
 	      else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
 	      else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
@@ -740,5 +752,8 @@
 	      else if (strcmp(name,"LevelsetObservation")==0) return LevelsetObservationEnum;
 	      else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum;
-	      else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
+         else stage=7;
+   }
+   if(stage==7){
+	      if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
 	      else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum;
 	      else if (strcmp(name,"MaskOceanLevelset")==0) return MaskOceanLevelsetEnum;
@@ -752,8 +767,5 @@
 	      else if (strcmp(name,"MaterialsRheologyEc")==0) return MaterialsRheologyEcEnum;
 	      else if (strcmp(name,"MaterialsRheologyEcbar")==0) return MaterialsRheologyEcbarEnum;
-         else stage=7;
-   }
-   if(stage==7){
-	      if (strcmp(name,"MaterialsRheologyEs")==0) return MaterialsRheologyEsEnum;
+	      else if (strcmp(name,"MaterialsRheologyEs")==0) return MaterialsRheologyEsEnum;
 	      else if (strcmp(name,"MaterialsRheologyEsbar")==0) return MaterialsRheologyEsbarEnum;
 	      else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
@@ -863,5 +875,8 @@
 	      else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum;
 	      else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
-	      else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
+         else stage=8;
+   }
+   if(stage==8){
+	      if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
 	      else if (strcmp(name,"SedimentHeadSubstep")==0) return SedimentHeadSubstepEnum;
 	      else if (strcmp(name,"SedimentHeadTransient")==0) return SedimentHeadTransientEnum;
@@ -875,8 +890,5 @@
 	      else if (strcmp(name,"SmbAccumulatedPrecipitation")==0) return SmbAccumulatedPrecipitationEnum;
 	      else if (strcmp(name,"SmbAccumulatedRain")==0) return SmbAccumulatedRainEnum;
-         else stage=8;
-   }
-   if(stage==8){
-	      if (strcmp(name,"SmbAccumulatedRefreeze")==0) return SmbAccumulatedRefreezeEnum;
+	      else if (strcmp(name,"SmbAccumulatedRefreeze")==0) return SmbAccumulatedRefreezeEnum;
 	      else if (strcmp(name,"SmbAccumulatedRunoff")==0) return SmbAccumulatedRunoffEnum;
 	      else if (strcmp(name,"SmbA")==0) return SmbAEnum;
@@ -886,4 +898,5 @@
 	      else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum;
 	      else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
+	      else if (strcmp(name,"SmbBasinsId")==0) return SmbBasinsIdEnum;
 	      else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
 	      else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
@@ -971,4 +984,5 @@
 	      else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum;
 	      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 if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum;
@@ -984,5 +998,8 @@
 	      else if (strcmp(name,"SolidearthExternalDisplacementNorthRate")==0) return SolidearthExternalDisplacementNorthRateEnum;
 	      else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum;
-	      else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
+         else stage=9;
+   }
+   if(stage==9){
+	      if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
 	      else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
 	      else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
@@ -998,8 +1015,5 @@
 	      else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
 	      else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
-         else stage=9;
-   }
-   if(stage==9){
-	      if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
+	      else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
 	      else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
 	      else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
@@ -1107,5 +1121,8 @@
 	      else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
 	      else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
-	      else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum;
+         else stage=10;
+   }
+   if(stage==10){
+	      if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum;
 	      else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
 	      else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
@@ -1121,8 +1138,5 @@
 	      else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
 	      else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
-         else stage=10;
-   }
-   if(stage==10){
-	      if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
+	      else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
 	      else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum;
 	      else if (strcmp(name,"Outputdefinition64")==0) return Outputdefinition64Enum;
@@ -1230,5 +1244,8 @@
 	      else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum;
 	      else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
-	      else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
+         else stage=11;
+   }
+   if(stage==11){
+	      if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
 	      else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
 	      else if (strcmp(name,"Dense")==0) return DenseEnum;
@@ -1244,8 +1261,5 @@
 	      else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
 	      else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
-         else stage=11;
-   }
-   if(stage==11){
-	      if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
+	      else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
 	      else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
 	      else if (strcmp(name,"Element")==0) return ElementEnum;
@@ -1353,5 +1367,8 @@
 	      else if (strcmp(name,"LoveLr")==0) return LoveLrEnum;
 	      else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum;
-	      else if (strcmp(name,"MINI")==0) return MINIEnum;
+         else stage=12;
+   }
+   if(stage==12){
+	      if (strcmp(name,"MINI")==0) return MINIEnum;
 	      else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
 	      else if (strcmp(name,"MantlePlumeGeothermalFlux")==0) return MantlePlumeGeothermalFluxEnum;
@@ -1367,8 +1384,5 @@
 	      else if (strcmp(name,"Matestar")==0) return MatestarEnum;
 	      else if (strcmp(name,"Matice")==0) return MaticeEnum;
-         else stage=12;
-   }
-   if(stage==12){
-	      if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
+	      else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
 	      else if (strcmp(name,"Mathydro")==0) return MathydroEnum;
 	      else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
@@ -1446,4 +1460,5 @@
 	      else if (strcmp(name,"SamplingSolution")==0) return SamplingSolutionEnum;
 	      else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
+	      else if (strcmp(name,"SMBautoregression")==0) return SMBautoregressionEnum;
 	      else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
 	      else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
@@ -1475,5 +1490,8 @@
 	      else if (strcmp(name,"SegInput")==0) return SegInputEnum;
 	      else if (strcmp(name,"Segment")==0) return SegmentEnum;
-	      else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
+         else stage=13;
+   }
+   if(stage==13){
+	      if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
 	      else if (strcmp(name,"Separate")==0) return SeparateEnum;
 	      else if (strcmp(name,"Seq")==0) return SeqEnum;
@@ -1490,8 +1508,5 @@
 	      else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
 	      else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
-         else stage=13;
-   }
-   if(stage==13){
-	      if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
+	      else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
 	      else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;
 	      else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
Index: /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 26476)
+++ /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 26477)
@@ -673,2 +673,18 @@
 
 } /*}}}*/
+void CholeskyRealPositiveDefinite(IssmDouble* Lchol, IssmDouble* A, int ndim) { /*{{{*/
+   /*CholeskyRealPositiveDefinite   computes lower triangular matrix of the Cholesky decomposition of A
+   Follows the Cholesky–Banachiewicz algorithm
+   Lchol should point to an IssmDouble* of same dimensions as A*/
+   for(int ii{0};ii<(ndim*ndim);ii++){Lchol[ii]=0;}; //ensure zero-initialization
+   for(int ii{0};ii<ndim;ii++){
+                for(int jj{0};jj<=ii;jj++){
+                        double sum{0.0};
+                        for(int kk{0};kk<jj;kk++) {
+                                sum += Lchol[ii*ndim+kk]*Lchol[jj*ndim+kk];
+                        }
+                        if(ii==jj){Lchol[ii*ndim+jj] = sqrt(A[ii*ndim+jj]-sum);}
+                        else{Lchol[ii*ndim+jj] = 1/Lchol[jj*ndim+jj] * (A[ii*ndim+jj]-sum);}
+                }
+        }
+} /*}}}*/
Index: /issm/trunk-jpl/src/c/shared/Matrix/matrix.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Matrix/matrix.h	(revision 26476)
+++ /issm/trunk-jpl/src/c/shared/Matrix/matrix.h	(revision 26477)
@@ -31,3 +31,4 @@
 void cellsplit(IssmDouble** pcell, int m, int i,IssmDouble scale);
 void cellecho(int numcells, int m, ...);
+void CholeskyRealPositiveDefinite(IssmDouble* Lchol, IssmDouble* A, int ndim);
 #endif //ifndef _MATRIXUTILS_H_
Index: /issm/trunk-jpl/src/c/shared/Random/random.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Random/random.cpp	(revision 26477)
+++ /issm/trunk-jpl/src/c/shared/Random/random.cpp	(revision 26477)
@@ -0,0 +1,62 @@
+/*!\file: random
+ * \brief random number generating functions
+ */ 
+
+/*Headers*/
+/*{{{*/
+#include <stdio.h>
+#include <sys/types.h>
+#include <math.h>
+#include <float.h>    /*  DBL_EPSILON  */
+#include <chrono>
+#include <cstdarg>
+#include <iostream>
+#include <random>
+
+#include "../Matrix/matrix.h"
+#include "../Exceptions/exceptions.h"
+#include "../MemOps/MemOps.h"
+#include "../io/io.h"
+/*}}}*/
+
+void univariateNormal(IssmDouble* prand, IssmDouble mean, IssmDouble sdev) { /*{{{*/
+	/*univariateNormal	generates a random value follwoing Normal distribution*/
+	unsigned seed = std::chrono::steady_clock::now().time_since_epoch().count(); //random seed using time_since_epoch
+   std::default_random_engine generator(seed); //generator of random numbers
+   std::normal_distribution<double> normdistri(mean,sdev); //Normal probability distribution
+	double tfunc;
+	*prand = normdistri(generator);
+} /*}}}*/
+void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix) { /*{{{*/
+   IssmDouble* sampleStandardNormal     = xNew<IssmDouble>(dim);
+   IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
+   IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
+   for(int ii{0};ii<dim;ii++){univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);}
+   CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
+   for(int ii{0};ii<dim;ii++){ //matrix by vector multiplication
+      double sum{0.0};
+      for(int jj{0};jj<dim;jj++){sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];} //entry-by-entry multiplication along matrix row
+      sampleMultivariateNormal[ii] = mean+sum; //assign value
+   }
+   *prand = sampleMultivariateNormal;
+   xDelete<IssmDouble>(sampleStandardNormal);
+   xDelete<IssmDouble>(Lchol);
+} /*}}}*/
+void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix) { /*{{{*/
+	IssmDouble* sampleStandardNormal     = xNew<IssmDouble>(dim);
+	IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
+	IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
+	for(int ii{0};ii<dim;ii++){univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);}
+	CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
+	for(int ii{0};ii<dim;ii++){ //matrix by vector multiplication
+		double sum{0.0};
+      for(int jj{0};jj<dim;jj++){sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];} //entry-by-entry multiplication along matrix row
+      sampleMultivariateNormal[ii] = mean[ii]+sum; //assign value
+	}
+	*prand = sampleMultivariateNormal;
+	xDelete<IssmDouble>(sampleStandardNormal);
+	xDelete<IssmDouble>(Lchol);
+} /*}}}*/
+
+
+
Index: /issm/trunk-jpl/src/c/shared/Random/random.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Random/random.h	(revision 26477)
+++ /issm/trunk-jpl/src/c/shared/Random/random.h	(revision 26477)
@@ -0,0 +1,12 @@
+/*!\file: random.h
+ * \brief prototypes for random.h
+ */ 
+
+#ifndef _RANDOM_H_
+#define _RANDOM_H_
+
+void univariateNormal(IssmDouble* prand, IssmDouble mean, IssmDouble sdev);
+void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix);
+void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix);
+
+#endif //ifndef _RANDOM_H_
Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp	(revision 26476)
+++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp	(revision 26477)
@@ -241,4 +241,5 @@
 		case 11: return SMBgradientscomponentsEnum;
 		case 12: return SMBsemicEnum;	
+		case 55: return SMBautoregressionEnum;
 		default: _error_("Marshalled SMB code \""<<enum_in<<"\" not supported yet");
 	}
Index: /issm/trunk-jpl/src/m/classes/SMBautoregression.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBautoregression.m	(revision 26477)
+++ /issm/trunk-jpl/src/m/classes/SMBautoregression.m	(revision 26477)
@@ -0,0 +1,130 @@
+%SMBautoregression Class definition
+%
+%   Usage:
+%      SMBautoregression=SMBautoregression();
+
+classdef SMBautoregression
+	properties (SetAccess=public)
+		num_basins        = 0;
+		beta0             = NaN;
+		beta1             = NaN;
+		ar_order          = 0;
+		ar_initialtime    = 0;
+		ar_timestep       = 0;
+		phi               = NaN;
+		covmat            = NaN;	
+		basin_id          = NaN;	
+		steps_per_step    = 1;
+		averaging         = 0;
+		requested_outputs = {};
+	end
+	methods
+		function self = SMBautoregression(varargin) % {{{
+			switch nargin
+				case 0
+					self=setdefaultparameters(self);
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function self = extrude(self,md) % {{{
+			%Nothing for now
+		end % }}}
+		function list = defaultoutputs(self,md) % {{{
+			list = {''};
+		end % }}}
+		function self = initialize(self,md) % {{{
+			if isnan(self.beta1)
+				self.beta1 = zeros(1,self.num_basins); %no trend in SMB
+				disp('      smb.beta1 (trend) not specified: value set to 0');
+			end
+			if (self.ar_order==0)
+				self.ar_order = 1; %dummy 1 value for autoregression
+				self.phi      = zeros(self.num_basins,self.ar_order); %autoregression coefficients all set to 0 
+				disp('      smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0');	
+			end
+			if (self.ar_initialtime==0)
+				self.ar_initialtime = md.timestepping.start_time; %autoregression model has no prescribed initial time
+				disp('      smb.ar_initialtime (initial time in the autoregressive model parameterization) not specified: set to md.timestepping.start_time');	
+			end
+			if (self.ar_timestep==0)
+				self.ar_timestep = md.timestepping.time_step; %autoregression model has no prescribed time step
+				disp('      smb.ar_timestep (timestep of autoregressive model) not specified: set to md.timestepping.time_step');	
+			end
+			if isnan(self.phi)
+				self.phi = zeros(self.num_basins,self.ar_order); %autoregression model of order 0 
+				disp('      smb.phi (lag coefficients) not specified: order of autoregressive model set to 0');
+			end
+			if isnan(self.covmat)
+				self.covmat = 1e-21*eye(self.num_basins); %no stochasticity and no covariance
+				disp('      smb.covmat not specified: stochasticity set to 0');
+			end
+		end % }}}
+		function self = setdefaultparameters(self) % {{{
+			self.ar_order    = 0.0; %autoregression model of order 0 
+		end % }}}
+		function md = checkconsistency(self,md,solution,analyses) % {{{
+
+			if ismember('MasstransportAnalysis',analyses),
+				md = checkfield(md,'fieldname','smb.num_basins','numel',1,'NaN',1,'Inf',1,'>',0);
+				md = checkfield(md,'fieldname','smb.basin_id','Inf',1,'>=',0,'<=',md.smb.num_basins,'size',[md.mesh.numberofelements,1]);
+				md = checkfield(md,'fieldname','smb.beta0','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); %scheme fails if passed as column vector
+				md = checkfield(md,'fieldname','smb.beta1','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); %scheme fails if passed as column vector
+				md = checkfield(md,'fieldname','smb.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0);
+				md = checkfield(md,'fieldname','smb.ar_initialtime','numel',1,'NaN',1,'Inf',1); 
+				md = checkfield(md,'fieldname','smb.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep
+				md = checkfield(md,'fieldname','smb.phi','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]);
+				md = checkfield(md,'fieldname','smb.covmat','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.num_basins]);
+			end
+			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
+			md = checkfield(md,'fieldname','smb.averaging','numel',[1],'values',[0 1 2]);
+			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
+		end % }}}
+		function disp(self) % {{{
+			disp(sprintf('   surface forcings parameters:'));
+			fielddisplay(self,'num_basins','number of different basins [unitless]');
+			fielddisplay(self,'basin_id','basin number assigned to each element [unitless]');
+			fielddisplay(self,'beta0','basin-specific intercept values [m ice eq./yr]');
+			fielddisplay(self,'beta1','basin-specific trend values [m ice eq. yr^(-2)]');
+			fielddisplay(self,'ar_order','order of the autoregressive model [unitless]');
+			fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]');
+			fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]');
+			fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]');
+			fielddisplay(self,'covmat','inter-basin covariance matrix for multivariate normal noise at each time step [m^2 ice eq. yr^(-2)]');
+			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%51s  0: Arithmetic (default)',' '));
+			disp(sprintf('%51s  1: Geometric',' '));
+			disp(sprintf('%51s  2: Harmonic',' '));
+			fielddisplay(self,'requested_outputs','additional outputs requested');
+
+		end % }}}
+		function marshall(self,prefix,md,fid) % {{{
+
+			yts=md.constants.yts;
+
+			WriteData(fid,prefix,'name','md.smb.model','data',55,'format','Integer');
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','num_basins','format','Integer');
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_order','format','Integer');
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_initialtime','format','Double','scale',yts);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_timestep','format','Double','scale',yts);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','basin_id','data',self.basin_id-1,'name','md.smb.basin_id','format','IntMat','mattype',2); %0-indexed
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta0','format','DoubleMat','name','md.smb.beta0','scale',1./yts,'yts',md.constants.yts);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',md.constants.yts);
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',md.constants.yts); 
+			WriteData(fid,prefix,'object',self,'class','smb','fieldname','covmat','format','DoubleMat','name','md.smb.covmat','scale',1./(yts^2),'yts',md.constants.yts);
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
+
+			%process requested outputs
+			outputs = self.requested_outputs;
+			pos  = find(ismember(outputs,'default'));
+			if ~isempty(pos),
+				outputs(pos) = [];                         %remove 'default' from outputs
+				outputs      = [outputs defaultoutputs(self,md)]; %add defaults
+			end
+			WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray');
+
+		end % }}}
+	end
+end