Index: issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 26603)
+++ issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 26604)
@@ -12,37 +12,8 @@
 
 void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
-	
-	/*intermediary: */
 	int finiteelement;
-	int         code,vector_layout;
-	IssmDouble *spcdata = NULL;
-	int         M,N;
-
-	/*Get finite element type for this analysis*/
 	iomodel->FindConstant(&finiteelement,"md.levelset.fe");
-
-	/*First of, find the record for the enum, and get code  of data type: */
-	iomodel->SetFilePointerToData(&code, &vector_layout,"md.levelset.spclevelset");
-	if(code!=7)_error_("expecting a IssmDouble vector for constraints md.levelset.spclevelset");
-	if(vector_layout!=1)_error_("expecting a nodal vector for constraints md.levelset.spclevelset");
-
-	/*Fetch vector:*/
-	iomodel->FetchData(&spcdata,&M,&N,"md.levelset.spclevelset");
-
-	/*Call IoModelToConstraintsx*/
-	if(N>1){
-		/*If it is a time series, most likely we are forcing the ice front position and do not want to have a Dynamic Constraint*/
-		_assert_(M==iomodel->numberofvertices+1);
-		IoModelToConstraintsx(constraints,iomodel,spcdata,M,N,LevelsetAnalysisEnum,finiteelement);
-	}
-	else{
-		/*This is not a time series, we probably have calving on, we need the levelset constraints to update as the levelset moves*/
-		_assert_(N==1);
-		_assert_(M==iomodel->numberofvertices);
-		IoModelToDynamicConstraintsx(constraints,iomodel,spcdata,M,N,LevelsetAnalysisEnum,finiteelement);
-	}
-
-	/*Clean up*/
-	xDelete<IssmDouble>(spcdata);
+	IoModelToDynamicConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,finiteelement);
+	//IoModelToConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,finiteelement);
 }
 /*}}}*/
@@ -83,9 +54,15 @@
 
 	/*Get moving front parameters*/
+	bool isstochastic;
 	int  calvinglaw;
 	iomodel->FindConstant(&calvinglaw,"md.calving.law");
+	iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing");
 	switch(calvinglaw){
 		case DefaultCalvingEnum:
 			iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",CalvingCalvingrateEnum);
+			if(isstochastic){
+				iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
+				iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",BaselineCalvingCalvingrateEnum);
+			}
 			break;
 		case CalvingLevermannEnum:
Index: issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26603)
+++ issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26604)
@@ -99,22 +99,75 @@
    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;
+void       Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type){/*{{{*/
+   
+   const int numvertices = this->GetNumberOfVertices();
+   int         basinid,M,N,arenum_type,basinenum_type,noiseenum_type,outenum_type;
+   IssmDouble  beta0_basin,beta1_basin,noiseterm;
    IssmDouble* phi_basin  = xNew<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);
+   Input*      noiseterm_input      = NULL;
+
+   /*Get field-specific enums*/
+	switch(enum_type){
+      case(SMBautoregressionEnum):
+			arenum_type    = SmbValuesAutoregressionEnum;
+			basinenum_type = SmbBasinsIdEnum;
+			noiseenum_type = SmbAutoregressionNoiseEnum;
+         outenum_type   = SmbMassBalanceEnum;
+			break;
+		case(FrontalForcingsRignotAutoregressionEnum):
+         arenum_type    = ThermalforcingValuesAutoregressionEnum;
+         basinenum_type = FrontalForcingsBasinIdEnum;
+         noiseenum_type = ThermalforcingAutoregressionNoiseEnum;
+         outenum_type   = FrontalForcingsThermalForcingEnum;
+			break;
+	}
+
+	/*Get noise and autoregressive terms*/
+	this->GetInputValue(&basinid,basinenum_type);
+	if(isfieldstochastic){
+      noiseterm_input = this->GetInput(noiseenum_type);
+      Gauss* gauss = this->NewGauss();
+      noiseterm_input->GetInputValue(&noiseterm,gauss);
+      delete gauss;
+   }
+	else noiseterm = 0.0;
+   this->inputs->GetArray(arenum_type,this->lid,&valuesautoregression,&M);
+
+
+	/*
+         this->GetInputValue(&basinid,SmbBasinsIdEnum);
+         if(isfieldstochastic){
+            noiseterm_input = this->GetInput(SmbAutoregressionNoiseEnum);
+            Gauss* gauss = this->NewGauss();
+            noiseterm_input->GetInputValue(&noiseterm,gauss);
+            delete gauss;
+         }
+         else noiseterm = 0.0;
+         this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
+         arenum_type  = SmbValuesAutoregressionEnum;
+         outenum_type = SmbMassBalanceEnum;
+         break;
+      case(FrontalForcingsRignotAutoregressionEnum):
+         this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum);
+         if(isfieldstochastic){
+            noiseterm_input = this->GetInput(ThermalforcingAutoregressionNoiseEnum);
+            Gauss* gauss = this->NewGauss();
+            noiseterm_input->GetInputValue(&noiseterm,gauss);
+            delete gauss;
+         }
+         else noiseterm = 0.0;
+         this->inputs->GetArray(ThermalforcingValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
+         arenum_type  = ThermalforcingValuesAutoregressionEnum;
+         outenum_type = FrontalForcingsThermalForcingEnum;
+         break;
+   }
+	*/
+
+	/*Get basin coefficients*/
    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*/
@@ -131,24 +184,22 @@
 
          /*Stochastic variable value*/
-         varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin;
+         varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noiseterm;
       }
-      /*Update autoregression TF values*/
+      /*Update autoregression values*/
       IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
       /*Assign newest values and shift older values*/
       for(int i=0;i<numvertices;i++) temparray[i] = varlist[i];
       for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = valuesautoregression[i];
-      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);
+      this->inputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder);
       xDelete<IssmDouble>(temparray);
    }
 
    /*Add input to element*/
-   if(enum_type==SMBautoregressionEnum)                   this->AddInput(SmbMassBalanceEnum,varlist,P1Enum);
-   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum);
+   this->AddInput(outenum_type,varlist,P1Enum);
 
    /*Cleanup*/
    xDelete<IssmDouble>(phi_basin);
    xDelete<IssmDouble>(varlist);
-   xDelete<IssmDouble>(valuesautoregression);
+   xDelete<IssmDouble>(valuesautoregression); 
 }/*}}}*/
 void       Element::ComputeLambdaS(){/*{{{*/
@@ -3778,5 +3829,4 @@
 	IssmDouble teValue=1.0;
 	IssmDouble aValue=0.0;
-	IssmDouble dulwrfValue=0.0;
 	IssmDouble szaValue=0.0;
 	IssmDouble cotValue=0.0;
@@ -3785,5 +3835,4 @@
 	IssmDouble dt,time,smb_dt;
 	int        aIdx=0;
-	int        eIdx=0;
 	int        denIdx=0;
 	int        dsnowIdx=0;
@@ -3815,9 +3864,7 @@
 	bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux;
 	bool isconstrainsurfaceT=false;
-	bool isdeltaLWup=false;
 	IssmDouble init_scaling=0.0;
 	IssmDouble thermo_scaling=1.0;
 	IssmDouble adThresh=1023.0;
-	IssmDouble teThresh=10;
 	/*}}}*/
 	/*Output variables:{{{ */
@@ -3869,5 +3916,4 @@
 	parameters->FindParam(&smb_dt,SmbDtEnum);                     /*time period for the smb solution,  usually smaller than the glaciological dt*/
 	parameters->FindParam(&aIdx,SmbAIdxEnum);
-	parameters->FindParam(&eIdx,SmbEIdxEnum);
 	parameters->FindParam(&denIdx,SmbDenIdxEnum);
 	parameters->FindParam(&swIdx,SmbSwIdxEnum);
@@ -3886,9 +3932,7 @@
 	parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum);
 	parameters->FindParam(&isconstrainsurfaceT,SmbIsconstrainsurfaceTEnum);
-	parameters->FindParam(&isdeltaLWup,SmbIsdeltaLWupEnum);
 	parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum);
 	parameters->FindParam(&thermo_scaling,SmbThermoDeltaTScalingEnum);
 	parameters->FindParam(&adThresh,SmbAdThreshEnum);
-	parameters->FindParam(&teThresh,SmbTeThreshEnum);
 	/*}}}*/
 	/*Retrieve inputs: {{{*/
@@ -4080,5 +4124,4 @@
 	Input *teValue_input= this->GetInput(SmbTeValueEnum,timeinputs); _assert_(teValue_input);
 	Input *aValue_input= this->GetInput(SmbAValueEnum,timeinputs); _assert_(aValue_input);
-	Input *dulwrfValue_input= this->GetInput(SmbDulwrfValueEnum,timeinputs); _assert_(dulwrfValue_input);
 	Input *szaValue_input= this->GetInput(SmbSzaValueEnum,timeinputs); _assert_(szaValue_input);
 	Input *cotValue_input= this->GetInput(SmbCotValueEnum,timeinputs); _assert_(cotValue_input);
@@ -4096,5 +4139,4 @@
 	pAir_input->GetInputValue(&pAir,gauss);  // screen level air pressure [Pa]
 	teValue_input->GetInputValue(&teValue,gauss);  // Emissivity [0-1]
-	dulwrfValue_input->GetInputValue(&dulwrfValue,gauss);  // LWup perturbation [W m-2]
 	aValue_input->GetInputValue(&aValue,gauss);  // Albedo [0 1]
 	szaValue_input->GetInputValue(&szaValue,gauss);  // Solar Zenith Angle [degree]
@@ -4122,5 +4164,5 @@
 	}
 	/*Thermal profile computation:*/
-	if(isthermal)thermo(&EC, &T, &ulw, re, dz, d, swf, dlw, Ta, V, eAir, pAir, eIdx, teValue, dulwrfValue, teThresh, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid(),isconstrainsurfaceT,isdeltaLWup);
+	if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid(),isconstrainsurfaceT);
 
 	/*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell.
Index: issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26603)
+++ issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26604)
@@ -69,5 +69,5 @@
 		bool               AnyFSet(void);
 		void					 AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type);
-      void               Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type);
+      void               Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type);
 		void               ComputeLambdaS(void);
 		void               ComputeNewDamage();
Index: issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp	(revision 26603)
+++ issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp	(revision 26604)
@@ -80,5 +80,6 @@
    /*Load parameters*/
 	bool isstochastic;
-   int M,N,Nphi,arorder,numbasins,my_rank;
+   bool istfstochastic = false;
+	int M,N,Nphi,arorder,numbasins,my_rank;
    femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
    femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
@@ -87,5 +88,4 @@
    IssmDouble* beta1      = NULL;
    IssmDouble* phi        = NULL;
-   IssmDouble* noiseterms = xNewZeroInit<IssmDouble>(numbasins);
 
 	femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
@@ -94,5 +94,4 @@
    femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
 
-	/*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/
 	femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
 	if(isstochastic){
@@ -102,8 +101,7 @@
 		femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
 		for(int i=0;i<numstochasticfields;i++){
-			if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum){
-				femmodel->parameters->FindParam(&noiseterms,&M,ThermalforcingAutoregressionNoiseEnum);  _assert_(M==numbasins);
-			}
+			if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum) istfstochastic = true;
 		}
+		xDelete<int>(stochasticfields);
 	}
    /*Time elapsed with respect to AR model initial time*/
@@ -113,5 +111,5 @@
    for(Object* &object:femmodel->elements->objects){
       Element* element = xDynamicCast<Element*>(object);
-      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,FrontalForcingsRignotAutoregressionEnum);
+      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum);
    }
 
@@ -120,4 +118,3 @@
    xDelete<IssmDouble>(beta1);
    xDelete<IssmDouble>(phi);
-   xDelete<IssmDouble>(noiseterms);
 }/*}}}*/
Index: issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 26603)
+++ issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 26604)
@@ -502,7 +502,8 @@
 	parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
 	if(isstochasticforcing){
-		int num_fields;
+		int num_fields,stochastic_dim;
 		char** fields;
 		parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum));
+		parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.defaultdimension",StochasticForcingDefaultDimensionEnum));
 		iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields");
 		if(num_fields<1) _error_("no stochasticforcing fields found");
@@ -515,6 +516,5 @@
 		parameters->AddObject(new IntVecParam(StochasticForcingFieldsEnum,stochasticforcing_enums,num_fields));
 		xDelete<int>(stochasticforcing_enums);
-
-      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.randomflag",StochasticForcingRandomflagEnum));
+		parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.randomflag",StochasticForcingRandomflagEnum));
 		iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.dimensions");
 		parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N));
Index: issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 26603)
+++ issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 26604)
@@ -36,5 +36,5 @@
       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*/
+		/*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
       IssmDouble* temparray = NULL;
       multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed);
@@ -46,15 +46,66 @@
 	int i=0;
    for(int j=0;j<numfields;j++){
-      int enum_type;
+      int dimenum_type,noiseenum_type;
       IssmDouble* noisefield = xNew<IssmDouble>(dimensions[j]);
       for(int k=0;k<dimensions[j];k++){
          noisefield[k]=noiseterms[i+k];
       }
-      
-      if(fields[j]==SMBautoregressionEnum)                        enum_type = SmbAutoregressionNoiseEnum;
-		else if(fields[j]==FrontalForcingsRignotAutoregressionEnum) enum_type = ThermalforcingAutoregressionNoiseEnum;
-		else _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.\n"); 
-      femmodel->parameters->SetParam(noisefield,dimensions[j],enum_type);
-      i=i+dimensions[j];
+     
+		int dimensionid;
+
+		/*Deal with the autoregressive models*/
+		if(fields[j]==SMBautoregressionEnum || fields[j]==FrontalForcingsRignotAutoregressionEnum){
+			switch(fields[j]){
+				case SMBautoregressionEnum:
+					dimenum_type   = SmbBasinsIdEnum;
+					noiseenum_type = SmbAutoregressionNoiseEnum;
+					break;
+				case FrontalForcingsRignotAutoregressionEnum:
+					dimenum_type   = FrontalForcingsBasinIdEnum;
+					noiseenum_type = ThermalforcingAutoregressionNoiseEnum;
+					break;
+			}
+			for(Object* &object:femmodel->elements->objects){
+            Element* element = xDynamicCast<Element*>(object);
+            int numvertices  = element->GetNumberOfVertices();
+            IssmDouble* noise_element = xNew<IssmDouble>(numvertices);
+            element->GetInputValue(&dimensionid,dimenum_type);
+            for(int i=0;i<numvertices;i++) noise_element[i] = noisefield[dimensionid];
+            element->AddInput(noiseenum_type,noise_element,P0Enum);
+            xDelete<IssmDouble>(noise_element);
+			}
+		}
+		else{
+			switch(fields[j]){
+				case SMBautoregressionEnum:
+				case FrontalForcingsRignotAutoregressionEnum:
+					/*Already done above*/
+					break;
+				case DefaultCalvingEnum:
+					/*Delete CalvingCalvingrateEnum at previous time step (required if it is transient)*/
+					femmodel->inputs->DeleteInput(CalvingCalvingrateEnum);
+					for(Object* &object:femmodel->elements->objects){
+						Element* element = xDynamicCast<Element*>(object);
+						int numvertices  = element->GetNumberOfVertices();
+						IssmDouble baselinecalvingrate;
+						IssmDouble calvingrate_tot[numvertices];
+						Input* baselinecalvingrate_input  = NULL;
+						baselinecalvingrate_input = element->GetInput(BaselineCalvingCalvingrateEnum); _assert_(baselinecalvingrate_input);
+						element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
+						Gauss* gauss = element->NewGauss();
+						for(int i=0;i<numvertices;i++){
+							gauss->GaussVertex(i);
+							baselinecalvingrate_input->GetInputValue(&baselinecalvingrate,gauss);
+							calvingrate_tot[i] = max(0.0,baselinecalvingrate+noisefield[dimensionid]);
+						}
+						element->AddInput(CalvingCalvingrateEnum,&calvingrate_tot[0],P1DGEnum);
+						delete gauss;
+					}
+					break;
+				default:
+					_error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.");
+			}
+		}
+		i=i+dimensions[j];
       xDelete<IssmDouble>(noisefield);
    }
Index: issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 26603)
+++ issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 26604)
@@ -197,4 +197,5 @@
 	/*Load parameters*/
 	bool isstochastic;
+	bool issmbstochastic = false;
 	int M,N,Nphi,arorder,numbasins,my_rank;
 	femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
@@ -210,6 +211,4 @@
 	femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
 
-	/*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/
-   IssmDouble* noiseterms = xNewZeroInit<IssmDouble>(numbasins);
 	femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
    if(isstochastic){
@@ -219,7 +218,5 @@
       femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
       for(int i=0;i<numstochasticfields;i++){
-         if(stochasticfields[i]==SMBautoregressionEnum){
-				femmodel->parameters->FindParam(&noiseterms,&M,SmbAutoregressionNoiseEnum);  _assert_(M==numbasins);
-			}
+         if(stochasticfields[i]==SMBautoregressionEnum) issmbstochastic = true;
 		}
 		xDelete<int>(stochasticfields);
@@ -231,5 +228,5 @@
 	for(Object* &object:femmodel->elements->objects){
 		Element* element = xDynamicCast<Element*>(object);
-		element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,SMBautoregressionEnum);
+		element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum);
 	}
 
@@ -238,5 +235,4 @@
 	xDelete<IssmDouble>(beta1);
 	xDelete<IssmDouble>(phi);
-	xDelete<IssmDouble>(noiseterms);
 }/*}}}*/
 void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
Index: issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26603)
+++ issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26604)
@@ -398,4 +398,5 @@
 syn keyword cConstant SolidearthSettingsOceanAreaScalingEnum
 syn keyword cConstant StochasticForcingCovarianceEnum
+syn keyword cConstant StochasticForcingDefaultDimensionEnum
 syn keyword cConstant StochasticForcingDimensionsEnum
 syn keyword cConstant StochasticForcingFieldsEnum
@@ -428,5 +429,4 @@
 syn keyword cConstant SmbAdThreshEnum
 syn keyword cConstant SmbAutoregressionInitialTimeEnum
-syn keyword cConstant SmbAutoregressionNoiseEnum
 syn keyword cConstant SmbAutoregressionTimestepEnum
 syn keyword cConstant SmbAutoregressiveOrderEnum
@@ -504,5 +504,4 @@
 syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum
 syn keyword cConstant StressbalanceShelfDampeningEnum
-syn keyword cConstant ThermalforcingAutoregressionNoiseEnum
 syn keyword cConstant ThermalIsdrainicecolumnEnum
 syn keyword cConstant ThermalIsdynamicbasalspcEnum
@@ -599,4 +598,5 @@
 syn keyword cConstant BasalStressEnum
 syn keyword cConstant BaseEnum
+syn keyword cConstant BaselineCalvingCalvingrateEnum
 syn keyword cConstant BaseOldEnum
 syn keyword cConstant BaseSlopeXEnum
@@ -613,4 +613,5 @@
 syn keyword cConstant BottomPressureOldEnum
 syn keyword cConstant CalvingCalvingrateEnum
+syn keyword cConstant CalvingCalvingrateNoiseEnum
 syn keyword cConstant CalvingHabFractionEnum
 syn keyword cConstant CalvingMeltingrateEnum
@@ -895,4 +896,5 @@
 syn keyword cConstant SmbAdiffiniEnum
 syn keyword cConstant SmbAiniEnum
+syn keyword cConstant SmbAutoregressionNoiseEnum
 syn keyword cConstant SmbBasinsIdEnum
 syn keyword cConstant SmbBMaxEnum
@@ -997,4 +999,5 @@
 syn keyword cConstant SolidearthExternalDisplacementUpRateEnum
 syn keyword cConstant SolidearthExternalGeoidRateEnum
+syn keyword cConstant StochasticForcingDefaultIdEnum
 syn keyword cConstant StrainRateeffectiveEnum
 syn keyword cConstant StrainRateparallelEnum
@@ -1032,4 +1035,5 @@
 syn keyword cConstant TemperaturePicardEnum
 syn keyword cConstant TemperatureSEMICEnum
+syn keyword cConstant ThermalforcingAutoregressionNoiseEnum
 syn keyword cConstant ThermalforcingValuesAutoregressionEnum
 syn keyword cConstant ThermalSpctemperatureEnum
Index: issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26603)
+++ issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26604)
@@ -392,4 +392,5 @@
 	SolidearthSettingsOceanAreaScalingEnum,
 	StochasticForcingCovarianceEnum,
+	StochasticForcingDefaultDimensionEnum,
 	StochasticForcingDimensionsEnum,
 	StochasticForcingFieldsEnum,
@@ -422,5 +423,4 @@
 	SmbAdThreshEnum,
 	SmbAutoregressionInitialTimeEnum,
-   SmbAutoregressionNoiseEnum,
 	SmbAutoregressionTimestepEnum,
    SmbAutoregressiveOrderEnum,
@@ -498,5 +498,4 @@
 	StressbalanceRiftPenaltyThresholdEnum,
 	StressbalanceShelfDampeningEnum,
-   ThermalforcingAutoregressionNoiseEnum,
 	ThermalIsdrainicecolumnEnum,
 	ThermalIsdynamicbasalspcEnum,
@@ -595,4 +594,5 @@
 	BasalStressEnum,
 	BaseEnum,
+	BaselineCalvingCalvingrateEnum,
 	BaseOldEnum,
 	BaseSlopeXEnum,
@@ -609,4 +609,5 @@
 	BottomPressureOldEnum,
 	CalvingCalvingrateEnum,
+	CalvingCalvingrateNoiseEnum,
 	CalvingHabFractionEnum,
 	CalvingMeltingrateEnum,
@@ -891,4 +892,5 @@
 	SmbAdiffiniEnum,
 	SmbAiniEnum,
+   SmbAutoregressionNoiseEnum,
 	SmbBasinsIdEnum,
 	SmbBMaxEnum,
@@ -994,4 +996,5 @@
 	SolidearthExternalDisplacementUpRateEnum,
 	SolidearthExternalGeoidRateEnum,
+	StochasticForcingDefaultIdEnum,
 	StrainRateeffectiveEnum,
 	StrainRateparallelEnum,
@@ -1029,4 +1032,5 @@
 	TemperaturePicardEnum,
 	TemperatureSEMICEnum,
+   ThermalforcingAutoregressionNoiseEnum,
 	ThermalforcingValuesAutoregressionEnum,	
 	ThermalSpctemperatureEnum,
Index: issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26603)
+++ issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26604)
@@ -400,4 +400,5 @@
 		case SolidearthSettingsOceanAreaScalingEnum : return "SolidearthSettingsOceanAreaScaling";
 		case StochasticForcingCovarianceEnum : return "StochasticForcingCovariance";
+		case StochasticForcingDefaultDimensionEnum : return "StochasticForcingDefaultDimension";
 		case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions";
 		case StochasticForcingFieldsEnum : return "StochasticForcingFields";
@@ -430,5 +431,4 @@
 		case SmbAdThreshEnum : return "SmbAdThresh";
 		case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime";
-		case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
 		case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep";
 		case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder";
@@ -506,5 +506,4 @@
 		case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold";
 		case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening";
-		case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise";
 		case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn";
 		case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc";
@@ -601,4 +600,5 @@
 		case BasalStressEnum : return "BasalStress";
 		case BaseEnum : return "Base";
+		case BaselineCalvingCalvingrateEnum : return "BaselineCalvingCalvingrate";
 		case BaseOldEnum : return "BaseOld";
 		case BaseSlopeXEnum : return "BaseSlopeX";
@@ -615,4 +615,5 @@
 		case BottomPressureOldEnum : return "BottomPressureOld";
 		case CalvingCalvingrateEnum : return "CalvingCalvingrate";
+		case CalvingCalvingrateNoiseEnum : return "CalvingCalvingrateNoise";
 		case CalvingHabFractionEnum : return "CalvingHabFraction";
 		case CalvingMeltingrateEnum : return "CalvingMeltingrate";
@@ -897,4 +898,5 @@
 		case SmbAdiffiniEnum : return "SmbAdiffini";
 		case SmbAiniEnum : return "SmbAini";
+		case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
 		case SmbBasinsIdEnum : return "SmbBasinsId";
 		case SmbBMaxEnum : return "SmbBMax";
@@ -999,4 +1001,5 @@
 		case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate";
 		case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate";
+		case StochasticForcingDefaultIdEnum : return "StochasticForcingDefaultId";
 		case StrainRateeffectiveEnum : return "StrainRateeffective";
 		case StrainRateparallelEnum : return "StrainRateparallel";
@@ -1034,4 +1037,5 @@
 		case TemperaturePicardEnum : return "TemperaturePicard";
 		case TemperatureSEMICEnum : return "TemperatureSEMIC";
+		case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise";
 		case ThermalforcingValuesAutoregressionEnum : return "ThermalforcingValuesAutoregression";
 		case ThermalSpctemperatureEnum : return "ThermalSpctemperature";
Index: issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26603)
+++ issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26604)
@@ -409,4 +409,5 @@
 	      else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum;
 	      else if (strcmp(name,"StochasticForcingCovariance")==0) return StochasticForcingCovarianceEnum;
+	      else if (strcmp(name,"StochasticForcingDefaultDimension")==0) return StochasticForcingDefaultDimensionEnum;
 	      else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum;
 	      else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum;
@@ -439,5 +440,4 @@
 	      else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
 	      else if (strcmp(name,"SmbAutoregressionInitialTime")==0) return SmbAutoregressionInitialTimeEnum;
-	      else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum;
 	      else if (strcmp(name,"SmbAutoregressionTimestep")==0) return SmbAutoregressionTimestepEnum;
 	      else if (strcmp(name,"SmbAutoregressiveOrder")==0) return SmbAutoregressiveOrderEnum;
@@ -518,5 +518,4 @@
 	      else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum;
 	      else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum;
-	      else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum;
 	      else if (strcmp(name,"ThermalIsdrainicecolumn")==0) return ThermalIsdrainicecolumnEnum;
 	      else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
@@ -613,4 +612,5 @@
 	      else if (strcmp(name,"BasalStress")==0) return BasalStressEnum;
 	      else if (strcmp(name,"Base")==0) return BaseEnum;
+	      else if (strcmp(name,"BaselineCalvingCalvingrate")==0) return BaselineCalvingCalvingrateEnum;
 	      else if (strcmp(name,"BaseOld")==0) return BaseOldEnum;
 	      else if (strcmp(name,"BaseSlopeX")==0) return BaseSlopeXEnum;
@@ -627,10 +627,11 @@
 	      else if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum;
 	      else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
+	      else if (strcmp(name,"CalvingCalvingrateNoise")==0) return CalvingCalvingrateNoiseEnum;
 	      else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
-	      else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
+	      if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
+	      else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
 	      else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum;
 	      else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
@@ -751,9 +752,9 @@
 	      else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
 	      else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
-	      else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum;
+	      if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
+	      else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum;
 	      else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum;
 	      else if (strcmp(name,"HydrologyTwsSpc")==0) return HydrologyTwsSpcEnum;
@@ -874,9 +875,9 @@
 	      else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
 	      else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
-	      else if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum;
+	      if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum;
+	      else if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum;
 	      else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum;
 	      else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum;
@@ -918,4 +919,5 @@
 	      else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum;
 	      else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
+	      else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum;
 	      else if (strcmp(name,"SmbBasinsId")==0) return SmbBasinsIdEnum;
 	      else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
@@ -996,10 +998,10 @@
 	      else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
 	      else if (strcmp(name,"SmbT")==0) return SmbTEnum;
-	      else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
-	      else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
+	      if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
+	      else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
+	      else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
 	      else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
 	      else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
@@ -1023,4 +1025,5 @@
 	      else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum;
 	      else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
+	      else if (strcmp(name,"StochasticForcingDefaultId")==0) return StochasticForcingDefaultIdEnum;
 	      else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
 	      else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
@@ -1058,4 +1061,5 @@
 	      else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
 	      else if (strcmp(name,"TemperatureSEMIC")==0) return TemperatureSEMICEnum;
+	      else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum;
 	      else if (strcmp(name,"ThermalforcingValuesAutoregression")==0) return ThermalforcingValuesAutoregressionEnum;
 	      else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum;
@@ -1117,12 +1121,12 @@
 	      else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
 	      else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
-	      else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
+         else stage=10;
+   }
+   if(stage==10){
+	      if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
 	      else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
 	      else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
 	      else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum;
-         else stage=10;
-   }
-   if(stage==10){
-	      if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
+	      else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
 	      else if (strcmp(name,"Outputdefinition30")==0) return Outputdefinition30Enum;
 	      else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
@@ -1240,12 +1244,12 @@
 	      else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
 	      else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
-	      else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
+         else stage=11;
+   }
+   if(stage==11){
+	      if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
 	      else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
 	      else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum;
 	      else if (strcmp(name,"Channel")==0) return ChannelEnum;
-         else stage=11;
-   }
-   if(stage==11){
-	      if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
+	      else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
 	      else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum;
 	      else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum;
@@ -1363,12 +1367,12 @@
 	      else if (strcmp(name,"IntParam")==0) return IntParamEnum;
 	      else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
-	      else if (strcmp(name,"Inputs")==0) return InputsEnum;
+         else stage=12;
+   }
+   if(stage==12){
+	      if (strcmp(name,"Inputs")==0) return InputsEnum;
 	      else if (strcmp(name,"Internal")==0) return InternalEnum;
 	      else if (strcmp(name,"Intersect")==0) return IntersectEnum;
 	      else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
-         else stage=12;
-   }
-   if(stage==12){
-	      if (strcmp(name,"J")==0) return JEnum;
+	      else if (strcmp(name,"J")==0) return JEnum;
 	      else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum;
 	      else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum;
@@ -1486,12 +1490,12 @@
 	      else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
 	      else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
-	      else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
+         else stage=13;
+   }
+   if(stage==13){
+	      if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
 	      else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
 	      else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
 	      else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
-         else stage=13;
-   }
-   if(stage==13){
-	      if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
+	      else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
 	      else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
 	      else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;
Index: issm/trunk-jpl/src/m/classes/stochasticforcing.m
===================================================================
--- issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 26603)
+++ issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 26604)
@@ -8,5 +8,6 @@
 		isstochasticforcing  = 0;
 		fields               = NaN;
-		dimensions           = NaN;
+		defaultdimension     = 0;
+      default_id           = NaN;
 		covariance           = NaN;
 		randomflag           = 1;
@@ -34,25 +35,59 @@
 
 			num_fields = numel(self.fields);
-			size_tot   = sum(self.dimensions);
+			
+			%Check that covariance matrix is positive definite
+			try
+				chol(self.covariance);
+			catch
+				error('md.stochasticforcing.covariance is not positive definite');
+			end
+
+			%Check that all fields agree with the corresponding md class and if any field needs the default params   
+         checkdefaults = false; %need to check defaults only if one of the field does not have its own dimensionality
+			for field=self.fields
+            %Checking agreement of classes
+            if(contains(field,'SMB'))
+               if~(isequal(class(md.smb),char(field)))
+                  error('md.smb does not agree with stochasticforcing field %s', char(field));
+               end
+            end
+            if(contains(field,'frontalforcings'))
+               if~(isequal(class(md.frontalforcings),char(field)))
+                  error('md.frontalforcings does not agree with stochasticforcing field %s', char(field));
+               end
+            end
+            %Checking for specific dimensions
+            if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression'))
+               checkdefaults = true; %field with non-specific dimensionality
+            end
+         end
+
+			%Retrieve sum of all the field dimensionalities
+			size_tot   = self.defaultdimension*num_fields;
+			indSMBar = -1; %about to check for index of SMBautoregression
+			indTFar  = -1; %about to check for index of FrontalForcingsRignotAutoregression
+			if any(contains(self.fields,'SMBautoregression'))
+            size_tot = size_tot-self.defaultdimension+md.smb.num_basins;
+				indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with TFar timestep (08Nov2021)
+         end
+         if any(contains(self.fields,'FrontalForcingsRignotAutoregression'))
+            size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins;
+				indTFar  = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with SMBar timestep (08Nov2021)
+         end
+
+			if(indSMBar~=-1 && indTFar~=-1) %both autoregressive models are used: check autoregressive time step consistency
+				if((md.smb.ar_timestep~=md.frontalforcings.ar_timestep) && any(self.covariance(1+sum(self.dimensions(1:indSMBar-1)):sum(self.dimensions(1:indSMBar)),1+sum(self.dimensions(1:indTFar-1)):sum(self.dimensions(1:indTFar)))))
+					error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance');
+				end
+			end
 
 			md = checkfield(md,'fieldname','stochasticforcing.isstochasticforcing','values',[0 1]);
 			md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings()); %VV check here 'cell' (19Oct2021)
-			md = checkfield(md,'fieldname','stochasticforcing.dimensions','NaN',1,'Inf',1,'>',0,'size',[1,num_fields]); %specific dimension for each field
 			md = checkfield(md,'fieldname','stochasticforcing.covariance','NaN',1,'Inf',1,'size',[size_tot,size_tot]); %global covariance matrix
 			md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]);
-
-			%Check that all fields agree with the corresponding md class
-			for field=self.fields
-				if(contains(field,'SMB'))
-					if~(isequal(class(md.smb),char(field)))
-						error('md.smb does not agree with stochasticforcing field %s', char(field));
-					end
-				end
-				if(contains(field,'frontalforcings'))
-					if~(isequal(class(md.frontalforcings),char(field)))
-						error('md.frontalforcings does not agree with stochasticforcing field %s', char(field));
-					end
-				end
-			end
+			if(checkdefaults) %need to check the defaults
+            md = checkfield(md,'fieldname','stochasticforcing.defaultdimension','numel',1,'NaN',1,'Inf',1,'>',0);
+            md = checkfield(md,'fieldname','stochasticforcing.default_id','Inf',1,'>=',0,'<=',self.defaultdimension,'size',[md.mesh.numberofelements,1]);
+         end
 		end % }}}
 		function disp(self) % {{{
@@ -60,5 +95,6 @@
 			fielddisplay(self,'isstochasticforcing','is stochasticity activated?');
 			fielddisplay(self,'fields','fields with stochasticity applied, ex: {''SMBautoregression''}, or {''FrontalForcingsRignotAutoregression''}');
-			fielddisplay(self,'dimensions','dimensionality of each field');
+			fielddisplay(self,'defaultdimension','dimensionality of the noise terms (does not apply to fields with their specific dimension)');
+         fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)');
 			fielddisplay(self,'covariance','covariance matrix for within- and between-fields covariance (units must be squared field units)');
 			fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)');
@@ -77,10 +113,25 @@
 				return
 			else
+
+				%Retrieve dimensionality of each field
+				dimensions = self.defaultdimension*ones(1,num_fields);
+				ind = 1;
+				for field=self.fields
+					%Checking for specific dimensions
+					if(strcmp(field,'SMBautoregression'))
+						dimensions(ind) = md.smb.num_basins;
+					end
+					if(strcmp(field,'FrontalForcingsRignotAutoregression'))
+						dimensions(ind) = md.frontalforcings.num_basins;
+					end
+					ind = ind+1;
+				end
+
 				%Scaling covariance matrix (scale column-by-column and row-by-row)
-				scaledfields = {'SMBautoregression'}; %list of fields that need scaling *1/yts
+				scaledfields = {'DefaultCalving','SMBautoregression'}; %list of fields that need scaling *1/yts
 				tempcovariance = self.covariance; %copy of covariance to avoid writing back in member variable
 				for i=1:num_fields
 					if any(strcmp(scaledfields,self.fields(i)))
-						inds = [1+sum(self.dimensions(1:i-1)):1:sum(self.dimensions(1:i))];
+						inds = [1+sum(dimensions(1:i-1)):1:sum(dimensions(1:i))];
 						for row=inds %scale rows corresponding to scaled field
 							tempcovariance(row,:) = 1./yts*tempcovariance(row,:);
@@ -93,5 +144,7 @@
 				WriteData(fid,prefix,'data',num_fields,'name','md.stochasticforcing.num_fields','format','Integer');
 				WriteData(fid,prefix,'object',self,'fieldname','fields','format','StringArray');
-				WriteData(fid,prefix,'object',self,'fieldname','dimensions','format','IntMat');
+				WriteData(fid,prefix,'data',dimensions,'name','md.stochasticforcing.dimensions','format','IntMat');
+            WriteData(fid,prefix,'object',self,'fieldname','default_id','data',self.default_id-1,'format','IntMat','mattype',2); %0-indexed
+				WriteData(fid,prefix,'object',self,'fieldname','defaultdimension','format','Integer');
 				WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat');
 				WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean');
@@ -105,6 +158,7 @@
 
    list = {...
-      'SMBautoregression',...
-      'FrontalForcingsRignotAutoregression'
+      'DefaultCalving',...
+		'FrontalForcingsRignotAutoregression',...
+      'SMBautoregression'
       };
 end % }}}
