Changeset 26608

Timestamp:

11/11/21 09:56:50 (3 years ago)

Author:

Mathieu Morlighem

Message:

BUG: reverting Vincent's commit

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");
-        IoModelToDynamicConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,finiteelement);
-        //IoModelToConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,finiteelement);
+
+        /*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);
 }
 /*}}}*/
@@ -54 +83 @@
 
         /*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,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;
+void       Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type){/*{{{*/
+        
+        const int numvertices = this->GetNumberOfVertices();
+   int         basinid,M,N;
+   IssmDouble  beta0_basin,beta1_basin,noise_basin;
    IssmDouble* phi_basin  = xNew<IssmDouble>(arorder);
    IssmDouble* varlist     = xNew<IssmDouble>(numvertices);
    IssmDouble* valuesautoregression = NULL;
-   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*/
+
+   /*Get Basin ID and Basin coefficients*/
+   if(enum_type==SMBautoregressionEnum)                   this->GetInputValue(&basinid,SmbBasinsIdEnum);
+   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum);
    for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii];
    beta0_basin   = beta0[basinid];
    beta1_basin   = beta1[basinid];
+   noise_basin   = noiseterms[basinid];
+   if(enum_type==SMBautoregressionEnum)                   this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
+   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->GetArray(ThermalforcingValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
 
    /*If not AR model timestep: take the old values of variable*/
@@ -184 +131 @@
 
          /*Stochastic variable value*/
-         varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noiseterm;
+         varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin;
       }
-      /*Update autoregression values*/
+      /*Update autoregression TF values*/
       IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
       /*Assign newest values and shift older values*/
       for(int i=0;i<numvertices;i++) temparray[i] = varlist[i];
       for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = valuesautoregression[i];
-      this->inputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder);
+      if(enum_type==SMBautoregressionEnum)                   this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder);
+      if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder);
       xDelete<IssmDouble>(temparray);
    }
 
    /*Add input to element*/
-   this->AddInput(outenum_type,varlist,P1Enum);
+   if(enum_type==SMBautoregressionEnum)                   this->AddInput(SmbMassBalanceEnum,varlist,P1Enum);
+   if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum);
 
    /*Cleanup*/
    xDelete<IssmDouble>(phi_basin);
    xDelete<IssmDouble>(varlist);
-   xDelete<IssmDouble>(valuesautoregression); 
+   xDelete<IssmDouble>(valuesautoregression);
 }/*}}}*/
 void       Element::ComputeLambdaS(){/*{{{*/
@@ -3829 +3778 @@
         IssmDouble teValue=1.0;
         IssmDouble aValue=0.0;
+        IssmDouble dulwrfValue=0.0;
         IssmDouble szaValue=0.0;
         IssmDouble cotValue=0.0;
@@ -3835 +3785 @@
         IssmDouble dt,time,smb_dt;
         int        aIdx=0;
+        int        eIdx=0;
         int        denIdx=0;
         int        dsnowIdx=0;
@@ -3864 +3815 @@
         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:{{{ */
@@ -3916 +3869 @@
         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);
@@ -3932 +3886 @@
         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: {{{*/
@@ -4124 +4080 @@
         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);
@@ -4139 +4096 @@
         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]
@@ -4164 +4122 @@
         }
         /*Thermal profile computation:*/
-        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);
+        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);
 
         /*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,bool isfieldstochastic,int enum_type);
+      void               Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type);
                 void               ComputeLambdaS(void);
                 void               ComputeNewDamage();

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

@@ -80 +80 @@
    /*Load parameters*/
         bool isstochastic;
-   bool istfstochastic = false;
-        int M,N,Nphi,arorder,numbasins,my_rank;
+   int M,N,Nphi,arorder,numbasins,my_rank;
    femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
    femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
@@ -88 +87 @@
    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){
@@ -101 +102 @@
                 femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
                 for(int i=0;i<numstochasticfields;i++){
-                        if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum) istfstochastic = true;
+                        if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum){
+                                femmodel->parameters->FindParam(&noiseterms,&M,ThermalforcingAutoregressionNoiseEnum);  _assert_(M==numbasins);
+                        }
                 }
-                xDelete<int>(stochasticfields);
         }
    /*Time elapsed with respect to AR model initial time*/
@@ -111 +113 @@
    for(Object* &object:femmodel->elements->objects){
       Element* element = xDynamicCast<Element*>(object);
-      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum);
+      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,FrontalForcingsRignotAutoregressionEnum);
    }
 
@@ -118 +120 @@
    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,stochastic_dim;
+                int num_fields;
                 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");
@@ -516 +515 @@
                 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 dimenum_type,noiseenum_type;
+      int enum_type;
       IssmDouble* noisefield = xNew<IssmDouble>(dimensions[j]);
       for(int k=0;k<dimensions[j];k++){
          noisefield[k]=noiseterms[i+k];
       }
-     
-                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];
+      
+      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];
       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);
@@ -211 +210 @@
         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){
@@ -218 +219 @@
       femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
       for(int i=0;i<numstochasticfields;i++){
-         if(stochasticfields[i]==SMBautoregressionEnum) issmbstochastic = true;
+         if(stochasticfields[i]==SMBautoregressionEnum){
+                                femmodel->parameters->FindParam(&noiseterms,&M,SmbAutoregressionNoiseEnum);  _assert_(M==numbasins);
+                        }
                 }
                 xDelete<int>(stochasticfields);
@@ -228 +231 @@
         for(Object* &object:femmodel->elements->objects){
                 Element* element = xDynamicCast<Element*>(object);
-                element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum);
+                element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,SMBautoregressionEnum);
         }
 
@@ -235 +238 @@
         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
@@ -429 +428 @@
 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
@@ -598 +599 @@
 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
@@ -896 +895 @@
 syn keyword cConstant SmbAdiffiniEnum
 syn keyword cConstant SmbAiniEnum
-syn keyword cConstant SmbAutoregressionNoiseEnum
 syn keyword cConstant SmbBasinsIdEnum
 syn keyword cConstant SmbBMaxEnum
@@ -999 +997 @@
 syn keyword cConstant SolidearthExternalDisplacementUpRateEnum
 syn keyword cConstant SolidearthExternalGeoidRateEnum
-syn keyword cConstant StochasticForcingDefaultIdEnum
 syn keyword cConstant StrainRateeffectiveEnum
 syn keyword cConstant StrainRateparallelEnum
@@ -1035 +1032 @@
 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,
@@ -423 +422 @@
         SmbAdThreshEnum,
         SmbAutoregressionInitialTimeEnum,
+   SmbAutoregressionNoiseEnum,
         SmbAutoregressionTimestepEnum,
    SmbAutoregressiveOrderEnum,
@@ -498 +498 @@
         StressbalanceRiftPenaltyThresholdEnum,
         StressbalanceShelfDampeningEnum,
+   ThermalforcingAutoregressionNoiseEnum,
         ThermalIsdrainicecolumnEnum,
         ThermalIsdynamicbasalspcEnum,
@@ -594 +595 @@
         BasalStressEnum,
         BaseEnum,
-        BaselineCalvingCalvingrateEnum,
         BaseOldEnum,
         BaseSlopeXEnum,
@@ -609 +609 @@
         BottomPressureOldEnum,
         CalvingCalvingrateEnum,
-        CalvingCalvingrateNoiseEnum,
         CalvingHabFractionEnum,
         CalvingMeltingrateEnum,
@@ -892 +891 @@
         SmbAdiffiniEnum,
         SmbAiniEnum,
-   SmbAutoregressionNoiseEnum,
         SmbBasinsIdEnum,
         SmbBMaxEnum,
@@ -996 +994 @@
         SolidearthExternalDisplacementUpRateEnum,
         SolidearthExternalGeoidRateEnum,
-        StochasticForcingDefaultIdEnum,
         StrainRateeffectiveEnum,
         StrainRateparallelEnum,
@@ -1032 +1029 @@
         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";
@@ -431 +430 @@
                 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";
@@ -600 +601 @@
                 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";
@@ -898 +897 @@
                 case SmbAdiffiniEnum : return "SmbAdiffini";
                 case SmbAiniEnum : return "SmbAini";
-                case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
                 case SmbBasinsIdEnum : return "SmbBasinsId";
                 case SmbBMaxEnum : return "SmbBMax";
@@ -1001 +999 @@
                 case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate";
                 case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate";
-                case StochasticForcingDefaultIdEnum : return "StochasticForcingDefaultId";
                 case StrainRateeffectiveEnum : return "StrainRateeffective";
                 case StrainRateparallelEnum : return "StrainRateparallel";
@@ -1037 +1034 @@
                 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;
@@ -440 +439 @@
               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;
@@ -612 +613 @@
               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,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
-              else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
+              if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
               else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum;
               else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
@@ -752 +751 @@
               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,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
-              else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum;
+              if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum;
               else if (strcmp(name,"HydrologyTws")==0) return HydrologyTwsEnum;
               else if (strcmp(name,"HydrologyTwsSpc")==0) return HydrologyTwsSpcEnum;
@@ -875 +874 @@
               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,"SealevelchangeG")==0) return SealevelchangeGEnum;
-              else if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum;
+              if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum;
               else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum;
               else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum;
@@ -919 +918 @@
               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;
@@ -998 +996 @@
               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,"SmbTa")==0) return SmbTaEnum;
-              else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
-              else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
+              if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
               else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
               else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
@@ -1025 +1023 @@
               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;
@@ -1061 +1058 @@
               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;
@@ -1121 +1117 @@
               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 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,"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 if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
+              if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
               else if (strcmp(name,"Outputdefinition30")==0) return Outputdefinition30Enum;
               else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
@@ -1244 +1240 @@
               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 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,"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 if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
+              if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
               else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum;
               else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum;
@@ -1367 +1363 @@
               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 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,"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 if (strcmp(name,"J")==0) return JEnum;
+              if (strcmp(name,"J")==0) return JEnum;
               else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum;
               else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum;
@@ -1490 +1486 @@
               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 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,"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 if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
+              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;
-                defaultdimension     = 0;
-      default_id           = NaN;
+                dimensions           = NaN;
                 covariance           = NaN;
                 randomflag           = 1;
@@ -35 +34 @@
 
                         num_fields = numel(self.fields);
-                        
-                        %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
+                        size_tot   = sum(self.dimensions);
 
                         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]);
-                        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
+
+                        %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
                 end % }}}
                 function disp(self) % {{{
@@ -95 +60 @@
                         fielddisplay(self,'isstochasticforcing','is stochasticity activated?');
                         fielddisplay(self,'fields','fields with stochasticity applied, ex: {''SMBautoregression''}, or {''FrontalForcingsRignotAutoregression''}');
-                        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,'dimensions','dimensionality of each field');
                         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)');
@@ -113 +77 @@
                                 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 = {'DefaultCalving','SMBautoregression'}; %list of fields that need scaling *1/yts
+                                scaledfields = {'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(dimensions(1:i-1)):1:sum(dimensions(1:i))];
+                                                inds = [1+sum(self.dimensions(1:i-1)):1:sum(self.dimensions(1:i))];
                                                 for row=inds %scale rows corresponding to scaled field
                                                         tempcovariance(row,:) = 1./yts*tempcovariance(row,:);
@@ -144 +93 @@
                                 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,'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,'object',self,'fieldname','dimensions','format','IntMat');
                                 WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat');
                                 WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean');
@@ -158 +105 @@
 
    list = {...
-      'DefaultCalving',...
-                'FrontalForcingsRignotAutoregression',...
-      'SMBautoregression'
+      'SMBautoregression',...
+      'FrontalForcingsRignotAutoregression'
       };
 end % }}}