source: issm/oecreview/Archive/25834-26739/ISSM-26614-26615.diff

Last change on this file was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 65.3 KB

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

     
    8282        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
    8383
    8484        /*Get moving front parameters*/
    85         int  calvinglaw;
    86         iomodel->FindConstant(&calvinglaw,"md.calving.law");
    87         switch(calvinglaw){
    88                 case DefaultCalvingEnum:
    89                         iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",CalvingCalvingrateEnum);
    90                         break;
     85        bool isstochastic;
     86   int  calvinglaw;
     87   iomodel->FindConstant(&calvinglaw,"md.calving.law");
     88   iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing");
     89   switch(calvinglaw){
     90      case DefaultCalvingEnum:
     91         iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",CalvingCalvingrateEnum);
     92         if(isstochastic){
     93            iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
     94            iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",BaselineCalvingCalvingrateEnum);
     95         }
     96         break;
    9197                case CalvingLevermannEnum:
    9298                        iomodel->FetchDataToInput(inputs,elements,"md.calving.coeff",CalvinglevermannCoeffEnum);
    9399                        break;

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

     
    9898   xDelete<IssmDouble>(varspin);
    9999   xDelete<IssmDouble>(phi_basin);     
    100100}/*}}}*/
    101 void       Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type){/*{{{*/
    102        
    103         const int numvertices = this->GetNumberOfVertices();
    104    int         basinid,M,N;
    105    IssmDouble  beta0_basin,beta1_basin,noise_basin;
     101void       Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type){/*{{{*/
     102
     103   const int numvertices = this->GetNumberOfVertices();
     104   int         basinid,M,N,arenum_type,basinenum_type,noiseenum_type,outenum_type;
     105   IssmDouble  beta0_basin,beta1_basin,noiseterm;
    106106   IssmDouble* phi_basin  = xNew<IssmDouble>(arorder);
    107107   IssmDouble* varlist     = xNew<IssmDouble>(numvertices);
    108108   IssmDouble* valuesautoregression = NULL;
     109   Input*      noiseterm_input      = NULL;
    109110
    110    /*Get Basin ID and Basin coefficients*/
    111    if(enum_type==SMBautoregressionEnum)                   this->GetInputValue(&basinid,SmbBasinsIdEnum);
    112    if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum);
     111   /*Get field-specific enums*/
     112   switch(enum_type){
     113      case(SMBautoregressionEnum):
     114         arenum_type    = SmbValuesAutoregressionEnum;
     115         basinenum_type = SmbBasinsIdEnum;
     116         noiseenum_type = SmbAutoregressionNoiseEnum;
     117         outenum_type   = SmbMassBalanceEnum;
     118         break;
     119      case(FrontalForcingsRignotAutoregressionEnum):
     120         arenum_type    = ThermalforcingValuesAutoregressionEnum;
     121         basinenum_type = FrontalForcingsBasinIdEnum;
     122         noiseenum_type = ThermalforcingAutoregressionNoiseEnum;
     123         outenum_type   = FrontalForcingsThermalForcingEnum;
     124         break;
     125   }
     126
     127   /*Get noise and autoregressive terms*/
     128   this->GetInputValue(&basinid,basinenum_type);
     129   if(isfieldstochastic){
     130      noiseterm_input = this->GetInput(noiseenum_type);
     131      Gauss* gauss = this->NewGauss();
     132      noiseterm_input->GetInputValue(&noiseterm,gauss);
     133      delete gauss;
     134   }
     135   else noiseterm = 0.0;
     136   this->inputs->GetArray(arenum_type,this->lid,&valuesautoregression,&M);
     137
     138   /*Get basin coefficients*/
    113139   for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii];
    114140   beta0_basin   = beta0[basinid];
    115141   beta1_basin   = beta1[basinid];
    116    noise_basin   = noiseterms[basinid];
    117    if(enum_type==SMBautoregressionEnum)                   this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
    118    if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->GetArray(ThermalforcingValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
    119142
    120    /*If not AR model timestep: take the old values of variable*/
     143        /*If not AR model timestep: take the old values of variable*/
    121144   if(isstepforar==false){
    122145      for(int i=0;i<numvertices;i++) varlist[i]=valuesautoregression[i];
    123146   }
    … …  
    130153         for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*valuesautoregression[v+ii*numvertices];
    131154
    132155         /*Stochastic variable value*/
    133          varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin;
     156         varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noiseterm;
    134157      }
    135       /*Update autoregression TF values*/
     158      /*Update autoregression values*/
    136159      IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
    137160      /*Assign newest values and shift older values*/
    138161      for(int i=0;i<numvertices;i++) temparray[i] = varlist[i];
    139162      for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = valuesautoregression[i];
    140       if(enum_type==SMBautoregressionEnum)                   this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder);
    141       if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder);
     163      this->inputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder);
    142164      xDelete<IssmDouble>(temparray);
    143165   }
    144166
    145167   /*Add input to element*/
    146    if(enum_type==SMBautoregressionEnum)                   this->AddInput(SmbMassBalanceEnum,varlist,P1Enum);
    147    if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum);
     168   this->AddInput(outenum_type,varlist,P1Enum);
    148169
    149170   /*Cleanup*/
    150171   xDelete<IssmDouble>(phi_basin);

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

     
    6868                /*bool               AnyActive(void);*/
    6969                bool               AnyFSet(void);
    7070                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);
    71       void               Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type);
     71      void               Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type);
    7272                void               ComputeLambdaS(void);
    7373                void               ComputeNewDamage();
    7474                void               ComputeStrainRate();

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

     
    7979
    8080   /*Load parameters*/
    8181        bool isstochastic;
    82    int M,N,Nphi,arorder,numbasins,my_rank;
     82   bool istfstochastic = false;
     83        int M,N,Nphi,arorder,numbasins,my_rank;
    8384   femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
    8485   femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
    8586   IssmDouble tinit_ar;
    … …  
    8687   IssmDouble* beta0      = NULL;
    8788   IssmDouble* beta1      = NULL;
    8889   IssmDouble* phi        = NULL;
    89    IssmDouble* noiseterms = xNewZeroInit<IssmDouble>(numbasins);
    9090
    9191        femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
    9292   femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum);    _assert_(M==numbasins);
    … …  
    9393   femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum);    _assert_(M==numbasins);
    9494   femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
    9595
    96         /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/
    9796        femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
    9897        if(isstochastic){
    9998                int  numstochasticfields;
    … …  
    101100                femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
    102101                femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
    103102                for(int i=0;i<numstochasticfields;i++){
    104                         if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum){
    105                                 femmodel->parameters->FindParam(&noiseterms,&M,ThermalforcingAutoregressionNoiseEnum);  _assert_(M==numbasins);
    106                         }
     103                        if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum) istfstochastic = true;
    107104                }
     105                xDelete<int>(stochasticfields);
    108106        }
    109107   /*Time elapsed with respect to AR model initial time*/
    110108   IssmDouble telapsed_ar = time-tinit_ar;
    … …  
    112110   /*Loop over each element to compute Thermal Forcing at vertices*/
    113111   for(Object* &object:femmodel->elements->objects){
    114112      Element* element = xDynamicCast<Element*>(object);
    115       element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,FrontalForcingsRignotAutoregressionEnum);
     113      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum);
    116114   }
    117115
    118116   /*Cleanup*/
    … …  
    119117   xDelete<IssmDouble>(beta0);
    120118   xDelete<IssmDouble>(beta1);
    121119   xDelete<IssmDouble>(phi);
    122    xDelete<IssmDouble>(noiseterms);
    123120}/*}}}*/

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

     
    499499        }
    500500
    501501        bool isstochasticforcing;
    502         parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
    503         if(isstochasticforcing){
    504                 int num_fields;
    505                 char** fields;
    506                 parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum));
    507                 iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields");
    508                 if(num_fields<1) _error_("no stochasticforcing fields found");
    509                 int* stochasticforcing_enums = xNew<int>(num_fields);
    510                 for(int i=0;i<num_fields;i++){
    511                         stochasticforcing_enums[i] = StringToEnumx(fields[i]);
    512                         xDelete<char>(fields[i]);
    513                 }
    514                 xDelete<char*>(fields);
    515                 parameters->AddObject(new IntVecParam(StochasticForcingFieldsEnum,stochasticforcing_enums,num_fields));
    516                 xDelete<int>(stochasticforcing_enums);
    517 
     502   parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
     503   if(isstochasticforcing){
     504      int num_fields,stochastic_dim;
     505      char** fields;
     506      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum));
     507      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.defaultdimension",StochasticForcingDefaultDimensionEnum));
     508      iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields");
     509      if(num_fields<1) _error_("no stochasticforcing fields found");
     510      int* stochasticforcing_enums = xNew<int>(num_fields);
     511      for(int i=0;i<num_fields;i++){
     512         stochasticforcing_enums[i] = StringToEnumx(fields[i]);
     513         xDelete<char>(fields[i]);
     514      }
     515      xDelete<char*>(fields);
     516      parameters->AddObject(new IntVecParam(StochasticForcingFieldsEnum,stochasticforcing_enums,num_fields));
     517      xDelete<int>(stochasticforcing_enums);
    518518      parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.randomflag",StochasticForcingRandomflagEnum));
    519                 iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.dimensions");
    520                 parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N));
    521                 xDelete<IssmDouble>(transparam);
    522                 iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.covariance");
    523                 parameters->AddObject(new DoubleMatParam(StochasticForcingCovarianceEnum,transparam,M,N));
    524                 xDelete<IssmDouble>(transparam);
    525         }
     519      iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.dimensions");
     520      parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N));
     521      xDelete<IssmDouble>(transparam);
     522      iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.covariance");
     523      parameters->AddObject(new DoubleMatParam(StochasticForcingCovarianceEnum,transparam,M,N));
     524      xDelete<IssmDouble>(transparam);
     525   }
    526526
    527527        /*Deal with mass flux segments: {{{*/
    528528        iomodel->FetchData(&qmu_mass_flux_present,"md.qmu.mass_flux_segments_present");

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

     
    3535                /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
    3636      if(randomflag) fixedseed=-1;
    3737      else fixedseed = reCast<int,IssmDouble>((time-starttime)/dt);
    38       /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
     38                /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
    3939      IssmDouble* temparray = NULL;
    4040      multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed);
    4141      for(int i=0;i<dimtot;i++) noiseterms[i]=temparray[i];
    … …  
    4545   
    4646        int i=0;
    4747   for(int j=0;j<numfields;j++){
    48       int enum_type;
     48      int dimenum_type,noiseenum_type;
    4949      IssmDouble* noisefield = xNew<IssmDouble>(dimensions[j]);
    5050      for(int k=0;k<dimensions[j];k++){
    5151         noisefield[k]=noiseterms[i+k];
    5252      }
    53      
    54       if(fields[j]==SMBautoregressionEnum)                        enum_type = SmbAutoregressionNoiseEnum;
    55                 else if(fields[j]==FrontalForcingsRignotAutoregressionEnum) enum_type = ThermalforcingAutoregressionNoiseEnum;
    56                 else _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.\n");
    57       femmodel->parameters->SetParam(noisefield,dimensions[j],enum_type);
    58       i=i+dimensions[j];
     53     
     54                int dimensionid;
     55
     56                /*Deal with the autoregressive models*/
     57                if(fields[j]==SMBautoregressionEnum || fields[j]==FrontalForcingsRignotAutoregressionEnum){
     58                        switch(fields[j]){
     59                                case SMBautoregressionEnum:
     60                                        dimenum_type   = SmbBasinsIdEnum;
     61                                        noiseenum_type = SmbAutoregressionNoiseEnum;
     62                                        break;
     63                                case FrontalForcingsRignotAutoregressionEnum:
     64                                        dimenum_type   = FrontalForcingsBasinIdEnum;
     65                                        noiseenum_type = ThermalforcingAutoregressionNoiseEnum;
     66                                        break;
     67                        }
     68                        for(Object* &object:femmodel->elements->objects){
     69            Element* element = xDynamicCast<Element*>(object);
     70            int numvertices  = element->GetNumberOfVertices();
     71            IssmDouble* noise_element = xNew<IssmDouble>(numvertices);
     72            element->GetInputValue(&dimensionid,dimenum_type);
     73            for(int i=0;i<numvertices;i++) noise_element[i] = noisefield[dimensionid];
     74            element->AddInput(noiseenum_type,noise_element,P0Enum);
     75            xDelete<IssmDouble>(noise_element);
     76                        }
     77                }
     78                else{
     79                        switch(fields[j]){
     80                                case SMBautoregressionEnum:
     81                                case FrontalForcingsRignotAutoregressionEnum:
     82                                        /*Already done above*/
     83                                        break;
     84                                case DefaultCalvingEnum:
     85                                        /*Delete CalvingCalvingrateEnum at previous time step (required if it is transient)*/
     86                                        femmodel->inputs->DeleteInput(CalvingCalvingrateEnum);
     87                                        for(Object* &object:femmodel->elements->objects){
     88                                                Element* element = xDynamicCast<Element*>(object);
     89                                                int numvertices  = element->GetNumberOfVertices();
     90                                                IssmDouble baselinecalvingrate;
     91                                                IssmDouble calvingrate_tot[numvertices];
     92                                                Input* baselinecalvingrate_input  = NULL;
     93                                                baselinecalvingrate_input = element->GetInput(BaselineCalvingCalvingrateEnum); _assert_(baselinecalvingrate_input);
     94                                                element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
     95                                                Gauss* gauss = element->NewGauss();
     96                                                for(int i=0;i<numvertices;i++){
     97                                                        gauss->GaussVertex(i);
     98                                                        baselinecalvingrate_input->GetInputValue(&baselinecalvingrate,gauss);
     99                                                        calvingrate_tot[i] = max(0.0,baselinecalvingrate+noisefield[dimensionid]);
     100                                                }
     101                                                element->AddInput(CalvingCalvingrateEnum,&calvingrate_tot[0],P1DGEnum);
     102                                                delete gauss;
     103                                        }
     104                                        break;
     105                                default:
     106                                        _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.");
     107                        }
     108                }
     109                i=i+dimensions[j];
    59110      xDelete<IssmDouble>(noisefield);
    60111   }
    61112

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

     
    176176}/*}}}*/
    177177void Smbautoregressionx(FemModel* femmodel){/*{{{*/
    178178
    179         /*Get time parameters*/
    180         IssmDouble time,dt,starttime,tstep_ar;
    181         femmodel->parameters->FindParam(&time,TimeEnum);
    182         femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    183         femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
     179   /*Get time parameters*/
     180   IssmDouble time,dt,starttime,tstep_ar;
     181   femmodel->parameters->FindParam(&time,TimeEnum);
     182   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     183   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    184184   femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
    185185
    186         /*Initialize module at first time step*/
    187         if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
    188         /*Determine if this is a time step for the AR model*/
    189         bool isstepforar = false;
     186   /*Initialize module at first time step*/
     187   if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
     188   /*Determine if this is a time step for the AR model*/
     189   bool isstepforar = false;
    190190
    191         #ifndef _HAVE_AD_
     191   #ifndef _HAVE_AD_
    192192   if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt) isstepforar = true;
    193         #else
    194         _error_("not implemented yet");
    195         #endif
     193   #else
     194   _error_("not implemented yet");
     195   #endif
    196196
    197         /*Load parameters*/
    198         bool isstochastic;
    199         int M,N,Nphi,arorder,numbasins,my_rank;
    200         femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
    201         femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
    202         IssmDouble tinit_ar;
    203         IssmDouble* beta0      = NULL;
    204         IssmDouble* beta1      = NULL;
    205         IssmDouble* phi        = NULL;
     197   /*Load parameters*/
     198   bool isstochastic;
     199   bool issmbstochastic = false;
     200   int M,N,Nphi,arorder,numbasins,my_rank;
     201   femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
     202   femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
     203   IssmDouble tinit_ar;
     204   IssmDouble* beta0      = NULL;
     205   IssmDouble* beta1      = NULL;
     206   IssmDouble* phi        = NULL;
    206207
    207         femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
    208         femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
    209         femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
    210         femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
     208   femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
     209   femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
     210   femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
     211   femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
    211212
    212         /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/
    213    IssmDouble* noiseterms = xNewZeroInit<IssmDouble>(numbasins);
    214         femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
     213   femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
    215214   if(isstochastic){
    216                 int  numstochasticfields;
     215      int  numstochasticfields;
    217216      int* stochasticfields;
    218217      femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
    219218      femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
    220219      for(int i=0;i<numstochasticfields;i++){
    221          if(stochasticfields[i]==SMBautoregressionEnum){
    222                                 femmodel->parameters->FindParam(&noiseterms,&M,SmbAutoregressionNoiseEnum);  _assert_(M==numbasins);
    223                         }
    224                 }
    225                 xDelete<int>(stochasticfields);
    226         }
    227         /*Time elapsed with respect to AR model initial time*/
    228         IssmDouble telapsed_ar = time-tinit_ar;
     220         if(stochasticfields[i]==SMBautoregressionEnum) issmbstochastic = true;
     221      }
     222      xDelete<int>(stochasticfields);
     223   }
     224   /*Time elapsed with respect to AR model initial time*/
     225   IssmDouble telapsed_ar = time-tinit_ar;
    229226
    230         /*Loop over each element to compute SMB at vertices*/
    231         for(Object* &object:femmodel->elements->objects){
    232                 Element* element = xDynamicCast<Element*>(object);
    233                 element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,SMBautoregressionEnum);
    234         }
     227   /*Loop over each element to compute SMB at vertices*/
     228   for(Object* &object:femmodel->elements->objects){
     229      Element* element = xDynamicCast<Element*>(object);
     230      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum);
     231   }
    235232
    236         /*Cleanup*/
    237         xDelete<IssmDouble>(beta0);
    238         xDelete<IssmDouble>(beta1);
    239         xDelete<IssmDouble>(phi);
    240         xDelete<IssmDouble>(noiseterms);
     233   /*Cleanup*/
     234   xDelete<IssmDouble>(beta0);
     235   xDelete<IssmDouble>(beta1);
     236   xDelete<IssmDouble>(phi);
    241237}/*}}}*/
    242238void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
    243239

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

     
    397397syn keyword cConstant SolidearthSettingsGrdOceanEnum
    398398syn keyword cConstant SolidearthSettingsOceanAreaScalingEnum
    399399syn keyword cConstant StochasticForcingCovarianceEnum
     400syn keyword cConstant StochasticForcingDefaultDimensionEnum
    400401syn keyword cConstant StochasticForcingDimensionsEnum
    401402syn keyword cConstant StochasticForcingFieldsEnum
    402403syn keyword cConstant StochasticForcingIsStochasticForcingEnum
    … …  
    427428syn keyword cConstant SmbAccurefEnum
    428429syn keyword cConstant SmbAdThreshEnum
    429430syn keyword cConstant SmbAutoregressionInitialTimeEnum
    430 syn keyword cConstant SmbAutoregressionNoiseEnum
    431431syn keyword cConstant SmbAutoregressionTimestepEnum
    432432syn keyword cConstant SmbAutoregressiveOrderEnum
    433433syn keyword cConstant SmbAveragingEnum
    … …  
    503503syn keyword cConstant StressbalanceRestolEnum
    504504syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum
    505505syn keyword cConstant StressbalanceShelfDampeningEnum
    506 syn keyword cConstant ThermalforcingAutoregressionNoiseEnum
    507506syn keyword cConstant ThermalIsdrainicecolumnEnum
    508507syn keyword cConstant ThermalIsdynamicbasalspcEnum
    509508syn keyword cConstant ThermalIsenthalpyEnum
    … …  
    601600syn keyword cConstant BaseOldEnum
    602601syn keyword cConstant BaseSlopeXEnum
    603602syn keyword cConstant BaseSlopeYEnum
     603syn keyword cConstant BaselineCalvingCalvingrateEnum
    604604syn keyword cConstant BedEnum
    605605syn keyword cConstant BedGRDEnum
    606606syn keyword cConstant BedEastEnum
    … …  
    894894syn keyword cConstant SmbAccumulationEnum
    895895syn keyword cConstant SmbAdiffiniEnum
    896896syn keyword cConstant SmbAiniEnum
     897syn keyword cConstant SmbAutoregressionNoiseEnum
    897898syn keyword cConstant SmbBasinsIdEnum
    898899syn keyword cConstant SmbBMaxEnum
    899900syn keyword cConstant SmbBMinEnum
    … …  
    996997syn keyword cConstant SolidearthExternalDisplacementNorthRateEnum
    997998syn keyword cConstant SolidearthExternalDisplacementUpRateEnum
    998999syn keyword cConstant SolidearthExternalGeoidRateEnum
     1000syn keyword cConstant StochasticForcingDefaultIdEnum
    9991001syn keyword cConstant StrainRateeffectiveEnum
    10001002syn keyword cConstant StrainRateparallelEnum
    10011003syn keyword cConstant StrainRateperpendicularEnum
    … …  
    10311033syn keyword cConstant TemperaturePDDEnum
    10321034syn keyword cConstant TemperaturePicardEnum
    10331035syn keyword cConstant TemperatureSEMICEnum
     1036syn keyword cConstant ThermalforcingAutoregressionNoiseEnum
    10341037syn keyword cConstant ThermalforcingValuesAutoregressionEnum
    10351038syn keyword cConstant ThermalSpctemperatureEnum
    10361039syn keyword cConstant ThicknessAbsGradientEnum

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

     
    391391        SolidearthSettingsGrdOceanEnum,
    392392        SolidearthSettingsOceanAreaScalingEnum,
    393393        StochasticForcingCovarianceEnum,
     394        StochasticForcingDefaultDimensionEnum,
    394395        StochasticForcingDimensionsEnum,
    395396        StochasticForcingFieldsEnum,
    396397        StochasticForcingIsStochasticForcingEnum,
    … …  
    421422        SmbAccurefEnum,
    422423        SmbAdThreshEnum,
    423424        SmbAutoregressionInitialTimeEnum,
    424    SmbAutoregressionNoiseEnum,
    425425        SmbAutoregressionTimestepEnum,
    426426   SmbAutoregressiveOrderEnum,
    427427        SmbAveragingEnum,
    … …  
    497497        StressbalanceRestolEnum,
    498498        StressbalanceRiftPenaltyThresholdEnum,
    499499        StressbalanceShelfDampeningEnum,
    500    ThermalforcingAutoregressionNoiseEnum,
    501500        ThermalIsdrainicecolumnEnum,
    502501        ThermalIsdynamicbasalspcEnum,
    503502        ThermalIsenthalpyEnum,
    … …  
    597596        BaseOldEnum,
    598597        BaseSlopeXEnum,
    599598        BaseSlopeYEnum,
     599        BaselineCalvingCalvingrateEnum,
    600600        BedEnum,
    601601        BedGRDEnum,
    602602        BedEastEnum,
    … …  
    890890        SmbAccumulationEnum,
    891891        SmbAdiffiniEnum,
    892892        SmbAiniEnum,
     893   SmbAutoregressionNoiseEnum,
    893894        SmbBasinsIdEnum,
    894895        SmbBMaxEnum,
    895896        SmbBMinEnum,
    … …  
    993994        SolidearthExternalDisplacementNorthRateEnum,
    994995        SolidearthExternalDisplacementUpRateEnum,
    995996        SolidearthExternalGeoidRateEnum,
     997        StochasticForcingDefaultIdEnum,
    996998        StrainRateeffectiveEnum,
    997999        StrainRateparallelEnum,
    9981000        StrainRateperpendicularEnum,
    … …  
    10281030        TemperaturePDDEnum,
    10291031        TemperaturePicardEnum,
    10301032        TemperatureSEMICEnum,
     1033   ThermalforcingAutoregressionNoiseEnum,
    10311034        ThermalforcingValuesAutoregressionEnum,
    10321035        ThermalSpctemperatureEnum,
    10331036        ThicknessAbsGradientEnum,

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

     
    399399                case SolidearthSettingsGrdOceanEnum : return "SolidearthSettingsGrdOcean";
    400400                case SolidearthSettingsOceanAreaScalingEnum : return "SolidearthSettingsOceanAreaScaling";
    401401                case StochasticForcingCovarianceEnum : return "StochasticForcingCovariance";
     402                case StochasticForcingDefaultDimensionEnum : return "StochasticForcingDefaultDimension";
    402403                case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions";
    403404                case StochasticForcingFieldsEnum : return "StochasticForcingFields";
    404405                case StochasticForcingIsStochasticForcingEnum : return "StochasticForcingIsStochasticForcing";
    … …  
    429430                case SmbAccurefEnum : return "SmbAccuref";
    430431                case SmbAdThreshEnum : return "SmbAdThresh";
    431432                case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime";
    432                 case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
    433433                case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep";
    434434                case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder";
    435435                case SmbAveragingEnum : return "SmbAveraging";
    … …  
    505505                case StressbalanceRestolEnum : return "StressbalanceRestol";
    506506                case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold";
    507507                case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening";
    508                 case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise";
    509508                case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn";
    510509                case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc";
    511510                case ThermalIsenthalpyEnum : return "ThermalIsenthalpy";
    … …  
    603602                case BaseOldEnum : return "BaseOld";
    604603                case BaseSlopeXEnum : return "BaseSlopeX";
    605604                case BaseSlopeYEnum : return "BaseSlopeY";
     605                case BaselineCalvingCalvingrateEnum : return "BaselineCalvingCalvingrate";
    606606                case BedEnum : return "Bed";
    607607                case BedGRDEnum : return "BedGRD";
    608608                case BedEastEnum : return "BedEast";
    … …  
    896896                case SmbAccumulationEnum : return "SmbAccumulation";
    897897                case SmbAdiffiniEnum : return "SmbAdiffini";
    898898                case SmbAiniEnum : return "SmbAini";
     899                case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
    899900                case SmbBasinsIdEnum : return "SmbBasinsId";
    900901                case SmbBMaxEnum : return "SmbBMax";
    901902                case SmbBMinEnum : return "SmbBMin";
    … …  
    998999                case SolidearthExternalDisplacementNorthRateEnum : return "SolidearthExternalDisplacementNorthRate";
    9991000                case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate";
    10001001                case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate";
     1002                case StochasticForcingDefaultIdEnum : return "StochasticForcingDefaultId";
    10011003                case StrainRateeffectiveEnum : return "StrainRateeffective";
    10021004                case StrainRateparallelEnum : return "StrainRateparallel";
    10031005                case StrainRateperpendicularEnum : return "StrainRateperpendicular";
    … …  
    10331035                case TemperaturePDDEnum : return "TemperaturePDD";
    10341036                case TemperaturePicardEnum : return "TemperaturePicard";
    10351037                case TemperatureSEMICEnum : return "TemperatureSEMIC";
     1038                case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise";
    10361039                case ThermalforcingValuesAutoregressionEnum : return "ThermalforcingValuesAutoregression";
    10371040                case ThermalSpctemperatureEnum : return "ThermalSpctemperature";
    10381041                case ThicknessAbsGradientEnum : return "ThicknessAbsGradient";

  • ../trunk-jpl/src/c/shared/Enum/Enumjl.vim

     
    169169syn keyword juliaConstC FrictionThresholdSpeedEnum
    170170syn keyword juliaConstC FrictionVoidRatioEnum
    171171syn keyword juliaConstC FrontalForcingsBasinIcefrontAreaEnum
     172syn keyword juliaConstC FrontalForcingsAutoregressionInitialTimeEnum
     173syn keyword juliaConstC FrontalForcingsAutoregressionTimestepEnum
     174syn keyword juliaConstC FrontalForcingsAutoregressiveOrderEnum
     175syn keyword juliaConstC FrontalForcingsBeta0Enum
     176syn keyword juliaConstC FrontalForcingsBeta1Enum
    172177syn keyword juliaConstC FrontalForcingsNumberofBasinsEnum
    173178syn keyword juliaConstC FrontalForcingsParamEnum
     179syn keyword juliaConstC FrontalForcingsPhiEnum
    174180syn keyword juliaConstC GrdModelEnum
    175181syn keyword juliaConstC GroundinglineFrictionInterpolationEnum
    176182syn keyword juliaConstC GroundinglineMeltInterpolationEnum
    … …  
    383389syn keyword juliaConstC SolidearthSettingsMaxiterEnum
    384390syn keyword juliaConstC SolidearthSettingsGrdOceanEnum
    385391syn keyword juliaConstC SolidearthSettingsOceanAreaScalingEnum
     392syn keyword juliaConstC StochasticForcingCovarianceEnum
     393syn keyword juliaConstC StochasticForcingDefaultDimensionEnum
     394syn keyword juliaConstC StochasticForcingDimensionsEnum
     395syn keyword juliaConstC StochasticForcingFieldsEnum
     396syn keyword juliaConstC StochasticForcingIsStochasticForcingEnum
     397syn keyword juliaConstC StochasticForcingNumFieldsEnum
     398syn keyword juliaConstC StochasticForcingRandomflagEnum
    386399syn keyword juliaConstC RotationalPolarMoiEnum
    387400syn keyword juliaConstC SolidearthSettingsReltolEnum
    388401syn keyword juliaConstC SealevelchangeRequestedOutputsEnum
    … …  
    413426syn keyword juliaConstC SmbAveragingEnum
    414427syn keyword juliaConstC SmbBeta0Enum
    415428syn keyword juliaConstC SmbBeta1Enum
    416 syn keyword juliaConstC SmbCovmatEnum
    417429syn keyword juliaConstC SmbDesfacEnum
    418430syn keyword juliaConstC SmbDpermilEnum
    419431syn keyword juliaConstC SmbDsnowIdxEnum
    … …  
    423435syn keyword juliaConstC SmbDenIdxEnum
    424436syn keyword juliaConstC SmbDtEnum
    425437syn keyword juliaConstC SmbEnum
     438syn keyword juliaConstC SmbEIdxEnum
    426439syn keyword juliaConstC SmbFEnum
    427440syn keyword juliaConstC SmbInitDensityScalingEnum
    428441syn keyword juliaConstC SmbIsaccumulationEnum
    … …  
    431444syn keyword juliaConstC SmbIsd18opdEnum
    432445syn keyword juliaConstC SmbIsdelta18oEnum
    433446syn keyword juliaConstC SmbIsdensificationEnum
     447syn keyword juliaConstC SmbIsdeltaLWupEnum
    434448syn keyword juliaConstC SmbIsfirnwarmingEnum
    435449syn keyword juliaConstC SmbIsgraingrowthEnum
    436450syn keyword juliaConstC SmbIsmeltEnum
    … …  
    446460syn keyword juliaConstC SmbNumRequestedOutputsEnum
    447461syn keyword juliaConstC SmbPfacEnum
    448462syn keyword juliaConstC SmbPhiEnum
    449 syn keyword juliaConstC SmbRandomflagEnum
    450463syn keyword juliaConstC SmbRdlEnum
    451464syn keyword juliaConstC SmbRequestedOutputsEnum
    452465syn keyword juliaConstC SmbRlapsEnum
    … …  
    459472syn keyword juliaConstC SmbSwIdxEnum
    460473syn keyword juliaConstC SmbT0dryEnum
    461474syn keyword juliaConstC SmbT0wetEnum
     475syn keyword juliaConstC SmbTeThreshEnum
    462476syn keyword juliaConstC SmbTdiffEnum
    463477syn keyword juliaConstC SmbThermoDeltaTScalingEnum
    464478syn keyword juliaConstC SmbTemperaturesReconstructedYearsEnum
    … …  
    579593syn keyword juliaConstC BaseOldEnum
    580594syn keyword juliaConstC BaseSlopeXEnum
    581595syn keyword juliaConstC BaseSlopeYEnum
     596syn keyword juliaConstC BaselineCalvingCalvingrateEnum
    582597syn keyword juliaConstC BedEnum
    583598syn keyword juliaConstC BedGRDEnum
    584599syn keyword juliaConstC BedEastEnum
    … …  
    872887syn keyword juliaConstC SmbAccumulationEnum
    873888syn keyword juliaConstC SmbAdiffiniEnum
    874889syn keyword juliaConstC SmbAiniEnum
     890syn keyword juliaConstC SmbAutoregressionNoiseEnum
    875891syn keyword juliaConstC SmbBasinsIdEnum
    876892syn keyword juliaConstC SmbBMaxEnum
    877893syn keyword juliaConstC SmbBMinEnum
    … …  
    893909syn keyword juliaConstC SmbDailywindspeedEnum
    894910syn keyword juliaConstC SmbDiniEnum
    895911syn keyword juliaConstC SmbDlwrfEnum
     912syn keyword juliaConstC SmbDulwrfValueEnum
    896913syn keyword juliaConstC SmbDswrfEnum
    897914syn keyword juliaConstC SmbDswdiffrfEnum
    898915syn keyword juliaConstC SmbDzAddEnum
    … …  
    973990syn keyword juliaConstC SolidearthExternalDisplacementNorthRateEnum
    974991syn keyword juliaConstC SolidearthExternalDisplacementUpRateEnum
    975992syn keyword juliaConstC SolidearthExternalGeoidRateEnum
     993syn keyword juliaConstC StochasticForcingDefaultIdEnum
    976994syn keyword juliaConstC StrainRateeffectiveEnum
    977995syn keyword juliaConstC StrainRateparallelEnum
    978996syn keyword juliaConstC StrainRateperpendicularEnum
    … …  
    10081026syn keyword juliaConstC TemperaturePDDEnum
    10091027syn keyword juliaConstC TemperaturePicardEnum
    10101028syn keyword juliaConstC TemperatureSEMICEnum
     1029syn keyword juliaConstC ThermalforcingAutoregressionNoiseEnum
     1030syn keyword juliaConstC ThermalforcingValuesAutoregressionEnum
    10111031syn keyword juliaConstC ThermalSpctemperatureEnum
    10121032syn keyword juliaConstC ThicknessAbsGradientEnum
    10131033syn keyword juliaConstC ThicknessAbsMisfitEnum
    … …  
    12541274syn keyword juliaConstC FreeSurfaceTopAnalysisEnum
    12551275syn keyword juliaConstC FrontalForcingsDefaultEnum
    12561276syn keyword juliaConstC FrontalForcingsRignotEnum
     1277syn keyword juliaConstC FrontalForcingsRignotAutoregressionEnum
    12571278syn keyword juliaConstC FsetEnum
    12581279syn keyword juliaConstC FullMeltOnPartiallyFloatingEnum
    12591280syn keyword juliaConstC GLheightadvectionAnalysisEnum

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

     
    408408              else if (strcmp(name,"SolidearthSettingsGrdOcean")==0) return SolidearthSettingsGrdOceanEnum;
    409409              else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum;
    410410              else if (strcmp(name,"StochasticForcingCovariance")==0) return StochasticForcingCovarianceEnum;
     411              else if (strcmp(name,"StochasticForcingDefaultDimension")==0) return StochasticForcingDefaultDimensionEnum;
    411412              else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum;
    412413              else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum;
    413414              else if (strcmp(name,"StochasticForcingIsStochasticForcing")==0) return StochasticForcingIsStochasticForcingEnum;
    … …  
    438439              else if (strcmp(name,"SmbAccuref")==0) return SmbAccurefEnum;
    439440              else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
    440441              else if (strcmp(name,"SmbAutoregressionInitialTime")==0) return SmbAutoregressionInitialTimeEnum;
    441               else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum;
    442442              else if (strcmp(name,"SmbAutoregressionTimestep")==0) return SmbAutoregressionTimestepEnum;
    443443              else if (strcmp(name,"SmbAutoregressiveOrder")==0) return SmbAutoregressiveOrderEnum;
    444444              else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum;
    … …  
    517517              else if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum;
    518518              else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum;
    519519              else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum;
    520               else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum;
    521520              else if (strcmp(name,"ThermalIsdrainicecolumn")==0) return ThermalIsdrainicecolumnEnum;
    522521              else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
    523522              else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
    … …  
    615614              else if (strcmp(name,"BaseOld")==0) return BaseOldEnum;
    616615              else if (strcmp(name,"BaseSlopeX")==0) return BaseSlopeXEnum;
    617616              else if (strcmp(name,"BaseSlopeY")==0) return BaseSlopeYEnum;
     617              else if (strcmp(name,"BaselineCalvingCalvingrate")==0) return BaselineCalvingCalvingrateEnum;
    618618              else if (strcmp(name,"Bed")==0) return BedEnum;
    619619              else if (strcmp(name,"BedGRD")==0) return BedGRDEnum;
    620620              else if (strcmp(name,"BedEast")==0) return BedEastEnum;
    … …  
    917917              else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum;
    918918              else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum;
    919919              else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
     920              else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum;
    920921              else if (strcmp(name,"SmbBasinsId")==0) return SmbBasinsIdEnum;
    921922              else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
    922923              else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
    … …  
    996997              else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
    997998              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
    998999              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    999               else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
     1003              if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
     1004              else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
    10041005              else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
    10051006              else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
    10061007              else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum;
    … …  
    10221023              else if (strcmp(name,"SolidearthExternalDisplacementNorthRate")==0) return SolidearthExternalDisplacementNorthRateEnum;
    10231024              else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum;
    10241025              else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
     1026              else if (strcmp(name,"StochasticForcingDefaultId")==0) return StochasticForcingDefaultIdEnum;
    10251027              else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
    10261028              else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
    10271029              else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum;
    … …  
    10571059              else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum;
    10581060              else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
    10591061              else if (strcmp(name,"TemperatureSEMIC")==0) return TemperatureSEMICEnum;
     1062              else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum;
    10601063              else if (strcmp(name,"ThermalforcingValuesAutoregression")==0) return ThermalforcingValuesAutoregressionEnum;
    10611064              else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum;
    10621065              else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
    … …  
    11171120              else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
    11181121              else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
    11191122              else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
    1120               else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
    1121               else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
    1122               else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
     1126              if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
     1127              else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
     1128              else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum;
     1129              else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
    11271130              else if (strcmp(name,"Outputdefinition30")==0) return Outputdefinition30Enum;
    11281131              else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
    11291132              else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum;
    … …  
    12401243              else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
    12411244              else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
    12421245              else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
    1243               else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
    1244               else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum;
    1245               else if (strcmp(name,"Channel")==0) return ChannelEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
     1249              if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
     1250              else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum;
     1251              else if (strcmp(name,"Channel")==0) return ChannelEnum;
     1252              else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
    12501253              else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum;
    12511254              else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum;
    12521255              else if (strcmp(name,"Closed")==0) return ClosedEnum;
    … …  
    13631366              else if (strcmp(name,"IntParam")==0) return IntParamEnum;
    13641367              else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
    13651368              else if (strcmp(name,"Inputs")==0) return InputsEnum;
    1366               else if (strcmp(name,"Internal")==0) return InternalEnum;
    1367               else if (strcmp(name,"Intersect")==0) return IntersectEnum;
    1368               else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"J")==0) return JEnum;
     1372              if (strcmp(name,"Internal")==0) return InternalEnum;
     1373              else if (strcmp(name,"Intersect")==0) return IntersectEnum;
     1374              else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
     1375              else if (strcmp(name,"J")==0) return JEnum;
    13731376              else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum;
    13741377              else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum;
    13751378              else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
    … …  
    14861489              else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
    14871490              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
    14881491              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    1489               else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
    1490               else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
    1491               else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
    14921492         else stage=13;
    14931493   }
    14941494   if(stage==13){
    1495               if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
     1495              if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
     1496              else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
     1497              else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
     1498              else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
    14961499              else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
    14971500              else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;
    14981501              else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;

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

     
    4848%Stochastic forcing
    4949md.stochasticforcing.isstochasticforcing = 1;
    5050md.stochasticforcing.fields              = [{'SMBautoregression'}];
    51 md.stochasticforcing.dimensions          = [md.smb.num_basins]; %dimension of each field
    5251md.stochasticforcing.covariance          = [[0.15 0.08 -0.02];[0.08 0.12 -0.05];[-0.02 -0.05 0.1]]; %global covariance among- and between-fields
    5352md.stochasticforcing.randomflag          = 0; %fixed random seeds
    5453

  • ../trunk-jpl/test/NightlyRun/test257.py

     
    5555# Stochastic forcing
    5656md.stochasticforcing.isstochasticforcing = 1
    5757md.stochasticforcing.fields = ['SMBautoregression']
    58 md.stochasticforcing.dimensions = [md.smb.num_basins] # dimension of each field
    5958md.stochasticforcing.covariance = np.array([[0.15, 0.08, -0.02], [0.08, 0.12, -0.05], [-0.02, -0.05, 0.1]]) # global covariance among- and between-fields
    6059md.stochasticforcing.randomflag = 0 # fixed random seeds
    6160

  • ../trunk-jpl/test/NightlyRun/test543.m

     
    3838%Stochastic forcing
    3939md.stochasticforcing.isstochasticforcing = 1;
    4040md.stochasticforcing.fields              = [{'FrontalForcingsRignotAutoregression'}];
    41 md.stochasticforcing.dimensions          = [md.frontalforcings.num_basins]; %dimension of each field
    4241md.stochasticforcing.covariance          = 1e-4*[[1.5,0.5];[0.5,0.4]]; %global covariance among- and between-fields
    4342md.stochasticforcing.randomflag          = 0; %determines true/false randomness
    4443

  • ../trunk-jpl/test/NightlyRun/test543.py

     
    4848# Stochastic forcing
    4949md.stochasticforcing.isstochasticforcing = 1
    5050md.stochasticforcing.fields = ['FrontalForcingsRignotAutoregression']
    51 md.stochasticforcing.dimensions = [md.frontalforcings.num_basins] # dimension of each field
    5251md.stochasticforcing.covariance = 1e-4 * np.array([[1.5, 0.5], [0.5, 0.4]]) # global covariance among- and between-fields
    5352md.stochasticforcing.randomflag = 0 # determines true/false randomness
    5453

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

     
    77        properties (SetAccess=public)
    88                isstochasticforcing  = 0;
    99                fields               = NaN;
    10                 dimensions           = NaN;
     10                defaultdimension     = 0;
     11      default_id           = NaN;
    1112                covariance           = NaN;
    1213                randomflag           = 1;
    1314        end
    … …  
    3334                        if ~self.isstochasticforcing, return; end
    3435
    3536                        num_fields = numel(self.fields);
    36                         size_tot   = sum(self.dimensions);
     37                       
     38                        %Check that covariance matrix is positive definite
     39                        try
     40                                chol(self.covariance);
     41                        catch
     42                                error('md.stochasticforcing.covariance is not positive definite');
     43                        end
    3744
     45                        %Check that all fields agree with the corresponding md class and if any field needs the default params   
     46         checkdefaults = false; %need to check defaults only if one of the field does not have its own dimensionality
     47                        for field=self.fields
     48            %Checking agreement of classes
     49            if(contains(field,'SMB'))
     50               if~(isequal(class(md.smb),char(field)))
     51                  error('md.smb does not agree with stochasticforcing field %s', char(field));
     52               end
     53            end
     54            if(contains(field,'frontalforcings'))
     55               if~(isequal(class(md.frontalforcings),char(field)))
     56                  error('md.frontalforcings does not agree with stochasticforcing field %s', char(field));
     57               end
     58            end
     59            %Checking for specific dimensions
     60            if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression'))
     61               checkdefaults = true; %field with non-specific dimensionality
     62            end
     63         end
     64
     65                        %Retrieve sum of all the field dimensionalities
     66                        size_tot   = self.defaultdimension*num_fields;
     67                        indSMBar = -1; %about to check for index of SMBautoregression
     68                        indTFar  = -1; %about to check for index of FrontalForcingsRignotAutoregression
     69                        if any(contains(self.fields,'SMBautoregression'))
     70            size_tot = size_tot-self.defaultdimension+md.smb.num_basins;
     71                                indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with TFar timestep (08Nov2021)
     72         end
     73         if any(contains(self.fields,'FrontalForcingsRignotAutoregression'))
     74            size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins;
     75                                indTFar  = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with SMBar timestep (08Nov2021)
     76         end
     77
     78                        if(indSMBar~=-1 && indTFar~=-1) %both autoregressive models are used: check autoregressive time step consistency
     79                                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))))~=0)
     80                                        error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance');
     81                                end
     82                        end
     83
    3884                        md = checkfield(md,'fieldname','stochasticforcing.isstochasticforcing','values',[0 1]);
    39                         md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings()); %VV check here 'cell' (19Oct2021)
    40                         md = checkfield(md,'fieldname','stochasticforcing.dimensions','NaN',1,'Inf',1,'>',0,'size',[1,num_fields]); %specific dimension for each field
     85                        md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings());
    4186                        md = checkfield(md,'fieldname','stochasticforcing.covariance','NaN',1,'Inf',1,'size',[size_tot,size_tot]); %global covariance matrix
    4287                        md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]);
    43 
    44                         %Check that all fields agree with the corresponding md class
    45                         for field=self.fields
    46                                 if(contains(field,'SMB'))
    47                                         if~(isequal(class(md.smb),char(field)))
    48                                                 error('md.smb does not agree with stochasticforcing field %s', char(field));
    49                                         end
    50                                 end
    51                                 if(contains(field,'frontalforcings'))
    52                                         if~(isequal(class(md.frontalforcings),char(field)))
    53                                                 error('md.frontalforcings does not agree with stochasticforcing field %s', char(field));
    54                                         end
    55                                 end
    56                         end
     88                        if(checkdefaults) %need to check the defaults
     89            md = checkfield(md,'fieldname','stochasticforcing.defaultdimension','numel',1,'NaN',1,'Inf',1,'>',0);
     90            md = checkfield(md,'fieldname','stochasticforcing.default_id','Inf',1,'>=',0,'<=',self.defaultdimension,'size',[md.mesh.numberofelements,1]);
     91         end
    5792                end % }}}
    5893                function disp(self) % {{{
    5994                        disp(sprintf('   stochasticforcing parameters:'));
    6095                        fielddisplay(self,'isstochasticforcing','is stochasticity activated?');
    6196                        fielddisplay(self,'fields','fields with stochasticity applied, ex: {''SMBautoregression''}, or {''FrontalForcingsRignotAutoregression''}');
    62                         fielddisplay(self,'dimensions','dimensionality of each field');
     97                        fielddisplay(self,'defaultdimension','dimensionality of the noise terms (does not apply to fields with their specific dimension)');
     98         fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)');
    6399                        fielddisplay(self,'covariance','covariance matrix for within- and between-fields covariance (units must be squared field units)');
    64100                        fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)');
    65101                        disp('Available fields:');
    … …  
    76112                        if ~self.isstochasticforcing
    77113                                return
    78114                        else
     115
     116                                %Retrieve dimensionality of each field
     117                                dimensions = self.defaultdimension*ones(1,num_fields);
     118                                ind = 1;
     119                                for field=self.fields
     120                                        %Checking for specific dimensions
     121                                        if(strcmp(field,'SMBautoregression'))
     122                                                dimensions(ind) = md.smb.num_basins;
     123                                        end
     124                                        if(strcmp(field,'FrontalForcingsRignotAutoregression'))
     125                                                dimensions(ind) = md.frontalforcings.num_basins;
     126                                        end
     127                                        ind = ind+1;
     128                                end
     129
    79130                                %Scaling covariance matrix (scale column-by-column and row-by-row)
    80                                 scaledfields = {'SMBautoregression'}; %list of fields that need scaling *1/yts
     131                                scaledfields = {'DefaultCalving','SMBautoregression'}; %list of fields that need scaling *1/yts
    81132                                tempcovariance = self.covariance; %copy of covariance to avoid writing back in member variable
    82133                                for i=1:num_fields
    83134                                        if any(strcmp(scaledfields,self.fields(i)))
    84                                                 inds = [1+sum(self.dimensions(1:i-1)):1:sum(self.dimensions(1:i))];
     135                                                inds = [1+sum(dimensions(1:i-1)):1:sum(dimensions(1:i))];
    85136                                                for row=inds %scale rows corresponding to scaled field
    86137                                                        tempcovariance(row,:) = 1./yts*tempcovariance(row,:);
    87138                                                end
    … …  
    92143                                end
    93144                                WriteData(fid,prefix,'data',num_fields,'name','md.stochasticforcing.num_fields','format','Integer');
    94145                                WriteData(fid,prefix,'object',self,'fieldname','fields','format','StringArray');
    95                                 WriteData(fid,prefix,'object',self,'fieldname','dimensions','format','IntMat');
     146                                WriteData(fid,prefix,'data',dimensions,'name','md.stochasticforcing.dimensions','format','IntMat');
     147            WriteData(fid,prefix,'object',self,'fieldname','default_id','data',self.default_id-1,'format','IntMat','mattype',2); %0-indexed
     148                                WriteData(fid,prefix,'object',self,'fieldname','defaultdimension','format','Integer');
    96149                                WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat');
    97150                                WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean');
    98151                        end
    … …  
    104157   % by the class md.stochasticforcing
    105158
    106159   list = {...
    107       'SMBautoregression',...
    108       'FrontalForcingsRignotAutoregression'
     160      'DefaultCalving',...
     161                'FrontalForcingsRignotAutoregression',...
     162      'SMBautoregression'
    109163      };
    110164end % }}}

  • ../trunk-jpl/src/m/classes/stochasticforcing.py

     
    1414
    1515    def __init__(self, *args):  # {{{
    1616        self.isstochasticforcing = 0
    17         self.fields = np.nan
    18         self.dimensions = np.nan
    19         self.covariance = np.nan
    20         self.randomflag = 1
     17        self.fields              = np.nan
     18        self.defaultdimension    = 0
     19        self.default_id          = np.nan
     20        self.covariance          = np.nan
     21        self.randomflag          = 1
    2122
    2223        if len(args) == 0:
    2324            self.setdefaultparameters()
    … …  
    2829        s = '   stochasticforcing parameters:\n'
    2930        s += '{}\n'.format(fielddisplay(self, 'isstochasticforcing', 'is stochasticity activated?'))
    3031        s += '{}\n'.format(fielddisplay(self, 'fields', 'fields with stochasticity applied, ex: [\'SMBautoregression\'], or [\'FrontalForcingsRignotAutoregression\']'))
     32        s += '{}\n'.format(fielddisplay(self, 'defaultdimension', 'dimensionality of the noise terms (does not apply to fields with their specific dimension)'))
     33        s += '{}\n'.format(fielddisplay(self, 'default_id', 'id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)'))
    3134        s += '{}\n'.format(fielddisplay(self, 'covariance', 'covariance matrix for within- and between-fields covariance (units must be squared field units)'))
    3235        s += '{}\n'.format(fielddisplay(self, 'randomflag', 'whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'))
    3336        s += 'Available fields:\n'
     37        s += '   DefaultCalving\n'
    3438        s += '   SMBautoregression\n'
    3539        s += '   FrontalForcingsRignotAutoregression (thermal forcing)\n'
    3640        return s
    … …  
    5054            return md
    5155
    5256        num_fields  = numel(self.fields)
    53         size_tot    = np.sum(self.dimensions)
    5457
    55         md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1])
    56         md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', supportedstochforcings()) # VV check here 'cell' (19Oct2021)
    57         md = checkfield(md, 'fieldname', 'stochasticforcing.dimensions', 'NaN', 1, 'Inf', 1, '>', 0, 'size', [num_fields]) # specific dimension for each field; NOTE: As opposed to MATLAB implementation, pass list
    58         md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix
    59         md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1])
     58        #Check that covariance matrix is positive definite
     59        try:
     60            np.linalg.cholesky(self.covariance);
     61        except:
     62            error('md.stochasticforcing.covariance is not positive definite');
    6063
    6164        # Check that all fields agree with the corresponding md class
     65        checkdefaults = False
    6266        for field in self.fields:
    6367            if (contains(field, 'SMB')):
    6468                if not (type(md.smb) == field):
    … …  
    6670            if (contains(field, 'frontalforcings')):
    6771                if not (type(md.frontalforcings) == field):
    6872                    error('md.frontalforcings does not agree with stochasticforcing field {}'.format(field))
     73            #Checking for specific dimensions
     74            if (field!='SMBautoregression' or field!='FrontalForcingsRignotAutoregression'):
     75                checkdefaults = True #field with non-specific dimensionality
     76
     77        #Retrieve sum of all the field dimensionalities
     78        size_tot   = self.defaultdimension*num_fields;
     79        indSMBar = -1; #about to check for index of SMBautoregression
     80        indTFar  = -1; #about to check for index of FrontalForcingsRignotAutoregression
     81        if ('SMBautoregression' in self.fields):
     82            size_tot = size_tot-self.defaultdimension+md.smb.num_basins
     83            indSMBar = np.where(self.fields=='SMBautoregression')[0][0]
     84        if ('FrontalForcingsRignotAutoregression' in self.fields):
     85            size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins
     86            indSMBar = np.where(self.fields=='FrontalForcingsRignotAutoregression')[0][0]
     87        if (indSMBar!=-1 and indTFar!=-1):
     88            if((md.smb.ar_timestep!=md.frontalforcings.ar_timestep) and np.any(self.covariance[np.sum(self.dimensions[0:indSMBar]).astype(int):np.sum(self.dimensions[0:indSMBar+1]).astype(int),np.sum(self.dimensions[0:indTFar]).astype(int):np.sum(self.dimensions[0:indTFar+1]).astype(int)]!=0)):
     89                error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance')
     90
     91        md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1])
     92        md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', supportedstochforcings())
     93        #md = checkfield(md, 'fieldname', 'stochasticforcing.dimensions', 'NaN', 1, 'Inf', 1, '>', 0, 'size', [num_fields]) # specific dimension for each field; NOTE: As opposed to MATLAB implementation, pass list
     94        md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix
     95        md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1])
     96        if (checkdefaults):
     97            md = checkfield(md, 'fieldname', 'stochasticforcing.defaultdimension', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0)
     98            md = checkfield(md, 'fieldname', 'stochasticforcing.default_id', 'Inf', 1, '>=', 0, '<=', self.defaultdimension, 'size', [md.mesh.numberofelements,1])
    6999        return md
    70100    # }}}
    71101
    … …  
    82112        if not self.isstochasticforcing:
    83113            return md
    84114        else:
     115            #Retrieve dimensionality of each field
     116            dimensions = self.defaultdimension*np.ones(num_fields)
     117            for ind,field in enumerate(self.fields):
     118                #Checking for specific dimensions
     119                if (field=='SMBautoregression'):
     120                    dimensions[ind] = md.smb.num_basins
     121                if (field=='FrontalForcingsRignotAutoregression'):
     122                    dimensions[ind] = md.frontalforcings.num_basins
     123
    85124            # Scaling covariance matrix (scale column-by-column and row-by-row)
    86             scaledfields = ['SMBautoregression'] # list of fields that need scaling * 1/yts
    87             tempcovariance = self.covariance
     125            scaledfields = ['DefaultCalving','SMBautoregression'] # list of fields that need scaling * 1/yts
     126            tempcovariance = np.copy(self.covariance)
    88127            for i in range(num_fields):
    89128                if self.fields[i] in scaledfields:
    90129                    inds = range(int(np.sum(self.dimensions[0:i])), int(np.sum(self.dimensions[0:i + 1])))
    91130                    for row in inds: # scale rows corresponding to scaled field
    92                         tempcovariance[row, :] = 1 / yts * self.covariance[row, :]
     131                        tempcovariance[row, :] = 1 / yts * tempcovariance[row, :]
    93132                    for col in inds: # scale columns corresponding to scaled field
    94                         tempcovariance[:, col] = 1 / yts * self.covariance[:, col]
     133                        tempcovariance[:, col] = 1 / yts * tempcovariance[:, col]
    95134            WriteData(fid, prefix, 'data', num_fields, 'name', 'md.stochasticforcing.num_fields', 'format', 'Integer')
    96135            WriteData(fid, prefix, 'object', self, 'fieldname', 'fields', 'format', 'StringArray')
    97             WriteData(fid, prefix, 'object', self, 'fieldname','dimensions', 'format', 'IntMat')
     136            WriteData(fid, prefix, 'data', 'dimensions', 'name', 'md.stochasticforcing.dimensions', 'format', 'IntMat')
     137            WriteData(fid, prefix, 'object', self, 'fieldname', 'default_id', 'format', 'IntMat') #12Nov2021 make sure this is zero-indexed!
     138            WriteData(fid, prefix, 'object', self, 'fieldname', 'defaultdimension', 'format', 'Integer')
    98139            WriteData(fid, prefix, 'data', tempcovariance, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat')
    99140            WriteData(fid, prefix, 'object', self, 'fieldname', 'randomflag', 'format', 'Boolean')
    100141    # }}}
    … …  
    103144    """ Defines list of fields supported  by the class stochasticforcings
    104145    """
    105146    return [
    106         'SMBautoregression',
    107         'FrontalForcingsRignotAutoregression'
     147        'DefaultCalving',
     148        'FrontalForcingsRignotAutoregression',
     149        'SMBautoregression'
    108150    ]
Note: See TracBrowser for help on using the repository browser.