Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 27353)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 27354)
@@ -153,9 +153,13 @@
          /*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,5 +249,6 @@
 			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,5 +278,31 @@
          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));
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 27353)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 27354)
@@ -99,5 +99,12 @@
          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,5 +133,5 @@
 			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,4 +142,5 @@
 		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,5 +182,10 @@
 			/*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,4 +2729,76 @@
 	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,5 +3929,5 @@
       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.;
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 27353)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 27354)
@@ -164,4 +164,5 @@
 		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();
Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp	(revision 27353)
+++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp	(revision 27354)
@@ -20,4 +20,7 @@
 		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,2 +108,70 @@
    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);
+}/*}}}*/
+
Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h	(revision 27353)
+++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h	(revision 27354)
@@ -7,4 +7,5 @@
 /* local prototypes: */
 void FrontalForcingsx(FemModel* femmodel);
+void Subglacialdischargearmax(FemModel* femmodel);
 void Thermalforcingarmax(FemModel* femmodel);
 
Index: /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 27353)
+++ /issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp	(revision 27354)
@@ -76,5 +76,5 @@
 
 		/*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,4 +90,8 @@
 					noiseenum_type = BasalforcingsDeepwaterMeltingRateNoiseEnum;
 					break;
+				case FrontalForcingsSubglacialDischargearmaEnum:
+					dimenum_type   = FrontalForcingsBasinIdEnum;
+					noiseenum_type = SubglacialdischargeARMANoiseEnum;
+					break;	
 			}
 			for(Object* &object:femmodel->elements->objects){
@@ -106,4 +110,5 @@
 				case FrontalForcingsRignotarmaEnum:
 				case BasalforcingsDeepwaterMeltingRatearmaEnum:
+				case FrontalForcingsSubglacialDischargearmaEnum:
 					/*Already done above*/
 					break;
Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 27353)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 27354)
@@ -221,4 +221,5 @@
 syn keyword cConstant FrontalForcingsARMAmonthtrendsEnum
 syn keyword cConstant FrontalForcingsARMApolyparamsEnum
+syn keyword cConstant FrontalForcingsIsDischargeARMAEnum
 syn keyword cConstant FrontalForcingsNumberofBasinsEnum
 syn keyword cConstant FrontalForcingsNumberofBreaksEnum
@@ -228,4 +229,14 @@
 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,4 +808,5 @@
 syn keyword cConstant FrictionfEnum
 syn keyword cConstant FrontalForcingsBasinIdEnum
+syn keyword cConstant FrontalForcingsSubglacialDischargearmaEnum
 syn keyword cConstant FrontalForcingsSubglacialDischargeEnum
 syn keyword cConstant GeometryHydrostaticRatioEnum
@@ -1109,4 +1121,7 @@
 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
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27353)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 27354)
@@ -215,4 +215,5 @@
    FrontalForcingsARMAmonthtrendsEnum,
    FrontalForcingsARMApolyparamsEnum,
+	FrontalForcingsIsDischargeARMAEnum,
 	FrontalForcingsNumberofBasinsEnum,
 	FrontalForcingsNumberofBreaksEnum,
@@ -222,4 +223,14 @@
    FrontalForcingsARMAarlagcoefsEnum,
    FrontalForcingsARMAmalagcoefsEnum,
+   FrontalForcingsSdarlagcoefsEnum,
+	FrontalForcingsSdARMATimestepEnum,
+   FrontalForcingsSdarOrderEnum,
+   FrontalForcingsSddatebreaksEnum,
+   FrontalForcingsSdmalagcoefsEnum,
+   FrontalForcingsSdmaOrderEnum,
+   FrontalForcingsSdMonthlyFracEnum,
+	FrontalForcingsSdNumberofBreaksEnum,
+	FrontalForcingsSdNumberofParamsEnum,
+   FrontalForcingsSdpolyparamsEnum,
 	GrdModelEnum,
 	GroundinglineFrictionInterpolationEnum,
@@ -793,4 +804,5 @@
 	FrictionfEnum,
 	FrontalForcingsBasinIdEnum,
+	FrontalForcingsSubglacialDischargearmaEnum,
 	FrontalForcingsSubglacialDischargeEnum,
 	GeometryHydrostaticRatioEnum,
@@ -1106,4 +1118,7 @@
 	StressTensoryzEnum,
 	StressTensorzzEnum,
+	SubglacialdischargeARMANoiseEnum,
+	SubglacialdischargeValuesAutoregressionEnum,
+	SubglacialdischargeValuesMovingaverageEnum,
 	SurfaceAbsMisfitEnum,
 	SurfaceAbsVelMisfitEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 27353)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 27354)
@@ -223,4 +223,5 @@
 		case FrontalForcingsARMAmonthtrendsEnum : return "FrontalForcingsARMAmonthtrends";
 		case FrontalForcingsARMApolyparamsEnum : return "FrontalForcingsARMApolyparams";
+		case FrontalForcingsIsDischargeARMAEnum : return "FrontalForcingsIsDischargeARMA";
 		case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins";
 		case FrontalForcingsNumberofBreaksEnum : return "FrontalForcingsNumberofBreaks";
@@ -230,4 +231,14 @@
 		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,4 +810,5 @@
 		case FrictionfEnum : return "Frictionf";
 		case FrontalForcingsBasinIdEnum : return "FrontalForcingsBasinId";
+		case FrontalForcingsSubglacialDischargearmaEnum : return "FrontalForcingsSubglacialDischargearma";
 		case FrontalForcingsSubglacialDischargeEnum : return "FrontalForcingsSubglacialDischarge";
 		case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio";
@@ -1111,4 +1123,7 @@
 		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";
Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 27353)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim	(revision 27354)
@@ -214,4 +214,5 @@
 syn keyword juliaConstC FrontalForcingsARMAmonthtrendsEnum
 syn keyword juliaConstC FrontalForcingsARMApolyparamsEnum
+syn keyword juliaConstC FrontalForcingsIsDischargeARMAEnum
 syn keyword juliaConstC FrontalForcingsNumberofBasinsEnum
 syn keyword juliaConstC FrontalForcingsNumberofBreaksEnum
@@ -221,4 +222,14 @@
 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,4 +801,5 @@
 syn keyword juliaConstC FrictionfEnum
 syn keyword juliaConstC FrontalForcingsBasinIdEnum
+syn keyword juliaConstC FrontalForcingsSubglacialDischargearmaEnum
 syn keyword juliaConstC FrontalForcingsSubglacialDischargeEnum
 syn keyword juliaConstC GeometryHydrostaticRatioEnum
@@ -1102,4 +1114,7 @@
 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
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 27353)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 27354)
@@ -226,4 +226,5 @@
 	      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,4 +234,14 @@
 	      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,5 +260,8 @@
 	      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,8 +274,5 @@
 	      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,5 +383,8 @@
 	      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,8 +397,5 @@
 	      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,5 +506,8 @@
 	      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,8 +520,5 @@
 	      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,5 +629,8 @@
 	      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,8 +643,5 @@
 	      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,5 +752,8 @@
 	      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,8 +766,5 @@
 	      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,4 +828,5 @@
 	      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,5 +875,8 @@
 	      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,8 +890,5 @@
 	      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,5 +998,8 @@
 	      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,8 +1013,5 @@
 	      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,5 +1121,8 @@
 	      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,8 +1136,5 @@
 	      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,4 +1150,7 @@
 	      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,5 +1244,8 @@
 	      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,8 +1262,5 @@
 	      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,5 +1367,8 @@
 	      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,8 +1385,5 @@
 	      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,5 +1490,8 @@
 	      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,8 +1508,5 @@
 	      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,5 +1613,8 @@
 	      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,8 +1631,5 @@
 	      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;
Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m	(revision 27353)
+++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m	(revision 27354)
@@ -22,4 +22,15 @@
 		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,5 +77,4 @@
          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,4 +137,41 @@
 		   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,5 +179,4 @@
 			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,5 +193,18 @@
          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,9 +287,54 @@
 			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,4 +349,19 @@
          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
Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py	(revision 27353)
+++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py	(revision 27354)
@@ -31,4 +31,15 @@
         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,5 +52,4 @@
         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,4 +61,16 @@
         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,4 +152,38 @@
             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,9 +257,46 @@
             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,4 +311,18 @@
         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)
     # }}}
 
Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 27353)
+++ /issm/trunk-jpl/src/m/classes/stochasticforcing.m	(revision 27354)
@@ -120,4 +120,5 @@
 			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,4 +136,11 @@
 				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,4 +173,28 @@
 					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,4 +248,7 @@
 					end
 					if(strcmp(field,'FrontalForcingsRignotarma'))
+						dimensions(ind) = md.frontalforcings.num_basins;
+					end
+					if(strcmp(field,'FrontalForcingsSubglacialDischargearma'))
 						dimensions(ind) = md.frontalforcings.num_basins;
 					end
@@ -274,4 +309,5 @@
 		'FrictionSchoofWaterPressure',...
 		'FrontalForcingsRignotarma',...
+		'FrontalForcingsSubglacialDischargearma',...
 		'SMBarma',...
 		'SMBforcing'
@@ -286,4 +322,5 @@
 		'frictionschoof',...
 		'frontalforcingsrignotarma',...
+		'frontalforcingsrignotarma',...
 		'SMBarma',...
 		'SMBforcing'
Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/stochasticforcing.py	(revision 27353)
+++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py	(revision 27354)
@@ -36,9 +36,13 @@
         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,4 +126,5 @@
         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,4 +138,9 @@
             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,4 +162,17 @@
             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,4 +210,6 @@
                     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,4 +262,5 @@
                      'FrictionSchoofWaterPressure': 'frictionschoof',
                      'FrontalForcingsRignotarma': 'frontalforcingsrignotarma',
+                     'FrontalForcingsSubglacialDischargearma': 'frontalforcingsrignotarma',
                      'SMBarma': 'SMBarma',
                      'SMBforcing': 'SMBforcing'}
Index: /issm/trunk-jpl/test/NightlyRun/test543.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test543.m	(revision 27353)
+++ /issm/trunk-jpl/test/NightlyRun/test543.m	(revision 27354)
@@ -67,4 +67,19 @@
 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,6 +99,22 @@
 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,9 +123,10 @@
 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;
Index: /issm/trunk-jpl/test/NightlyRun/test543.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test543.py	(revision 27353)
+++ /issm/trunk-jpl/test/NightlyRun/test543.py	(revision 27354)
@@ -72,4 +72,19 @@
 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,4 +104,18 @@
 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,10 +132,12 @@
 
 #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,5 +145,5 @@
 # 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