Changeset 26604 for issm/trunk-jpl/src

Timestamp:

11/11/21 09:33:47 (3 years ago)

Author:

vverjans

Message:

CHG: making stochasticforcing.m more generic, extending stochasticforcing capability to DefaultCalving

Location:

issm/trunk-jpl/src

Files:

• c/analyses/LevelsetAnalysis.cpp (modified) (2 diffs)
• c/classes/Elements/Element.cpp (modified) (10 diffs)
• c/classes/Elements/Element.h (modified) (1 diff)
• c/modules/FrontalForcingsx/FrontalForcingsx.cpp (modified) (6 diffs)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (2 diffs)
• c/modules/StochasticForcingx/StochasticForcingx.cpp (modified) (2 diffs)
• c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (5 diffs)
• c/shared/Enum/Enum.vim (modified) (8 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (8 diffs)
• c/shared/Enum/EnumToStringx.cpp (modified) (8 diffs)
• c/shared/Enum/StringToEnumx.cpp (modified) (15 diffs)
• m/classes/stochasticforcing.m (modified) (6 diffs)

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

@@ -12 +12 @@
 
 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 +54 @@
 
         /*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:

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

@@ -99 +99 @@
    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 +184 @@
 
          /*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 +3829 @@
         IssmDouble teValue=1.0;
         IssmDouble aValue=0.0;
-        IssmDouble dulwrfValue=0.0;
         IssmDouble szaValue=0.0;
         IssmDouble cotValue=0.0;
@@ -3785 +3835 @@
         IssmDouble dt,time,smb_dt;
         int        aIdx=0;
-        int        eIdx=0;
         int        denIdx=0;
         int        dsnowIdx=0;
@@ -3815 +3864 @@
         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 +3916 @@
         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 +3932 @@
         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 +4124 @@
         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 +4139 @@
         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 +4164 @@
         }
         /*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.

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

@@ -69 +69 @@
                 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();

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

@@ -80 +80 @@
    /*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 +88 @@
    IssmDouble* beta1      = NULL;
    IssmDouble* phi        = NULL;
-   IssmDouble* noiseterms = xNewZeroInit<IssmDouble>(numbasins);
 
         femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
@@ -94 +94 @@
    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 +101 @@
                 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 +111 @@
    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 +118 @@
    xDelete<IssmDouble>(beta1);
    xDelete<IssmDouble>(phi);
-   xDelete<IssmDouble>(noiseterms);
 }/*}}}*/

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

@@ -502 +502 @@
         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 +516 @@
                 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));

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

@@ -36 +36 @@
       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 +46 @@
         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);
    }

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

@@ -197 +197 @@
         /*Load parameters*/
         bool isstochastic;
+        bool issmbstochastic = false;
         int M,N,Nphi,arorder,numbasins,my_rank;
         femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
@@ -210 +211 @@
         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 +218 @@
       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 +228 @@
         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 +235 @@
         xDelete<IssmDouble>(beta1);
         xDelete<IssmDouble>(phi);
-        xDelete<IssmDouble>(noiseterms);
 }/*}}}*/
 void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/

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

@@ -398 +398 @@
 syn keyword cConstant SolidearthSettingsOceanAreaScalingEnum
 syn keyword cConstant StochasticForcingCovarianceEnum
+syn keyword cConstant StochasticForcingDefaultDimensionEnum
 syn keyword cConstant StochasticForcingDimensionsEnum
 syn keyword cConstant StochasticForcingFieldsEnum
@@ -428 +429 @@
 syn keyword cConstant SmbAdThreshEnum
 syn keyword cConstant SmbAutoregressionInitialTimeEnum
-syn keyword cConstant SmbAutoregressionNoiseEnum
 syn keyword cConstant SmbAutoregressionTimestepEnum
 syn keyword cConstant SmbAutoregressiveOrderEnum
@@ -504 +504 @@
 syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum
 syn keyword cConstant StressbalanceShelfDampeningEnum
-syn keyword cConstant ThermalforcingAutoregressionNoiseEnum
 syn keyword cConstant ThermalIsdrainicecolumnEnum
 syn keyword cConstant ThermalIsdynamicbasalspcEnum
@@ -599 +598 @@
 syn keyword cConstant BasalStressEnum
 syn keyword cConstant BaseEnum
+syn keyword cConstant BaselineCalvingCalvingrateEnum
 syn keyword cConstant BaseOldEnum
 syn keyword cConstant BaseSlopeXEnum
@@ -613 +613 @@
 syn keyword cConstant BottomPressureOldEnum
 syn keyword cConstant CalvingCalvingrateEnum
+syn keyword cConstant CalvingCalvingrateNoiseEnum
 syn keyword cConstant CalvingHabFractionEnum
 syn keyword cConstant CalvingMeltingrateEnum
@@ -895 +896 @@
 syn keyword cConstant SmbAdiffiniEnum
 syn keyword cConstant SmbAiniEnum
+syn keyword cConstant SmbAutoregressionNoiseEnum
 syn keyword cConstant SmbBasinsIdEnum
 syn keyword cConstant SmbBMaxEnum
@@ -997 +999 @@
 syn keyword cConstant SolidearthExternalDisplacementUpRateEnum
 syn keyword cConstant SolidearthExternalGeoidRateEnum
+syn keyword cConstant StochasticForcingDefaultIdEnum
 syn keyword cConstant StrainRateeffectiveEnum
 syn keyword cConstant StrainRateparallelEnum
@@ -1032 +1035 @@
 syn keyword cConstant TemperaturePicardEnum
 syn keyword cConstant TemperatureSEMICEnum
+syn keyword cConstant ThermalforcingAutoregressionNoiseEnum
 syn keyword cConstant ThermalforcingValuesAutoregressionEnum
 syn keyword cConstant ThermalSpctemperatureEnum

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

@@ -392 +392 @@
         SolidearthSettingsOceanAreaScalingEnum,
         StochasticForcingCovarianceEnum,
+        StochasticForcingDefaultDimensionEnum,
         StochasticForcingDimensionsEnum,
         StochasticForcingFieldsEnum,
@@ -422 +423 @@
         SmbAdThreshEnum,
         SmbAutoregressionInitialTimeEnum,
-   SmbAutoregressionNoiseEnum,
         SmbAutoregressionTimestepEnum,
    SmbAutoregressiveOrderEnum,
@@ -498 +498 @@
         StressbalanceRiftPenaltyThresholdEnum,
         StressbalanceShelfDampeningEnum,
-   ThermalforcingAutoregressionNoiseEnum,
         ThermalIsdrainicecolumnEnum,
         ThermalIsdynamicbasalspcEnum,
@@ -595 +594 @@
         BasalStressEnum,
         BaseEnum,
+        BaselineCalvingCalvingrateEnum,
         BaseOldEnum,
         BaseSlopeXEnum,
@@ -609 +609 @@
         BottomPressureOldEnum,
         CalvingCalvingrateEnum,
+        CalvingCalvingrateNoiseEnum,
         CalvingHabFractionEnum,
         CalvingMeltingrateEnum,
@@ -891 +892 @@
         SmbAdiffiniEnum,
         SmbAiniEnum,
+   SmbAutoregressionNoiseEnum,
         SmbBasinsIdEnum,
         SmbBMaxEnum,
@@ -994 +996 @@
         SolidearthExternalDisplacementUpRateEnum,
         SolidearthExternalGeoidRateEnum,
+        StochasticForcingDefaultIdEnum,
         StrainRateeffectiveEnum,
         StrainRateparallelEnum,
@@ -1029 +1032 @@
         TemperaturePicardEnum,
         TemperatureSEMICEnum,
+   ThermalforcingAutoregressionNoiseEnum,
         ThermalforcingValuesAutoregressionEnum, 
         ThermalSpctemperatureEnum,

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

@@ -400 +400 @@
                 case SolidearthSettingsOceanAreaScalingEnum : return "SolidearthSettingsOceanAreaScaling";
                 case StochasticForcingCovarianceEnum : return "StochasticForcingCovariance";
+                case StochasticForcingDefaultDimensionEnum : return "StochasticForcingDefaultDimension";
                 case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions";
                 case StochasticForcingFieldsEnum : return "StochasticForcingFields";
@@ -430 +431 @@
                 case SmbAdThreshEnum : return "SmbAdThresh";
                 case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime";
-                case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
                 case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep";
                 case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder";
@@ -506 +506 @@
                 case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold";
                 case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening";
-                case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise";
                 case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn";
                 case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc";
@@ -601 +600 @@
                 case BasalStressEnum : return "BasalStress";
                 case BaseEnum : return "Base";
+                case BaselineCalvingCalvingrateEnum : return "BaselineCalvingCalvingrate";
                 case BaseOldEnum : return "BaseOld";
                 case BaseSlopeXEnum : return "BaseSlopeX";
@@ -615 +615 @@
                 case BottomPressureOldEnum : return "BottomPressureOld";
                 case CalvingCalvingrateEnum : return "CalvingCalvingrate";
+                case CalvingCalvingrateNoiseEnum : return "CalvingCalvingrateNoise";
                 case CalvingHabFractionEnum : return "CalvingHabFraction";
                 case CalvingMeltingrateEnum : return "CalvingMeltingrate";
@@ -897 +898 @@
                 case SmbAdiffiniEnum : return "SmbAdiffini";
                 case SmbAiniEnum : return "SmbAini";
+                case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
                 case SmbBasinsIdEnum : return "SmbBasinsId";
                 case SmbBMaxEnum : return "SmbBMax";
@@ -999 +1001 @@
                 case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate";
                 case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate";
+                case StochasticForcingDefaultIdEnum : return "StochasticForcingDefaultId";
                 case StrainRateeffectiveEnum : return "StrainRateeffective";
                 case StrainRateparallelEnum : return "StrainRateparallel";
@@ -1034 +1037 @@
                 case TemperaturePicardEnum : return "TemperaturePicard";
                 case TemperatureSEMICEnum : return "TemperatureSEMIC";
+                case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise";
                 case ThermalforcingValuesAutoregressionEnum : return "ThermalforcingValuesAutoregression";
                 case ThermalSpctemperatureEnum : return "ThermalSpctemperature";

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

@@ -409 +409 @@
               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 +440 @@
               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 +518 @@
               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 +612 @@
               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 +627 @@
               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 +752 @@
               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 +875 @@
               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 +919 @@
               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 +998 @@
               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 +1025 @@
               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 +1061 @@
               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 +1121 @@
               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 +1244 @@
               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 +1367 @@
               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 +1490 @@
               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;

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

@@ -8 +8 @@
                 isstochasticforcing  = 0;
                 fields               = NaN;
-                dimensions           = NaN;
+                defaultdimension     = 0;
+      default_id           = NaN;
                 covariance           = NaN;
                 randomflag           = 1;
@@ -34 +35 @@
 
                         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 +95 @@
                         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 +113 @@
                                 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 +144 @@
                                 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 +158 @@
 
    list = {...
-      'SMBautoregression',...
-      'FrontalForcingsRignotAutoregression'
+      'DefaultCalving',...
+                'FrontalForcingsRignotAutoregression',...
+      'SMBautoregression'
       };
 end % }}}