Changeset 27354

Timestamp:

11/04/22 04:51:01 (2 years ago)

Author:

vverjans

Message:

NEW added autoregression and stochastic capabilities for subglacial discharge

Location:

issm/trunk-jpl

Files:

• src/c/analyses/LevelsetAnalysis.cpp (modified) (3 diffs)
• src/c/classes/Elements/Element.cpp (modified) (6 diffs)
• src/c/classes/Elements/Element.h (modified) (1 diff)
• src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (modified) (2 diffs)
• src/c/modules/FrontalForcingsx/FrontalForcingsx.h (modified) (1 diff)
• src/c/modules/StochasticForcingx/StochasticForcingx.cpp (modified) (3 diffs)
• src/c/shared/Enum/Enum.vim (modified) (4 diffs)
• src/c/shared/Enum/EnumDefinitions.h (modified) (4 diffs)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (4 diffs)
• src/c/shared/Enum/Enumjl.vim (modified) (4 diffs)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (28 diffs)
• src/m/classes/frontalforcingsrignotarma.m (modified) (7 diffs)
• src/m/classes/frontalforcingsrignotarma.py (modified) (6 diffs)
• src/m/classes/stochasticforcing.m (modified) (6 diffs)
• src/m/classes/stochasticforcing.py (modified) (6 diffs)
• test/Archives/Archive543.arch (modified) ( previous)
• test/Archives/Archive544.arch (modified) ( previous)
• test/NightlyRun/test543.m (modified) (3 diffs)
• test/NightlyRun/test543.py (modified) (4 diffs)

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

@@ -153 +153 @@
          /*Retrieve thermal forcing only in the case of non-arma FrontalForcingsRignot*/
          iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.thermalforcing",ThermalForcingEnum);
-         /* Do not break here, still retrieve basin_ID,subglacial_discharge, etc.*/
-      case FrontalForcingsRignotarmaEnum:
          iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.basin_id",FrontalForcingsBasinIdEnum);
          iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.subglacial_discharge",FrontalForcingsSubglacialDischargeEnum);
-         break; 
+                        break;
+                case FrontalForcingsRignotarmaEnum:
+                        bool isdischargearma;
+                        iomodel->FindConstant(&isdischargearma,"md.frontalforcings.isdischargearma");
+         iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.basin_id",FrontalForcingsBasinIdEnum);
+         if(isdischargearma==false) iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.subglacial_discharge",FrontalForcingsSubglacialDischargeEnum);
+                        break;  
                 default:
                         _error_("Frontal forcings"<<EnumToStringx(melt_parameterization)<<" not supported yet");
@@ -245 +249 @@
                         break;
                 case FrontalForcingsRignotarmaEnum:
-                        /*Retrieve autoregressive parameters*/
+                        parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_basins",FrontalForcingsNumberofBasinsEnum));
+                        /*Retrieve thermal forcing parameters*/
                         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_params",FrontalForcingsNumberofParamsEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_breaks",FrontalForcingsNumberofBreaksEnum));
@@ -273 +278 @@
          parameters->AddObject(new DoubleMatParam(FrontalForcingsARMAmonthtrendsEnum,transparam,M,N));
          xDelete<IssmDouble>(transparam);
-                        /*Do not break here, generic FrontalForcingsRignot parameters still to be retrieved*/
+                        parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.isdischargearma",FrontalForcingsIsDischargeARMAEnum));
+                        /*Retrieve subglacial discharge parameters */
+                        bool isdischargearma;
+                        parameters->FindParam(&isdischargearma,FrontalForcingsIsDischargeARMAEnum);
+                        if(isdischargearma==true){
+                                parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.sd_num_params",FrontalForcingsSdNumberofParamsEnum));
+                                parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.sd_num_breaks",FrontalForcingsSdNumberofBreaksEnum));
+           parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.sd_ar_order",FrontalForcingsSdarOrderEnum));
+           parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.sd_ma_order",FrontalForcingsSdmaOrderEnum));
+           parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.sd_arma_timestep",FrontalForcingsSdARMATimestepEnum));
+           iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.sd_polynomialparams");
+           parameters->AddObject(new DoubleMatParam(FrontalForcingsSdpolyparamsEnum,transparam,M,N));
+           xDelete<IssmDouble>(transparam);
+           iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.sd_datebreaks");
+           parameters->AddObject(new DoubleMatParam(FrontalForcingsSddatebreaksEnum,transparam,M,N));
+           xDelete<IssmDouble>(transparam);
+           iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.sd_arlag_coefs");
+           parameters->AddObject(new DoubleMatParam(FrontalForcingsSdarlagcoefsEnum,transparam,M,N));
+           xDelete<IssmDouble>(transparam);
+           iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.sd_malag_coefs");
+           parameters->AddObject(new DoubleMatParam(FrontalForcingsSdmalagcoefsEnum,transparam,M,N));
+           xDelete<IssmDouble>(transparam);
+                                iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.sd_monthlyfrac");
+           parameters->AddObject(new DoubleMatParam(FrontalForcingsSdMonthlyFracEnum,transparam,M,N));
+           xDelete<IssmDouble>(transparam);
+                        }
+                        break;
                 case FrontalForcingsRignotEnum:
                         parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_basins",FrontalForcingsNumberofBasinsEnum));

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

@@ -99 +99 @@
          outenum_type   = BasalforcingsSpatialDeepwaterMeltingRateEnum;
          break;
-   }
+                case(FrontalForcingsSubglacialDischargearmaEnum):
+         arenum_type    = SubglacialdischargeValuesAutoregressionEnum;
+         maenum_type    = SubglacialdischargeValuesMovingaverageEnum;
+         basinenum_type = FrontalForcingsBasinIdEnum;
+         noiseenum_type = SubglacialdischargeARMANoiseEnum;
+         outenum_type   = FrontalForcingsSubglacialDischargeEnum;
+         break;
+        }
 
         /*Get time parameters*/
@@ -126 +133 @@
                         indperiod = 0;
                         for(int i=0;i<numbreaks;i++){
-                                if(tatstep>datebreaks_basin[i]) indperiod = i+1;
+                                if(tatstep>=datebreaks_basin[i]) indperiod = i+1;
                         }
                         /*Compute polynomial with parameters of indperiod*/
@@ -135 +142 @@
                 else for(int j=0;j<numparams;j++) sumpoly[s] = sumpoly[s]+polyparams_basin[j*numperiods]*pow(tatstep,j);
         }
+
         /*Initialze autoregressive and moving-average values at first time step*/
         if(time<=starttime+dt){
@@ -174 +182 @@
                         /*Stochastic variable value*/
          varlist[v] = sumpoly[0]+autoregressionterm+movingaverageterm+noiseterm;
-      }
+      
+                        /*Impose zero-bound*/
+                        if(outenum_type == ThermalForcingEnum || outenum_type == FrontalForcingsSubglacialDischargeEnum) varlist[v] = max(varlist[v],0.0);
+
+                }
+
       /*Update autoregression and moving-average values*/
       IssmDouble* temparrayar = xNew<IssmDouble>(numvertices*arorder);
@@ -2716 +2729 @@
         this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
 
+}/*}}}*/
+void       Element::MonthlyFractionBasin(IssmDouble* monthlyfrac,int enum_type){/*{{{*/
+        
+        /*Variable declaration*/
+   const int numvertices = this->GetNumberOfVertices();
+        int basinid,mindex,mindexnext,basinenum_type,varenum_type,indperiod;
+   IssmDouble time,dt,fracyear,fracyearnext,fracmonth,fracmonthnext,yts; 
+   IssmDouble monthsteps[12]  = {0.,1./12,2./12,3./12,4./12,5./12,6./12,7./12,8./12,9./12,10./12,11./12};
+   IssmDouble* monthlyfrac_b  = xNew<IssmDouble>(12);
+   IssmDouble* monthlyrate_b  = xNew<IssmDouble>(12);
+        IssmDouble* fracdtinmonth  = xNew<IssmDouble>(12);
+        IssmDouble* rateinmonth    = xNew<IssmDouble>(numvertices*12);
+        IssmDouble* varlistinput   = xNew<IssmDouble>(numvertices);
+        IssmDouble* varlist        = xNewZeroInit<IssmDouble>(numvertices);
+
+        /*Get field-specific enums*/
+   switch(enum_type){
+      case(FrontalForcingsSubglacialDischargearmaEnum):
+         basinenum_type = FrontalForcingsBasinIdEnum;
+         varenum_type   = FrontalForcingsSubglacialDischargeEnum;
+         break;
+        }
+        
+        /*Evaluate the month index now and at (now-timestepjump)*/
+        this->parameters->FindParam(&yts,ConstantsYtsEnum);
+        this->parameters->FindParam(&time,TimeEnum);
+   this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); _assert_(dt<yts);
+        fracyear     = time/yts-floor(time/yts);
+        fracyearnext = (time+dt)/yts-floor((time+dt)/yts);
+        for(int i=1;i<12;i++){
+                if(fracyear>=monthsteps[i-1])     mindex     = i-1;
+                if(fracyearnext>=monthsteps[i-1]) mindexnext = i-1;
+        }
+        if(fracyear>=monthsteps[11])         mindex     = 11;
+        if(fracyearnext>=monthsteps[11])     mindexnext = 11;
+
+        /*Calculate fraction of the time step spent in each month*/
+        for(int i=0;i<12;i++){
+                if(mindex<i && mindexnext>i)                            fracdtinmonth[i] = 1.0/dt*yts/12.0;
+                else if(mindex<i && mindexnext<i && mindexnext<mindex)  fracdtinmonth[i] = 1.0/dt*yts/12.0;
+                else if(mindex>i && mindexnext<mindex && mindexnext>i)  fracdtinmonth[i] = 1.0/dt*yts/12.0;
+                else if(mindex>i && mindexnext<mindex && mindexnext==i) fracdtinmonth[i] = 1.0/dt*yts*(fracyearnext-monthsteps[i]);
+                else if(mindex==i && mindexnext==i)                     fracdtinmonth[i] = 1.0/dt*yts*(fracyearnext-fracyear); 
+                else if(mindex==i && mindexnext!=mindex)                fracdtinmonth[i] = 1.0/dt*yts*(1.0/12-(fracyear-monthsteps[i]));
+                else if(mindexnext==i && mindex!=mindexnext)            fracdtinmonth[i] = 1.0/dt*yts*(fracyearnext-monthsteps[i]);
+                else                                                      fracdtinmonth[i] = 0.0;
+        }
+
+        /*Get basin-specific parameters of the element*/
+   this->GetInputValue(&basinid,basinenum_type);
+        for(int i=0;i<12;i++) monthlyfrac_b[i]   = monthlyfrac[basinid*12+i];
+
+        /*Retrieve input*/
+        this->GetInputListOnVertices(varlistinput,varenum_type);
+
+        /*Calculate monthly rate for each month and weight-average it for application over dt*/
+        for(int v=0;v<numvertices;v++){
+                for(int i=0;i<12;i++){
+                        rateinmonth[v*12+i] = varlistinput[v]*monthlyfrac_b[i]*12;
+                        varlist[v]          = varlist[v]+fracdtinmonth[i]*rateinmonth[v*12+i];
+                }
+        }
+        /*Update input*/
+   this->AddInput(varenum_type,varlist,P1Enum);
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(fracdtinmonth);
+        xDelete<IssmDouble>(rateinmonth);
+        xDelete<IssmDouble>(monthlyfrac_b);
+        xDelete<IssmDouble>(monthlyrate_b);
+        xDelete<IssmDouble>(varlist);
+        xDelete<IssmDouble>(varlistinput);
 }/*}}}*/
 void       Element::MonthlyPiecewiseLinearEffectBasin(int nummonthbreaks,IssmDouble* monthlyintercepts,IssmDouble* monthlytrends,IssmDouble* monthlydatebreaks,int enum_type){/*{{{*/
@@ -3844 +3929 @@
       bed_input->GetInputValue(&bed,&gauss);
       qsg_input->GetInputValue(&qsg,&gauss);
-      TF_input->GetInputValue(&TF,&gauss);
+                TF_input->GetInputValue(&TF,&gauss);
 
       if(basin_icefront_area[basinid]==0.) meltrates[iv]=0.;

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

@@ -164 +164 @@
                 void               MigrateGroundingLine(IssmDouble* sheet_ungrounding);
                 void               MismipFloatingiceMeltingRate();
+                void               MonthlyFractionBasin(IssmDouble* monthlyfrac,int enum_type); 
                 void               MonthlyPiecewiseLinearEffectBasin(int nummonthbreaks,IssmDouble* monthlyintercepts,IssmDouble* monthlytrends,IssmDouble* monthlydatebreaks,int enum_type); 
                 void               BeckmannGoosseFloatingiceMeltingRate();

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

@@ -20 +20 @@
                 case FrontalForcingsRignotarmaEnum:
                         Thermalforcingarmax(femmodel);
+                        bool isdischargearma;
+                        femmodel->parameters->FindParam(&isdischargearma,FrontalForcingsIsDischargeARMAEnum);
+                        if(isdischargearma==true) Subglacialdischargearmax(femmodel);
                         /*Do not break here, call IcefrontAreax(),RignotMeltParameterizationx()*/
                 case FrontalForcingsRignotEnum:
@@ -105 +108 @@
    xDelete<IssmDouble>(monthtrends);
 }/*}}}*/
+void Subglacialdischargearmax(FemModel* femmodel){/*{{{*/
+
+        /*Get time parameters*/
+   IssmDouble time,dt,starttime,tstep_arma;
+   femmodel->parameters->FindParam(&time,TimeEnum);
+   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+   femmodel->parameters->FindParam(&tstep_arma,FrontalForcingsSdARMATimestepEnum);
+
+   /*Determine if this is a time step for the ARMA model*/
+   bool isstepforarma = false;
+
+   #ifndef _HAVE_AD_
+   if((fmod(time,tstep_arma)<fmod((time-dt),tstep_arma)) || (time<=starttime+dt) || tstep_arma==dt) isstepforarma = true;
+   #else
+   _error_("not implemented yet");
+   #endif
+
+   /*Load parameters*/
+        bool isstochastic;
+   bool isdischargestochastic = false;
+        int M,N,arorder,maorder,numbasins,numparams,numbreaks,my_rank;
+   femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
+   femmodel->parameters->FindParam(&numparams,FrontalForcingsSdNumberofParamsEnum);
+   femmodel->parameters->FindParam(&numbreaks,FrontalForcingsSdNumberofBreaksEnum);
+   femmodel->parameters->FindParam(&arorder,FrontalForcingsSdarOrderEnum);
+   femmodel->parameters->FindParam(&maorder,FrontalForcingsSdmaOrderEnum);
+   IssmDouble* datebreaks        = NULL;
+   IssmDouble* arlagcoefs        = NULL;
+   IssmDouble* malagcoefs        = NULL;
+   IssmDouble* monthlyfrac       = NULL;
+        IssmDouble* polyparams        = NULL;
+
+   femmodel->parameters->FindParam(&datebreaks,&M,&N,FrontalForcingsSddatebreaksEnum);            _assert_(M==numbasins); _assert_(N==max(numbreaks,1));        
+   femmodel->parameters->FindParam(&polyparams,&M,&N,FrontalForcingsSdpolyparamsEnum);            _assert_(M==numbasins); _assert_(N==(numbreaks+1)*numparams);        
+   femmodel->parameters->FindParam(&arlagcoefs,&M,&N,FrontalForcingsSdarlagcoefsEnum);            _assert_(M==numbasins); _assert_(N==arorder);
+   femmodel->parameters->FindParam(&malagcoefs,&M,&N,FrontalForcingsSdmalagcoefsEnum);            _assert_(M==numbasins); _assert_(N==maorder);
+   femmodel->parameters->FindParam(&monthlyfrac,&M,&N,FrontalForcingsSdMonthlyFracEnum);          _assert_(M==numbasins); _assert_(N==12); 
+
+        femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
+        if(isstochastic){
+                int  numstochasticfields;
+                int* stochasticfields;
+                femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
+                femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
+                for(int i=0;i<numstochasticfields;i++){
+                        if(stochasticfields[i]==FrontalForcingsSubglacialDischargearmaEnum) isdischargestochastic = true;
+                }
+                xDelete<int>(stochasticfields);
+        }
+
+   /*Loop over each element to compute Subglacial Discharge at vertices*/
+   for(Object* &object:femmodel->elements->objects){
+      Element* element = xDynamicCast<Element*>(object);
+                /*Compute ARMA*/
+      element->ArmaProcess(isstepforarma,arorder,maorder,numparams,numbreaks,tstep_arma,polyparams,arlagcoefs,malagcoefs,datebreaks,isdischargestochastic,FrontalForcingsSubglacialDischargearmaEnum);
+                /*Scale with monthly fractions*/
+                element->MonthlyFractionBasin(monthlyfrac,FrontalForcingsSubglacialDischargeEnum);
+        }
+
+   /*Cleanup*/
+   xDelete<IssmDouble>(arlagcoefs);
+   xDelete<IssmDouble>(malagcoefs);
+   xDelete<IssmDouble>(monthlyfrac);
+   xDelete<IssmDouble>(polyparams);
+   xDelete<IssmDouble>(datebreaks);
+}/*}}}*/
+

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

@@ -7 +7 @@
 /* local prototypes: */
 void FrontalForcingsx(FemModel* femmodel);
+void Subglacialdischargearmax(FemModel* femmodel);
 void Thermalforcingarmax(FemModel* femmodel);
 

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

@@ -76 +76 @@
 
                 /*Deal with the ARMA models*/
-                if(fields[j]==SMBarmaEnum || fields[j]==FrontalForcingsRignotarmaEnum || fields[j]==BasalforcingsDeepwaterMeltingRatearmaEnum){
+                if(fields[j]==SMBarmaEnum || fields[j]==FrontalForcingsRignotarmaEnum || fields[j]==BasalforcingsDeepwaterMeltingRatearmaEnum || fields[j]==FrontalForcingsSubglacialDischargearmaEnum){
                         switch(fields[j]){
                                 case SMBarmaEnum:
@@ -90 +90 @@
                                         noiseenum_type = BasalforcingsDeepwaterMeltingRateNoiseEnum;
                                         break;
+                                case FrontalForcingsSubglacialDischargearmaEnum:
+                                        dimenum_type   = FrontalForcingsBasinIdEnum;
+                                        noiseenum_type = SubglacialdischargeARMANoiseEnum;
+                                        break;  
                         }
                         for(Object* &object:femmodel->elements->objects){
@@ -106 +110 @@
                                 case FrontalForcingsRignotarmaEnum:
                                 case BasalforcingsDeepwaterMeltingRatearmaEnum:
+                                case FrontalForcingsSubglacialDischargearmaEnum:
                                         /*Already done above*/
                                         break;

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

@@ -221 +221 @@
 syn keyword cConstant FrontalForcingsARMAmonthtrendsEnum
 syn keyword cConstant FrontalForcingsARMApolyparamsEnum
+syn keyword cConstant FrontalForcingsIsDischargeARMAEnum
 syn keyword cConstant FrontalForcingsNumberofBasinsEnum
 syn keyword cConstant FrontalForcingsNumberofBreaksEnum
@@ -228 +229 @@
 syn keyword cConstant FrontalForcingsARMAarlagcoefsEnum
 syn keyword cConstant FrontalForcingsARMAmalagcoefsEnum
+syn keyword cConstant FrontalForcingsSdarlagcoefsEnum
+syn keyword cConstant FrontalForcingsSdARMATimestepEnum
+syn keyword cConstant FrontalForcingsSdarOrderEnum
+syn keyword cConstant FrontalForcingsSddatebreaksEnum
+syn keyword cConstant FrontalForcingsSdmalagcoefsEnum
+syn keyword cConstant FrontalForcingsSdmaOrderEnum
+syn keyword cConstant FrontalForcingsSdMonthlyFracEnum
+syn keyword cConstant FrontalForcingsSdNumberofBreaksEnum
+syn keyword cConstant FrontalForcingsSdNumberofParamsEnum
+syn keyword cConstant FrontalForcingsSdpolyparamsEnum
 syn keyword cConstant GrdModelEnum
 syn keyword cConstant GroundinglineFrictionInterpolationEnum
@@ -797 +808 @@
 syn keyword cConstant FrictionfEnum
 syn keyword cConstant FrontalForcingsBasinIdEnum
+syn keyword cConstant FrontalForcingsSubglacialDischargearmaEnum
 syn keyword cConstant FrontalForcingsSubglacialDischargeEnum
 syn keyword cConstant GeometryHydrostaticRatioEnum
@@ -1109 +1121 @@
 syn keyword cConstant StressTensoryzEnum
 syn keyword cConstant StressTensorzzEnum
+syn keyword cConstant SubglacialdischargeARMANoiseEnum
+syn keyword cConstant SubglacialdischargeValuesAutoregressionEnum
+syn keyword cConstant SubglacialdischargeValuesMovingaverageEnum
 syn keyword cConstant SurfaceAbsMisfitEnum
 syn keyword cConstant SurfaceAbsVelMisfitEnum

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

@@ -215 +215 @@
    FrontalForcingsARMAmonthtrendsEnum,
    FrontalForcingsARMApolyparamsEnum,
+        FrontalForcingsIsDischargeARMAEnum,
         FrontalForcingsNumberofBasinsEnum,
         FrontalForcingsNumberofBreaksEnum,
@@ -222 +223 @@
    FrontalForcingsARMAarlagcoefsEnum,
    FrontalForcingsARMAmalagcoefsEnum,
+   FrontalForcingsSdarlagcoefsEnum,
+        FrontalForcingsSdARMATimestepEnum,
+   FrontalForcingsSdarOrderEnum,
+   FrontalForcingsSddatebreaksEnum,
+   FrontalForcingsSdmalagcoefsEnum,
+   FrontalForcingsSdmaOrderEnum,
+   FrontalForcingsSdMonthlyFracEnum,
+        FrontalForcingsSdNumberofBreaksEnum,
+        FrontalForcingsSdNumberofParamsEnum,
+   FrontalForcingsSdpolyparamsEnum,
         GrdModelEnum,
         GroundinglineFrictionInterpolationEnum,
@@ -793 +804 @@
         FrictionfEnum,
         FrontalForcingsBasinIdEnum,
+        FrontalForcingsSubglacialDischargearmaEnum,
         FrontalForcingsSubglacialDischargeEnum,
         GeometryHydrostaticRatioEnum,
@@ -1106 +1118 @@
         StressTensoryzEnum,
         StressTensorzzEnum,
+        SubglacialdischargeARMANoiseEnum,
+        SubglacialdischargeValuesAutoregressionEnum,
+        SubglacialdischargeValuesMovingaverageEnum,
         SurfaceAbsMisfitEnum,
         SurfaceAbsVelMisfitEnum,

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

@@ -223 +223 @@
                 case FrontalForcingsARMAmonthtrendsEnum : return "FrontalForcingsARMAmonthtrends";
                 case FrontalForcingsARMApolyparamsEnum : return "FrontalForcingsARMApolyparams";
+                case FrontalForcingsIsDischargeARMAEnum : return "FrontalForcingsIsDischargeARMA";
                 case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins";
                 case FrontalForcingsNumberofBreaksEnum : return "FrontalForcingsNumberofBreaks";
@@ -230 +231 @@
                 case FrontalForcingsARMAarlagcoefsEnum : return "FrontalForcingsARMAarlagcoefs";
                 case FrontalForcingsARMAmalagcoefsEnum : return "FrontalForcingsARMAmalagcoefs";
+                case FrontalForcingsSdarlagcoefsEnum : return "FrontalForcingsSdarlagcoefs";
+                case FrontalForcingsSdARMATimestepEnum : return "FrontalForcingsSdARMATimestep";
+                case FrontalForcingsSdarOrderEnum : return "FrontalForcingsSdarOrder";
+                case FrontalForcingsSddatebreaksEnum : return "FrontalForcingsSddatebreaks";
+                case FrontalForcingsSdmalagcoefsEnum : return "FrontalForcingsSdmalagcoefs";
+                case FrontalForcingsSdmaOrderEnum : return "FrontalForcingsSdmaOrder";
+                case FrontalForcingsSdMonthlyFracEnum : return "FrontalForcingsSdMonthlyFrac";
+                case FrontalForcingsSdNumberofBreaksEnum : return "FrontalForcingsSdNumberofBreaks";
+                case FrontalForcingsSdNumberofParamsEnum : return "FrontalForcingsSdNumberofParams";
+                case FrontalForcingsSdpolyparamsEnum : return "FrontalForcingsSdpolyparams";
                 case GrdModelEnum : return "GrdModel";
                 case GroundinglineFrictionInterpolationEnum : return "GroundinglineFrictionInterpolation";
@@ -799 +810 @@
                 case FrictionfEnum : return "Frictionf";
                 case FrontalForcingsBasinIdEnum : return "FrontalForcingsBasinId";
+                case FrontalForcingsSubglacialDischargearmaEnum : return "FrontalForcingsSubglacialDischargearma";
                 case FrontalForcingsSubglacialDischargeEnum : return "FrontalForcingsSubglacialDischarge";
                 case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio";
@@ -1111 +1123 @@
                 case StressTensoryzEnum : return "StressTensoryz";
                 case StressTensorzzEnum : return "StressTensorzz";
+                case SubglacialdischargeARMANoiseEnum : return "SubglacialdischargeARMANoise";
+                case SubglacialdischargeValuesAutoregressionEnum : return "SubglacialdischargeValuesAutoregression";
+                case SubglacialdischargeValuesMovingaverageEnum : return "SubglacialdischargeValuesMovingaverage";
                 case SurfaceAbsMisfitEnum : return "SurfaceAbsMisfit";
                 case SurfaceAbsVelMisfitEnum : return "SurfaceAbsVelMisfit";

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

@@ -214 +214 @@
 syn keyword juliaConstC FrontalForcingsARMAmonthtrendsEnum
 syn keyword juliaConstC FrontalForcingsARMApolyparamsEnum
+syn keyword juliaConstC FrontalForcingsIsDischargeARMAEnum
 syn keyword juliaConstC FrontalForcingsNumberofBasinsEnum
 syn keyword juliaConstC FrontalForcingsNumberofBreaksEnum
@@ -221 +222 @@
 syn keyword juliaConstC FrontalForcingsARMAarlagcoefsEnum
 syn keyword juliaConstC FrontalForcingsARMAmalagcoefsEnum
+syn keyword juliaConstC FrontalForcingsSdarlagcoefsEnum
+syn keyword juliaConstC FrontalForcingsSdARMATimestepEnum
+syn keyword juliaConstC FrontalForcingsSdarOrderEnum
+syn keyword juliaConstC FrontalForcingsSddatebreaksEnum
+syn keyword juliaConstC FrontalForcingsSdmalagcoefsEnum
+syn keyword juliaConstC FrontalForcingsSdmaOrderEnum
+syn keyword juliaConstC FrontalForcingsSdMonthlyFracEnum
+syn keyword juliaConstC FrontalForcingsSdNumberofBreaksEnum
+syn keyword juliaConstC FrontalForcingsSdNumberofParamsEnum
+syn keyword juliaConstC FrontalForcingsSdpolyparamsEnum
 syn keyword juliaConstC GrdModelEnum
 syn keyword juliaConstC GroundinglineFrictionInterpolationEnum
@@ -790 +801 @@
 syn keyword juliaConstC FrictionfEnum
 syn keyword juliaConstC FrontalForcingsBasinIdEnum
+syn keyword juliaConstC FrontalForcingsSubglacialDischargearmaEnum
 syn keyword juliaConstC FrontalForcingsSubglacialDischargeEnum
 syn keyword juliaConstC GeometryHydrostaticRatioEnum
@@ -1102 +1114 @@
 syn keyword juliaConstC StressTensoryzEnum
 syn keyword juliaConstC StressTensorzzEnum
+syn keyword juliaConstC SubglacialdischargeARMANoiseEnum
+syn keyword juliaConstC SubglacialdischargeValuesAutoregressionEnum
+syn keyword juliaConstC SubglacialdischargeValuesMovingaverageEnum
 syn keyword juliaConstC SurfaceAbsMisfitEnum
 syn keyword juliaConstC SurfaceAbsVelMisfitEnum

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

@@ -226 +226 @@
               else if (strcmp(name,"FrontalForcingsARMAmonthtrends")==0) return FrontalForcingsARMAmonthtrendsEnum;
               else if (strcmp(name,"FrontalForcingsARMApolyparams")==0) return FrontalForcingsARMApolyparamsEnum;
+              else if (strcmp(name,"FrontalForcingsIsDischargeARMA")==0) return FrontalForcingsIsDischargeARMAEnum;
               else if (strcmp(name,"FrontalForcingsNumberofBasins")==0) return FrontalForcingsNumberofBasinsEnum;
               else if (strcmp(name,"FrontalForcingsNumberofBreaks")==0) return FrontalForcingsNumberofBreaksEnum;
@@ -233 +234 @@
               else if (strcmp(name,"FrontalForcingsARMAarlagcoefs")==0) return FrontalForcingsARMAarlagcoefsEnum;
               else if (strcmp(name,"FrontalForcingsARMAmalagcoefs")==0) return FrontalForcingsARMAmalagcoefsEnum;
+              else if (strcmp(name,"FrontalForcingsSdarlagcoefs")==0) return FrontalForcingsSdarlagcoefsEnum;
+              else if (strcmp(name,"FrontalForcingsSdARMATimestep")==0) return FrontalForcingsSdARMATimestepEnum;
+              else if (strcmp(name,"FrontalForcingsSdarOrder")==0) return FrontalForcingsSdarOrderEnum;
+              else if (strcmp(name,"FrontalForcingsSddatebreaks")==0) return FrontalForcingsSddatebreaksEnum;
+              else if (strcmp(name,"FrontalForcingsSdmalagcoefs")==0) return FrontalForcingsSdmalagcoefsEnum;
+              else if (strcmp(name,"FrontalForcingsSdmaOrder")==0) return FrontalForcingsSdmaOrderEnum;
+              else if (strcmp(name,"FrontalForcingsSdMonthlyFrac")==0) return FrontalForcingsSdMonthlyFracEnum;
+              else if (strcmp(name,"FrontalForcingsSdNumberofBreaks")==0) return FrontalForcingsSdNumberofBreaksEnum;
+              else if (strcmp(name,"FrontalForcingsSdNumberofParams")==0) return FrontalForcingsSdNumberofParamsEnum;
+              else if (strcmp(name,"FrontalForcingsSdpolyparams")==0) return FrontalForcingsSdpolyparamsEnum;
               else if (strcmp(name,"GrdModel")==0) return GrdModelEnum;
               else if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum;
@@ -249 +260 @@
               else if (strcmp(name,"HydrologyNumRequestedOutputs")==0) return HydrologyNumRequestedOutputsEnum;
               else if (strcmp(name,"HydrologyPressureMeltCoefficient")==0) return HydrologyPressureMeltCoefficientEnum;
-              else if (strcmp(name,"HydrologyRelaxation")==0) return HydrologyRelaxationEnum;
+         else stage=3;
+   }
+   if(stage==3){
+              if (strcmp(name,"HydrologyRelaxation")==0) return HydrologyRelaxationEnum;
               else if (strcmp(name,"HydrologyRequestedOutputs")==0) return HydrologyRequestedOutputsEnum;
               else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
@@ -260 +274 @@
               else if (strcmp(name,"HydrologydcEplMaxThickness")==0) return HydrologydcEplMaxThicknessEnum;
               else if (strcmp(name,"HydrologydcEplPoreWaterMass")==0) return HydrologydcEplPoreWaterMassEnum;
-         else stage=3;
-   }
-   if(stage==3){
-              if (strcmp(name,"HydrologydcEplThickComp")==0) return HydrologydcEplThickCompEnum;
+              else if (strcmp(name,"HydrologydcEplThickComp")==0) return HydrologydcEplThickCompEnum;
               else if (strcmp(name,"HydrologydcEplflipLock")==0) return HydrologydcEplflipLockEnum;
               else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum;
@@ -372 +383 @@
               else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
               else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum;
-              else if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum;
+         else stage=4;
+   }
+   if(stage==4){
+              if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum;
               else if (strcmp(name,"MeshNumberoflayers")==0) return MeshNumberoflayersEnum;
               else if (strcmp(name,"MeshNumberofvertices")==0) return MeshNumberofverticesEnum;
@@ -383 +397 @@
               else if (strcmp(name,"OceanGridNx")==0) return OceanGridNxEnum;
               else if (strcmp(name,"OceanGridNy")==0) return OceanGridNyEnum;
-         else stage=4;
-   }
-   if(stage==4){
-              if (strcmp(name,"OceanGridX")==0) return OceanGridXEnum;
+              else if (strcmp(name,"OceanGridX")==0) return OceanGridXEnum;
               else if (strcmp(name,"OceanGridY")==0) return OceanGridYEnum;
               else if (strcmp(name,"OutputBufferPointer")==0) return OutputBufferPointerEnum;
@@ -495 +506 @@
               else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum;
               else if (strcmp(name,"SolidearthSettingsMaxSHCoeff")==0) return SolidearthSettingsMaxSHCoeffEnum;
-              else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
+         else stage=5;
+   }
+   if(stage==5){
+              if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
               else if (strcmp(name,"SettingsNumResultsOnNodes")==0) return SettingsNumResultsOnNodesEnum;
               else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum;
@@ -506 +520 @@
               else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum;
               else if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum;
-         else stage=5;
-   }
-   if(stage==5){
-              if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum;
+              else if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum;
               else if (strcmp(name,"SmbAccugrad")==0) return SmbAccugradEnum;
               else if (strcmp(name,"SmbAccuref")==0) return SmbAccurefEnum;
@@ -618 +629 @@
               else if (strcmp(name,"TimesteppingCycleForcing")==0) return TimesteppingCycleForcingEnum;
               else if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum;
-              else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
+         else stage=6;
+   }
+   if(stage==6){
+              if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
               else if (strcmp(name,"TimesteppingTimeStepMax")==0) return TimesteppingTimeStepMaxEnum;
               else if (strcmp(name,"TimesteppingTimeStepMin")==0) return TimesteppingTimeStepMinEnum;
@@ -629 +643 @@
               else if (strcmp(name,"TransientAmrFrequency")==0) return TransientAmrFrequencyEnum;
               else if (strcmp(name,"TransientIsage")==0) return TransientIsageEnum;
-         else stage=6;
-   }
-   if(stage==6){
-              if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum;
+              else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum;
               else if (strcmp(name,"TransientIsdebris")==0) return TransientIsdebrisEnum;
               else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum;
@@ -741 +752 @@
               else if (strcmp(name,"DamageD")==0) return DamageDEnum;
               else if (strcmp(name,"DamageDOld")==0) return DamageDOldEnum;
-              else if (strcmp(name,"DamageDbar")==0) return DamageDbarEnum;
+         else stage=7;
+   }
+   if(stage==7){
+              if (strcmp(name,"DamageDbar")==0) return DamageDbarEnum;
               else if (strcmp(name,"DamageDbarOld")==0) return DamageDbarOldEnum;
               else if (strcmp(name,"DamageF")==0) return DamageFEnum;
@@ -752 +766 @@
               else if (strcmp(name,"DeltaDsl")==0) return DeltaDslEnum;
               else if (strcmp(name,"DslOld")==0) return DslOldEnum;
-         else stage=7;
-   }
-   if(stage==7){
-              if (strcmp(name,"Dsl")==0) return DslEnum;
+              else if (strcmp(name,"Dsl")==0) return DslEnum;
               else if (strcmp(name,"DeltaStr")==0) return DeltaStrEnum;
               else if (strcmp(name,"StrOld")==0) return StrOldEnum;
@@ -817 +828 @@
               else if (strcmp(name,"Frictionf")==0) return FrictionfEnum;
               else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum;
+              else if (strcmp(name,"FrontalForcingsSubglacialDischargearma")==0) return FrontalForcingsSubglacialDischargearmaEnum;
               else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum;
               else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
@@ -863 +875 @@
               else if (strcmp(name,"HydrologyMaskNodeActivation")==0) return HydrologyMaskNodeActivationEnum;
               else if (strcmp(name,"Ice")==0) return IceEnum;
-              else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
+         else stage=8;
+   }
+   if(stage==8){
+              if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
               else if (strcmp(name,"Input")==0) return InputEnum;
               else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
@@ -875 +890 @@
               else if (strcmp(name,"LevelsetObservation")==0) return LevelsetObservationEnum;
               else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum;
-         else stage=8;
-   }
-   if(stage==8){
-              if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
+              else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
               else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum;
               else if (strcmp(name,"MaskOceanLevelset")==0) return MaskOceanLevelsetEnum;
@@ -986 +998 @@
               else if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum;
               else if (strcmp(name,"SealevelchangeViscousSG")==0) return SealevelchangeViscousSGEnum;
-              else if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum;
+         else stage=9;
+   }
+   if(stage==9){
+              if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum;
               else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum;
               else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum;
@@ -998 +1013 @@
               else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum;
               else if (strcmp(name,"SigmaVM")==0) return SigmaVMEnum;
-         else stage=9;
-   }
-   if(stage==9){
-              if (strcmp(name,"SmbAccumulatedEC")==0) return SmbAccumulatedECEnum;
+              else if (strcmp(name,"SmbAccumulatedEC")==0) return SmbAccumulatedECEnum;
               else if (strcmp(name,"SmbAccumulatedMassBalance")==0) return SmbAccumulatedMassBalanceEnum;
               else if (strcmp(name,"SmbAccumulatedMelt")==0) return SmbAccumulatedMeltEnum;
@@ -1109 +1121 @@
               else if (strcmp(name,"SmbW")==0) return SmbWEnum;
               else if (strcmp(name,"SmbWAdd")==0) return SmbWAddEnum;
-              else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
+         else stage=10;
+   }
+   if(stage==10){
+              if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
               else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum;
               else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum;
@@ -1121 +1136 @@
               else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
               else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
-         else stage=10;
-   }
-   if(stage==10){
-              if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum;
+              else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum;
               else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum;
               else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum;
@@ -1138 +1150 @@
               else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
               else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
+              else if (strcmp(name,"SubglacialdischargeARMANoise")==0) return SubglacialdischargeARMANoiseEnum;
+              else if (strcmp(name,"SubglacialdischargeValuesAutoregression")==0) return SubglacialdischargeValuesAutoregressionEnum;
+              else if (strcmp(name,"SubglacialdischargeValuesMovingaverage")==0) return SubglacialdischargeValuesMovingaverageEnum;
               else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
               else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
@@ -1229 +1244 @@
               else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
               else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum;
-              else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
+         else stage=11;
+   }
+   if(stage==11){
+              if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
               else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
               else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
@@ -1244 +1262 @@
               else if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum;
               else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum;
-         else stage=11;
-   }
-   if(stage==11){
-              if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum;
+              else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum;
               else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
               else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
@@ -1352 +1367 @@
               else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum;
               else if (strcmp(name,"Channel")==0) return ChannelEnum;
-              else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
+         else stage=12;
+   }
+   if(stage==12){
+              if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
               else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum;
               else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum;
@@ -1367 +1385 @@
               else if (strcmp(name,"ControlInputValues")==0) return ControlInputValuesEnum;
               else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum;
-         else stage=12;
-   }
-   if(stage==12){
-              if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
+              else if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
               else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum;
               else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
@@ -1475 +1490 @@
               else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
               else if (strcmp(name,"IntParam")==0) return IntParamEnum;
-              else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
+         else stage=13;
+   }
+   if(stage==13){
+              if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
               else if (strcmp(name,"Inputs")==0) return InputsEnum;
               else if (strcmp(name,"Internal")==0) return InternalEnum;
@@ -1490 +1508 @@
               else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
               else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
-         else stage=13;
-   }
-   if(stage==13){
-              if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum;
+              else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum;
               else if (strcmp(name,"LinearFloatingMeltRatearma")==0) return LinearFloatingMeltRatearmaEnum;
               else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
@@ -1598 +1613 @@
               else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum;
               else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum;
-              else if (strcmp(name,"Regular")==0) return RegularEnum;
+         else stage=14;
+   }
+   if(stage==14){
+              if (strcmp(name,"Regular")==0) return RegularEnum;
               else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
               else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
@@ -1613 +1631 @@
               else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
               else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
-         else stage=14;
-   }
-   if(stage==14){
-              if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
+              else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
               else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;
               else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;

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

@@ -22 +22 @@
                 basin_id                 = NaN;
                 subglacial_discharge     = NaN;
+                isdischargearma          = 0;
+                sd_ar_order              = 0;
+                sd_ma_order              = 0;
+                sd_arma_timestep         = 0;
+                sd_arlag_coefs           = NaN;
+                sd_malag_coefs           = NaN;
+                sd_monthlyfrac           = NaN;
+                sd_num_breaks            = 0;
+                sd_num_params            = 0;
+                sd_polynomialparams      = NaN;
+                sd_datebreaks            = NaN;
         end
         methods
@@ -66 +77 @@
          md = checkfield(md,'fieldname','frontalforcings.num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0);
          md = checkfield(md,'fieldname','frontalforcings.basin_id','Inf',1,'>=',0,'<=',md.frontalforcings.num_basins,'size',[md.mesh.numberofelements 1]);
-         md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'NaN',1,'Inf',1,'timeseries',1);
          if(nbas>1 && nbrk>=1 && nprm>1)
             md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm);
@@ -127 +137 @@
                    end
 
+                        %%% Checking subglacial discharge %%%
+                        md = checkfield(md,'fieldname','frontalforcings.isdischargearma','values',[0 1]);
+                        if(~self.isdischargearma)
+                                md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'NaN',1,'Inf',1,'timeseries',1);
+                        else
+                                sdnbrk  = md.frontalforcings.sd_num_breaks; 
+                                sdnprm  = md.frontalforcings.sd_num_params;
+                                md = checkfield(md,'fieldname','frontalforcings.sd_ar_order','numel',1,'NaN',1,'Inf',1,'>=',0);
+                                md = checkfield(md,'fieldname','frontalforcings.sd_ma_order','numel',1,'NaN',1,'Inf',1,'>=',0);
+                md = checkfield(md,'fieldname','frontalforcings.sd_arma_timestep','numel',1,'NaN',1,'Inf',1,'>=',max(1,md.timestepping.time_step)); %ARMA time step cannot be finer than ISSM timestep and annual timestep
+                md = checkfield(md,'fieldname','frontalforcings.sd_arlag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.sd_ar_order]);
+                md = checkfield(md,'fieldname','frontalforcings.sd_malag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.sd_ma_order]);
+                md = checkfield(md,'fieldname','frontalforcings.sd_monthlyfrac','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,12]);
+                                if(any(abs(sum(self.sd_monthlyfrac,2)-1)>1e-3))
+                                        error('the 12 entries for each basin of md.frontalforcings.sd_monthlyfrac should add up to 1');
+                                end
+                 md = checkfield(md,'fieldname','frontalforcings.sd_num_params','numel',1,'NaN',1,'Inf',1,'>',0);
+                      md = checkfield(md,'fieldname','frontalforcings.sd_num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0);
+                 if(nbas>1 && sdnbrk>=1 && sdnprm>1)
+                    md = checkfield(md,'fieldname','frontalforcings.sd_polynomialparams','NaN',1,'Inf',1,'size',[nbas,sdnbrk+1,sdnprm],'numel',nbas*(sdnbrk+1)*sdnprm);
+                 elseif(nbas==1)
+                    md = checkfield(md,'fieldname','frontalforcings.sd_polynomialparams','NaN',1,'Inf',1,'size',[sdnprm,sdnbrk+1],'numel',nbas*(sdnbrk+1)*sdnprm);
+                 elseif(sdnbrk==0)
+                    md = checkfield(md,'fieldname','frontalforcings.sd_polynomialparams','NaN',1,'Inf',1,'size',[nbas,sdnprm],'numel',nbas*(sdnbrk+1)*sdnprm);
+                 elseif(sdnprm==1)
+                    md = checkfield(md,'fieldname','frontalforcings.sd_polynomialparams','NaN',1,'Inf',1,'size',[nbas,sdnbrk+1],'numel',nbas*(sdnbrk+1)*sdnprm);
+                 end
+                 if(sdnbrk>0)
+                    md = checkfield(md,'fieldname','frontalforcings.sd_datebreaks','NaN',1,'Inf',1,'size',[nbas,sdnbrk]);
+                 elseif(numel(md.frontalforcings.sd_datebreaks)==0 || all(isnan(md.frontalforcings.sd_datebreaks)))
+                    ;
+                 else
+                    error('md.frontalforcings.sd_num_breaks is 0 but md.frontalforcings.sd_datebreaks is not empty');
+                 end
+
+                        end
+
       end % }}}
                 function disp(self) % {{{
@@ -132 +179 @@
                         fielddisplay(self,'num_basins','number of different basins');
          fielddisplay(self,'basin_id','basin number assigned to each element [unitless]');
-         fielddisplay(self,'subglacial_discharge','sum of subglacial discharge for each basin [m/d]');
                         fielddisplay(self,'num_breaks','number of different breakpoints in the piecewise-polynomial (separating num_breaks+1 periods)');
                         fielddisplay(self,'num_params','number of different parameters in the piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)');
@@ -147 +193 @@
          fielddisplay(self,'monthlyvals_numbreaks','number of breakpoints in the piecewise-linear functions of monthly values');
          fielddisplay(self,'monthlyvals_datebreaks','dates at which the breakpoints in the piecewise-linear monthly values occur (1 row per basin)');
-                end % }}}
+         fielddisplay(self,'isdischargearma','whether an ARMA model is also used for the subglacial discharge (if 0: subglacial_discharge is used, if 1: sd_ parameters are used)');
+         fielddisplay(self,'subglacial_discharge','sum of subglacial discharge for each basin [m^3/d]');
+                        disp(sprintf('%51s  if isdischargearma is 1, sd_variables are used (sd arma model variable: sum of subglacial discharge for each basin [m^3/d])',' '));
+         fielddisplay(self,'sd_ar_order','order of the subglacial discharge autoregressive model [unitless]');
+         fielddisplay(self,'sd_ma_order','order of the subglacial discharge moving-average model [unitless]');
+         fielddisplay(self,'sd_arma_timestep','time resolution of the subglacial discharge autoregressive model [yr]');
+         fielddisplay(self,'sd_arlag_coefs','basin-specific vectors of AR lag coefficients for subglacial discharge [unitless]');
+         fielddisplay(self,'sd_malag_coefs','basin-specific vectors of MA lag coefficients for subglacial discharge [unitless]');
+         fielddisplay(self,'sd_monthlyfrac','basin-specific vectors of 12 values with fraction of the annual discharge occuring every month [unitless]');
+                        fielddisplay(self,'sd_num_params','number of different parameters in the subglacial discharge piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)');
+                        fielddisplay(self,'sd_num_breaks','number of different breakpoints in the subglacial discharge piecewise-polynomial (separating sd_num_breaks+1 periods)');
+         fielddisplay(self,'sd_datebreaks','dates at which the breakpoints in the piecewise polynomial occur (1 row per basin) [yr]');
+         fielddisplay(self,'sd_polynomialparams','coefficients for the sd_polynomial (const,trend,quadratic,etc.),dim1 for basins,dim2 for periods,dim3 for orders');
+                end % }}}
                 function marshall(self,prefix,md,fid) % {{{
                         yts=md.constants.yts;
@@ -228 +287 @@
                         end
 
+                        %%% Deal with the subglacial discharge polynomial %%%
+                        if(self.isdischargearma)
+                                sdnprm  = md.frontalforcings.sd_num_params;
+                                sdnper  = md.frontalforcings.sd_num_breaks+1;
+                                sdpolyparamsScaled   = md.frontalforcings.sd_polynomialparams;
+                 sdpolyparams2dScaled = zeros(nbas,sdnper*sdnprm);
+                 if(sdnprm>1)
+                    % Case 3D %
+                    if(nbas>1 && sdnper>1)
+                       for(ii=[1:sdnprm])
+                          sdpolyparamsScaled(:,:,ii) = sdpolyparamsScaled(:,:,ii)*((1/yts)^(ii-1));
+                       end
+                       % Fit in 2D array %
+                       for(ii=[1:sdnprm])
+                          jj = 1+(ii-1)*sdnper;
+                          sdpolyparams2dScaled(:,jj:jj+sdnper-1) = sdpolyparamsScaled(:,:,ii);
+                       end
+                    % Case 2D and higher-order params at increasing row index %
+                    elseif(nbas==1)
+                       for(ii=[1:sdnprm])
+                          sdpolyparamsScaled(ii,:) = sdpolyparamsScaled(ii,:)*((1/yts)^(ii-1));
+                       end
+                       % Fit in row array %
+                       for(ii=[1:sdnprm])
+                          jj = 1+(ii-1)*sdnper;
+                          sdpolyparams2dScaled(1,jj:jj+sdnper-1) = sdpolyparamsScaled(ii,:);
+                       end
+                    % Case 2D and higher-order params at incrasing column index %
+                    elseif(sdnper==1)
+                       for(ii=[1:sdnprm])
+                          sdpolyparamsScaled(:,ii) = sdpolyparamsScaled(:,ii)*((1/yts)^(ii-1));
+                       end
+                       % 2D array is already in correct format %
+                                                sdpolyparams2dScaled = sdpolyparamsScaled;
+                    end
+                 else
+                                        % 2D array is already in correct format and no need for scaling %
+               sdpolyparams2dScaled = sdpolyparamsScaled;
+                 end
+                 if(sdnper==1) %a single period (no break date)
+                    sd_dbreaks = zeros(nbas,1); %dummy
+                 else
+                    sd_dbreaks = md.frontalforcings.sd_datebreaks;
+                 end
+                        end
+
                         WriteData(fid,prefix,'name','md.frontalforcings.parameterization','data',3,'format','Integer');
                         WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','num_basins','format','Integer');
                         WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','num_breaks','format','Integer');
                         WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','num_params','format','Integer');
-                        WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','subglacial_discharge','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
          WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_order','format','Integer');
          WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ma_order','format','Integer');
@@ -245 +349 @@
          WriteData(fid,prefix,'data',interceptsM,'name','md.frontalforcings.monthlyvals_intercepts','format','DoubleMat');
          WriteData(fid,prefix,'data',trendsM,'name','md.frontalforcings.monthlyvals_trends','format','DoubleMat','scale',1/yts);
+                        WriteData(fid,prefix,'object',self,'fieldname','isdischargearma','format','Boolean');
+                        if(self.isdischargearma==0)
+                                WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','subglacial_discharge','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
+                        else
+                                WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_num_breaks','format','Integer');
+                                WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_num_params','format','Integer');
+                 WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_ar_order','format','Integer');
+                 WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_ma_order','format','Integer');
+                 WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_arma_timestep','format','Double','scale',yts);
+                                WriteData(fid,prefix,'data',sdpolyparams2dScaled,'name','md.frontalforcings.sd_polynomialparams','format','DoubleMat');
+                 WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_arlag_coefs','format','DoubleMat','name','md.frontalforcings.sd_arlag_coefs','yts',yts);
+                 WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_malag_coefs','format','DoubleMat','name','md.frontalforcings.sd_malag_coefs','yts',yts);
+                 WriteData(fid,prefix,'data',sd_dbreaks,'name','md.frontalforcings.sd_datebreaks','format','DoubleMat','scale',yts);
+                 WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_monthlyfrac','format','DoubleMat','name','md.frontalforcings.sd_monthlyfrac','yts',yts);
+                        end
                 end % }}}
         end

issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py

@@ -31 +31 @@
         self.basin_id = np.nan
         self.subglacial_discharge = np.nan
+        self.isdischargearma = 0
+        self.sd_ar_order = 0
+        self.sd_ma_order = 0
+        self.sd_arma_timestep = 0
+        self.sd_arlag_coefs = np.nan
+        self.sd_malag_coefs = np.nan
+        self.sd_monthlyfrac = np.nan
+        self.sd_num_breaks  = 0
+        self.sd_num_params  = 0
+        self.sd_polynomialparams = np.nan
+        self.sd_datebreaks = np.nan
 
         if len(args) == 0:
@@ -41 +52 @@
         s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]'))
         s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]'))
-        s += '{}\n'.format(fielddisplay(self, 'subglacial_discharge', 'sum of subglacial discharge for each basin [m/d]'))
         s += '{}\n'.format(fielddisplay(self, 'num_breaks', 'number of different breakpoints in the piecewise-polynomial (separating num_breaks+1 periods)'))
         s += '{}\n'.format(fielddisplay(self, 'num_params', 'number of different parameters in the piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)'))
@@ -51 +61 @@
         s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of AR lag coefficients [unitless]'))
         s += '{}\n'.format(fielddisplay(self, 'malag_coefs', 'basin-specific vectors of MA lag coefficients [unitless]'))
+        s += '{}\n'.format(fielddisplay(self, 'isdischargearma','whether an ARMA model is also used for the subglacial discharge (if 0: subglacial_discharge is used, if 1: sd_ parameters are used)'))
+        s += '{}\n'.format(fielddisplay(self, 'subglacial_discharge', 'sum of subglacial discharge for each basin [m/d]'))
+        s += '{}\n'.format(fielddisplay(self, 'sd_ar_order','order of the subglacial discharge autoregressive model [unitless]'))
+        s += '{}\n'.format(fielddisplay(self, 'sd_ma_order','order of the subglacial discharge moving-average model [unitless]'))
+        s += '{}\n'.format(fielddisplay(self, 'sd_arma_timestep','time resolution of the subglacial discharge autoregressive model [yr]'))
+        s += '{}\n'.format(fielddisplay(self, 'sd_arlag_coefs','basin-specific vectors of AR lag coefficients for subglacial discharge [unitless]'))
+        s += '{}\n'.format(fielddisplay(self, 'sd_malag_coefs','basin-specific vectors of MA lag coefficients for subglacial discharge [unitless]'))
+        s += '{}\n'.format(fielddisplay(self, 'sd_monthlyfrac','basin-specific vectors of 12 values with fraction of the annual discharge occuring every month [unitless]'))
+        s += '{}\n'.format(fielddisplay(self, 'sd_num_params','number of different parameters in the subglacial discharge piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)'))
+        s += '{}\n'.format(fielddisplay(self, 'sd_num_breaks','number of different breakpoints in the subglacial discharge piecewise-polynomial (separating sd_num_breaks+1 periods)'))
+        s += '{}\n'.format(fielddisplay(self, 'sd_datebreaks','dates at which the breakpoints in the piecewise polynomial occur (1 row per basin) [yr]'))
+        s += '{}\n'.format(fielddisplay(self, 'sd_polynomialparams','coefficients for the sd_polynomial (const,trend,quadratic,etc.),dim1 for basins,dim2 for periods,dim3 for orders'))
         return s
     #}}}
@@ -130 +152 @@
             raise RuntimeError('md.frontalforcings.monthlyvals_numbreaks is 0 but md.frontalforcings.monthlyvals_datebreaks is not empty')
 
+        ### Chacking subglacial discharge ###
+        md = checkfield(md, 'fieldname', 'frontalforcings.isdischargearma', 'values', [0, 1])
+        if(self.isdischargearma==0):
+            md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'NaN',1,'Inf',1,'timeseries',1)
+        else:
+            sdnbrk  = md.frontalforcings.sd_num_breaks
+            sdnprm  = md.frontalforcings.sd_num_params
+            md = checkfield(md,'fieldname','frontalforcings.sd_ar_order','numel',1,'NaN',1,'Inf',1,'>=',0)
+            md = checkfield(md,'fieldname','frontalforcings.sd_ma_order','numel',1,'NaN',1,'Inf',1,'>=',0)
+            md = checkfield(md,'fieldname','frontalforcings.sd_arma_timestep','numel',1,'NaN',1,'Inf',1,'>=',max(1,md.timestepping.time_step)) #ARMA time step cannot be finer than ISSM timestep and annual timestep
+            md = checkfield(md,'fieldname','frontalforcings.sd_arlag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.sd_ar_order])
+            md = checkfield(md,'fieldname','frontalforcings.sd_malag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.sd_ma_order])
+            md = checkfield(md,'fieldname','frontalforcings.sd_monthlyfrac','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,12])
+            if(np.any(abs(np.sum(self.sd_monthlyfrac,axis=1)-1)>1e-3)):
+                raise RuntimeError('the 12 entries for each basin of md.frontalforcings.sd_monthlyfrac should add up to 1')
+            md = checkfield(md,'fieldname','frontalforcings.sd_num_params','numel',1,'NaN',1,'Inf',1,'>',0)
+            md = checkfield(md,'fieldname','frontalforcings.sd_num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0)
+            if len(np.shape(self.sd_polynomialparams)) == 1:
+                self.sd_polynomialparams = np.array([[self.sd_polynomialparams]])
+            if(nbas>1 and sdnbrk>=1 and sdnprm>1):
+                md = checkfield(md,'fieldname','frontalforcings.sd_polynomialparams','NaN',1,'Inf',1,'size',[nbas,sdnbrk+1,sdnprm],'numel',nbas*(sdnbrk+1)*sdnprm)
+            elif(nbas==1):
+                md = checkfield(md,'fieldname','frontalforcings.sd_polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(sdnbrk+1)*sdnprm)
+            elif(sdnbrk==0):
+                md = checkfield(md,'fieldname','frontalforcings.sd_polynomialparams','NaN',1,'Inf',1,'size',[nbas,sdnprm],'numel',nbas*(sdnbrk+1)*sdnprm)
+            elif(sdnprm==1):
+                md = checkfield(md,'fieldname','frontalforcings.sd_polynomialparams','NaN',1,'Inf',1,'size',[nbas,sdnbrk],'numel',nbas*(sdnbrk+1)*sdnprm)
+            if(sdnbrk>0):
+                md = checkfield(md, 'fieldname', 'frontalforcings.sd_datebreaks', 'NaN', 1, 'Inf', 1, 'size', [nbas,sdnbrk])
+            elif(np.size(md.frontalforcings.sd_datebreaks)==0 or np.all(np.isnan(md.frontalforcings.sd_datebreaks))):
+                pass
+            else:
+                raise RuntimeError('md.frontalforcings.sd_num_breaks is 0 but md.frontalforcings.sd_datebreaks is not empty')
+
         return md
     # }}}
@@ -201 +257 @@
             dMbreaks = np.copy(md.frontalforcings.monthlyvals_datebreaks)
 
+        ### Deal with the subglacial discharge polynomial ###
+        if(self.isdischargearma):
+            sdnprm  = md.frontalforcings.sd_num_params
+            sdnper  = md.frontalforcings.sd_num_breaks+1
+            sdpolyparamsScaled   = np.copy(md.frontalforcings.sd_polynomialparams)
+            sdpolyparams2dScaled = np.zeros((nbas,sdnper*sdnprm))
+            if(sdnprm>1):
+                # Case 3D #
+                if(nbas>1 and sdnper>1):
+                    for ii in range(sdnprm):
+                        sdpolyparamsScaled[:,:,ii] = sdpolyparamsScaled[:,:,ii]*(1/yts)**ii
+                    # Fit in 2D array #
+                    for ii in range(sdnprm):
+                        sdpolyparams2dScaled[:,ii*sdnper:(ii+1)*sdnper] = 1*sdpolyparamsScaled[:,:,ii]
+                # Case 2D and higher-order params at increasing row index #
+                elif(nbas==1):
+                    for ii in range(sdnprm):
+                        sdpolyparamsScaled[ii,:] = sdpolyparamsScaled[ii,:]*(1/yts)**ii
+                    # Fit in row array #
+                    for ii in range(nprm):
+                        sdpolyparams2dScaled[0,ii*sdnper:(ii+1)*sdnper] = 1*sdpolyparamsScaled[ii,:]
+                # Case 2D and higher-order params at incrasing column index #
+                elif(sdnper==1):
+                    for ii in range(sdnprm):
+                        sdpolyparamsScaled[:,ii] = sdpolyparamsScaled[:,ii]*(1/yts)**ii
+                    # 2D array is already in correct format #
+                    sdpolyparams2dScaled = np.copy(sdpolyparamsScaled)
+            else:
+                # 2D array is already in correct format and no need for scaling #
+                sdpolyparams2dScaled = np.copy(sdpolyparamsScaled)
+            if(sdnper==1):
+                sd_dbreaks = np.zeros((nbas,1))
+            else:
+                sd_dbreaks = np.copy(md.frontalforcings.sd_datebreaks)
+
+
+
+
         WriteData(fid, prefix, 'name', 'md.frontalforcings.parameterization', 'data', 3, 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'num_basins', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'num_breaks', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'num_params', 'format', 'Integer')
-        WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'subglacial_discharge', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ar_order', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ma_order', 'format', 'Integer')
@@ -218 +311 @@
         WriteData(fid,prefix,'data',interceptsM,'name','md.frontalforcings.monthlyvals_intercepts','format','DoubleMat')
         WriteData(fid,prefix,'data',trendsM,'name','md.frontalforcings.monthlyvals_trends','format','DoubleMat','scale',1/yts)
+        WriteData(fid,prefix,'object',self,'fieldname','isdischargearma','format','Boolean')
+        if(self.isdischargearma==0):
+            WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'subglacial_discharge', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
+        else:
+            WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_num_breaks','format','Integer')
+            WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_num_params','format','Integer')
+            WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_ar_order','format','Integer')
+            WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_ma_order','format','Integer')
+            WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_arma_timestep','format','Double','scale',yts)
+            WriteData(fid,prefix,'data',sdpolyparams2dScaled,'name','md.frontalforcings.sd_polynomialparams','format','DoubleMat')
+            WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_arlag_coefs','format','DoubleMat','name','md.frontalforcings.sd_arlag_coefs','yts',yts)
+            WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_malag_coefs','format','DoubleMat','name','md.frontalforcings.sd_malag_coefs','yts',yts)
+            WriteData(fid,prefix,'data',sd_dbreaks,'name','md.frontalforcings.sd_datebreaks','format','DoubleMat','scale',yts)
+            WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','sd_monthlyfrac','format','DoubleMat','name','md.frontalforcings.sd_monthlyfrac','yts',yts)
     # }}}
 

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

@@ -120 +120 @@
                         indSMBarma = -1; %about to check for index of SMBarma
                         indTFarma  = -1; %about to check for index of FrontalForcingsRignotarma
+                        indSdarma  = -1; %about to check for index of SubglacialDischargearma
                         indBDWarma  = -1; %about to check for index of BasalforcingsDeepwaterMeltingRatearma
                         if any(contains(self.fields,'SMBarma'))
@@ -135 +136 @@
                                 end
                         end
+                        if any(contains(self.fields,'FrontalForcingsSubglacialDischargearma'))
+                                indSdarma       = find(contains(self.fields,'FrontalForcingsSubglacialDischargearma')); %index of Sdarma, now check for consistency with other arma timesteps 
+                                dimensions(indSdarma) = md.frontalforcings.num_basins;
+                                if(md.frontalforcings.sd_arma_timestep<self.stochastictimestep)
+                                        error('FrontalForcingsSubglacialDischargearma cannot have a timestep shorter than stochastictimestep');
+                                end
+                        end
                         if any(contains(self.fields,'BasalforcingsDeepwaterMeltingRatearma'))
                                 indBDWarma      = find(contains(self.fields,'BasalforcingsDeepwaterMeltingRatearma')); %index of BDWarma, now check for consistency with other arma timesteps 
@@ -165 +173 @@
                                         if any(crossentries~=0)
                                                 error('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance');
+                                        end
+                                end
+                        end
+                        if(indSMBarma~=-1 && indSdarma~=-1) %both ARMA models are used: check ARMA time step consistency
+                                if(md.smb.arma_timestep~=md.frontalforcings.sd_arma_timestep)
+                                        crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma))),1,[]);
+                                        if any(crossentries~=0)
+                                                error('SMBarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance');
+                                        end
+                                end
+                        end
+                        if(indSdarma~=-1 && indBDWarma~=-1) %both ARMA models are used: check ARMA time step consistency
+                                if(md.frontalforcings.sd_arma_timestep~=md.basalforcings.arma_timestep)
+                                        crossentries = reshape(self.covariance(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]);
+                                        if any(crossentries~=0)
+                                                error('FrontalForcingsSubglacialDischargearma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance');
+                                        end
+                                end
+                        end
+                        if(indSdarma~=-1 && indTFarma~=-1) %both ARMA models are used: check ARMA time step consistency
+                                if(md.frontalforcings.sd_arma_timestep~=md.frontalforcings.arma_timestep)
+                                        crossentries = reshape(self.covariance(1+sum(dimensions(1:indSdarma-1)):sum(dimensions(1:indSdarma)),1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma))),1,[]);
+                                        if any(crossentries~=0)
+                                                error('FrontalForcingsRignotarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance');
                                         end
                                 end
@@ -216 +248 @@
                                         end
                                         if(strcmp(field,'FrontalForcingsRignotarma'))
+                                                dimensions(ind) = md.frontalforcings.num_basins;
+                                        end
+                                        if(strcmp(field,'FrontalForcingsSubglacialDischargearma'))
                                                 dimensions(ind) = md.frontalforcings.num_basins;
                                         end
@@ -274 +309 @@
                 'FrictionSchoofWaterPressure',...
                 'FrontalForcingsRignotarma',...
+                'FrontalForcingsSubglacialDischargearma',...
                 'SMBarma',...
                 'SMBforcing'
@@ -286 +322 @@
                 'frictionschoof',...
                 'frontalforcingsrignotarma',...
+                'frontalforcingsrignotarma',...
                 'SMBarma',...
                 'SMBforcing'

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

@@ -36 +36 @@
         s += '{}\n'.format(fielddisplay(self, 'randomflag', 'whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'))
         s += 'Available fields:\n'
+        s += '   BasalforcingsDeepwaterMeltingRatearma\n'
         s += '   BasalforcingsSpatialDeepwaterMeltingRate\n'
         s += '   DefaultCalving\n'
         s += '   FloatingMeltRate\n'
         s += '   FrictionWaterPressure\n'
+        s += '   FrictionCoulombWaterPressure\n'
+        s += '   FrictionSchoofWaterPressure\n'
         s += '   FrontalForcingsRignotarma (thermal forcing)\n'
+        s += '   FrontalForcingsSubglacialDischargearma\n'
         s += '   SMBarma\n'
         s += '   SMBforcing\n'
@@ -122 +126 @@
         indSMBarma   = -1  # About to check for index of SMBarma
         indTFarma    = -1  # About to check for index of FrontalForcingsRignotarma
+        indSdarma    = -1  # About to check for index of FrontalForcingsSubglacialDischargearma
         indBDWarma   = -1  # About to check for index of BasalforcingsDeepwaterMeltingRatearma
         if ('SMBarma' in self.fields):
@@ -133 +138 @@
             if(md.frontalforcings.arma_timestep<self.stochastictimestep):
                 raise TypeError('FrontalForcingsRignotarma cannot have a timestep shorter than stochastictimestep')
+        if ('FrontalForcingsSubglacialDischargearma' in self.fields):
+            indSdarma = self.fields.index('FrontalForcingsSubglacialDischargearma')  # Index of Sdarma, now check for consistency with other timesteps
+            dimensions[indSdarma] = md.frontalforcings.num_basins
+            if(md.frontalforcings.sd_arma_timestep<self.stochastictimestep):
+                raise TypeError('FrontalForcingsSubglacialDischargearma cannot have a timestep shorter than stochastictimestep')
         if ('BasalforcingsDeepwaterMeltingRatearma' in self.fields):
             indBDWarma = self.fields.index('BasalforcingsDeepwaterMeltingRatearma')  # Index of BDWarma, now check for consistency with other timesteps
@@ -152 +162 @@
             if((md.frontalforcings.arma_timestep != md.basalforcings.arma_timestep) and np.any(covsum != 0)):
                 raise IOError('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
+        if (indSMBarma != -1 and indSdarma != -1):  # Both ARMA models are used: check ARMA time step consistency
+            covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int)]
+            if((md.smb.arma_timestep != md.frontalforcings.sd_arma_timestep) and np.any(covsum != 0)):
+                raise IOError('SMBarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance')
+        if (indSdarma != -1 and indBDWarma != -1):  # Both ARMA models are used: check ARMA time step consistency
+            covsum = self.covariance[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)]
+            if((md.frontalforcings.sd_arma_timestep != md.basalforcings.arma_timestep) and np.any(covsum != 0)):
+                raise IOError('FrontalForcingsSubglacialDischargearma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
+        if (indSdarma != -1 and indTFarma != -1):  # Both ARMA models are used: check ARMA time step consistency
+            covsum = self.covariance[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int)]
+            if((md.frontalforcings.sd_arma_timestep != md.frontalforcings.arma_timestep) and np.any(covsum != 0)):
+                raise IOError('FrontalForcingsSubglacialDischargearma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance')
+
 
         md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1])
@@ -187 +210 @@
                     dimensions[ind] = md.smb.num_basins
                 if (field == 'FrontalForcingsRignotarma'):
+                    dimensions[ind] = md.frontalforcings.num_basins
+                if (field == 'FrontalForcingsSubglacialDischargearma'):
                     dimensions[ind] = md.frontalforcings.num_basins
                 if (field == 'BasalforcingsDeepwaterMeltingRatearma'):
@@ -237 +262 @@
                      'FrictionSchoofWaterPressure': 'frictionschoof',
                      'FrontalForcingsRignotarma': 'frontalforcingsrignotarma',
+                     'FrontalForcingsSubglacialDischargearma': 'frontalforcingsrignotarma',
                      'SMBarma': 'SMBarma',
                      'SMBforcing': 'SMBforcing'}

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

@@ -67 +67 @@
 trendsM          = cat(3,trendsMp0,trendsMp1);
 datebreaksM      = [1;1];
+% Subglacial discharge parameters %
+isdischargearma            = 1;
+sd_ar_order                = 1;
+sd_ma_order                = 1;
+sd_num_breaks              = 1;
+sd_num_params              = 2;
+sd_arma_timestep           = 1;
+sd_arlag_coefs             = [0.95;0.95];
+sd_malag_coefs             = [0;0];
+sd_datebreaks              = [1;1];
+sd_monthlyfrac             = [0,0,0,0,0,0,0.5,0.5,0,0,0,0;
+                              0,0,0,0,0,0,0.5,0.5,0,0,0,0];
+sd_const                   = [50000,70000.0;8000,10000.0];
+sd_trend                   = [0,10000;0,0];
+sd_polyparam               = cat(3,sd_const,sd_trend);
 
 md.frontalforcings.num_basins              = nb_tf;
@@ -84 +99 @@
 md.frontalforcings.monthlyvals_trends      = trendsM;
 md.frontalforcings.monthlyvals_datebreaks  = datebreaksM;
+md.frontalforcings.isdischargearma         = isdischargearma;
+if(isdischargearma==0)
+        md.frontalforcings.subglacial_discharge    = 0.01*ones(md.mesh.numberofvertices,1);
+else
+    md.frontalforcings.sd_num_breaks         = sd_num_breaks;
+    md.frontalforcings.sd_num_params         = sd_num_params;
+    md.frontalforcings.sd_ar_order           = sd_ar_order;
+    md.frontalforcings.sd_ma_order           = sd_ma_order;
+    md.frontalforcings.sd_arma_timestep      = sd_arma_timestep;
+    md.frontalforcings.sd_arlag_coefs        = sd_arlag_coefs;
+    md.frontalforcings.sd_malag_coefs        = sd_malag_coefs;
+    md.frontalforcings.sd_datebreaks         = sd_datebreaks;
+    md.frontalforcings.sd_monthlyfrac        = sd_monthlyfrac;
+    md.frontalforcings.sd_polynomialparams   = sd_polyparam;
+end
 % Floating Ice Melt parameters
 md.basalforcings.floatingice_melting_rate = 0.1*ones(md.mesh.numberofvertices,1);
+
 
 % Covariance matrix
@@ -92 +123 @@
 covclv(1,1) = 1/10*covclv(1,1);
 covflmlt    = 0.05*eye(nb_flmlt);
-covglob     = blkdiag(covtf,covclv,covflmlt);
+covsd       = 1e3*eye(nb_tf);
+covglob     = blkdiag(covtf,covclv,covflmlt,covsd);
 
 %Stochastic forcing
 md.stochasticforcing.isstochasticforcing = 1;
-md.stochasticforcing.fields              = [{'FrontalForcingsRignotarma'},{'DefaultCalving'},{'FloatingMeltRate'}];
+md.stochasticforcing.fields              = [{'FrontalForcingsRignotarma'},{'DefaultCalving'},{'FloatingMeltRate'},{'FrontalForcingsSubglacialDischargearma'}];
 md.stochasticforcing.defaultdimension    = 2;
 md.stochasticforcing.default_id          = idb_df;

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

@@ -72 +72 @@
 trendsM          = np.transpose(np.stack((trendsMp0,trendsMp1)),(1,2,0)) 
 datebreaksM      = np.array([[1],[1]]) 
+# Subglacial discharge params #
+isdischargearma            = 1
+sd_ar_order                = 1
+sd_ma_order                = 1
+sd_num_breaks              = 1
+sd_num_params              = 2
+sd_arma_timestep           = 1
+sd_arlag_coefs             = np.array([[0.95],[0.95]])
+sd_malag_coefs             = np.array([[0.0],[0.0]])
+sd_datebreaks              = np.array([[1.0],[1.0]])
+sd_monthlyfrac             = np.array([[0,0,0,0,0,0,0.5,0.5,0,0,0,0],[0,0,0,0,0,0,0.5,0.5,0,0,0,0]])
+sd_const                   = np.array([[50000,70000],[8000,10000.0]])
+sd_trend                   = np.array([[0.0,10000],[0,0]])
+sd_polyparam               = np.transpose(np.stack((sd_const,sd_trend)),(1,2,0))
+
 
 md.frontalforcings.num_basins = nb_tf
@@ -89 +104 @@
 md.frontalforcings.monthlyvals_trends      = trendsM
 md.frontalforcings.monthlyvals_datebreaks  = datebreaksM
+md.frontalforcings.isdischargearma         = isdischargearma
+if(isdischargearma==0):
+   md.frontalforcings.subglacial_discharge    = 0.01*ones(md.mesh.numberofvertices,1)
+else:
+    md.frontalforcings.sd_num_breaks         = sd_num_breaks
+    md.frontalforcings.sd_num_params         = sd_num_params
+    md.frontalforcings.sd_ar_order           = sd_ar_order
+    md.frontalforcings.sd_ma_order           = sd_ma_order
+    md.frontalforcings.sd_arma_timestep      = sd_arma_timestep
+    md.frontalforcings.sd_arlag_coefs        = sd_arlag_coefs
+    md.frontalforcings.sd_malag_coefs        = sd_malag_coefs
+    md.frontalforcings.sd_datebreaks         = sd_datebreaks
+    md.frontalforcings.sd_monthlyfrac        = sd_monthlyfrac
+    md.frontalforcings.sd_polynomialparams   = sd_polyparam
 #Floating Ice Melt parameters
 md.basalforcings.floatingice_melting_rate = 0.1 * np.ones((md.mesh.numberofvertices,))
@@ -103 +132 @@
 
 #Hard-coding covariance matrix because python is complaining
-covglob = np.array([[1e-4, 0., 0., 0., 0., 0.],
-                    [0., 1e-4, 0., 0., 0., 0.],
-                    [0., 0., 1e-2, 0., 0., 0.],
-                    [0., 0., 0., 1e-1, 0., 0.],
-                    [0., 0., 0., 0., 0.05, 0.],
-                    [0., 0., 0., 0., 0., 0.05]])
+covglob = np.array([[1e-4, 0., 0., 0., 0., 0.,0.,0.],
+                    [0., 1e-4, 0., 0., 0., 0.,0.,0.],
+                    [0., 0., 1e-2, 0., 0., 0.,0.,0.],
+                    [0., 0., 0., 1e-1, 0., 0.,0.,0.],
+                    [0., 0., 0., 0., 0.05, 0.,0.,0.],
+                    [0., 0., 0., 0., 0., 0.05,0.,0.],
+                    [0., 0., 0., 0., 0., 0., 1e3,0.],
+                    [0., 0., 0., 0., 0., 0.,0., 1e3]])
 #testchol = np.linalg.cholesky(covglob)
 #print(testchol)
@@ -114 +145 @@
 # Stochastic forcing
 md.stochasticforcing.isstochasticforcing = 1
-md.stochasticforcing.fields = ['FrontalForcingsRignotarma', 'DefaultCalving', 'FloatingMeltRate']
+md.stochasticforcing.fields = ['FrontalForcingsRignotarma', 'DefaultCalving', 'FloatingMeltRate','FrontalForcingsSubglacialDischargearma']
 md.stochasticforcing.defaultdimension = 2
 md.stochasticforcing.default_id = idb_df