Index: ../trunk-jpl/src/c/cores/smb_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/smb_core.cpp	(revision 24789)
+++ ../trunk-jpl/src/c/cores/smb_core.cpp	(revision 24790)
@@ -54,7 +54,7 @@
 
 	/*if yes compute necessary intermediaries and start looping*/
 	if (dtslices>1){
-		int        substep;
+		int        substep,smb_averaging;
 		IssmDouble global_time,subtime,yts;
 		IssmDouble dt,subdt;
 
@@ -61,6 +61,7 @@
 		femmodel->parameters->FindParam(&global_time,TimeEnum);
 		femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
 		femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
+		femmodel->parameters->FindParam(&smb_averaging,SmbAveragingEnum);
 
 		subtime=global_time-dt; //getting the time back to the start of the timestep
 		subdt=dt/dtslices; //computing substep from dt and a divider
@@ -82,7 +83,7 @@
 		}
 		delete analysis;
       /*averaging the transient input*/
-		femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput);
+		femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,smb_averaging);
 		/*and reset timesteping variables to original*/
 		femmodel->parameters->SetParam(global_time,TimeEnum);
 		femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
Index: ../trunk-jpl/src/c/cores/hydrology_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 24789)
+++ ../trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 24790)
@@ -62,7 +62,7 @@
 			femmodel->parameters->FindParam(&dtslices,HydrologyStepsPerStepEnum);
 
 			if(dtslices>1){
-				int        substep, numaveragedinput;
+				int        substep, numaveragedinput, hydro_averaging;
 				IssmDouble global_time, subtime, yts;
 				IssmDouble dt, subdt;
 
@@ -69,6 +69,7 @@
             femmodel->parameters->FindParam(&global_time,TimeEnum);
             femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
             femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
+				femmodel->parameters->FindParam(&hydro_averaging,HydrologyAveragingEnum);
 
 				subtime=global_time-dt; //getting the time back to the start of the timestep
 				subdt=dt/dtslices; //computing hydro dt from dt and a divider
@@ -116,7 +117,7 @@
 					femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
 				}
 				/*averaging the stack*/
-				femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput);
+				femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging);
 
 				/*And reseting to global time*/
 				femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
Index: ../trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.h	(revision 24789)
+++ ../trunk-jpl/src/c/classes/FemModel.h	(revision 24790)
@@ -175,7 +175,7 @@
 		void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count);
 		void InitTransientInputx(int* transientinput_enum,int numoutputs);
 		void StackTransientInputx(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs);
-		void AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs);
+		void AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method);
 		void UpdateConstraintsx(void);
 		int  UpdateVertexPositionsx(void);
 
Index: ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.h
===================================================================
--- ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.h	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.h	(revision 24790)
@@ -37,7 +37,9 @@
 		TriaInput2* GetTriaInput(){return this;};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss);
 		void Scale(IssmDouble scalar);
+		void Pow(IssmDouble scalar);
 		void AXPY(Input2* xinput,IssmDouble scalar);
+		void PointWiseMult(Input2* xinput);
 		void Serve(int numindices,int* indices);
 		void Serve(int row,int numindices);
 		void ServeCollapsed(int row,int id0,int in1);
Index: ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp	(revision 24790)
@@ -154,8 +154,8 @@
 void PentaInput2::SetInput(int interp_in,int row,IssmDouble value_in){/*{{{*/
 
 	_assert_(this);
-	_assert_(row>=0); 
-	_assert_(row<this->M); 
+	_assert_(row>=0);
+	_assert_(row<this->M);
 	_assert_(this->N==1);
 
 	this->values[row] = value_in;
@@ -169,8 +169,8 @@
 		_assert_(this->N==1);
 		for(int i=0;i<numindices;i++){
 			int row = indices[i];
-			_assert_(row>=0); 
-			_assert_(row<this->M); 
+			_assert_(row>=0);
+			_assert_(row<this->M);
 			this->values[row] = values_in[i];
 		}
 	}
@@ -178,8 +178,8 @@
 		this->Reset(interp_in);
 		for(int i=0;i<numindices;i++){
 			int row = indices[i];
-			_assert_(row>=0); 
-			_assert_(row<this->M); 
+			_assert_(row>=0);
+			_assert_(row<this->M);
 			this->values[row] = values_in[i];
 		}
 	}
@@ -212,8 +212,8 @@
 
 	for(int i=0;i<numindices;i++){
 		int row = indices[i];
-		_assert_(row>=0); 
-		_assert_(row<this->M); 
+		_assert_(row>=0);
+		_assert_(row<this->M);
 		this->element_values[i] = this->values[row];
 	}
 
@@ -388,6 +388,12 @@
 	for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*this->element_values[i];
 }
 /*}}}*/
+void PentaInput2::Pow(IssmDouble alpha){/*{{{*/
+
+	for(int i=0;i<this->M*this->N;i++) this->values[i] = pow(this->values[i],alpha);
+	for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = pow(this->element_values[i],alpha);
+}
+/*}}}*/
 void PentaInput2::AXPY(Input2* xinput,IssmDouble alpha){/*{{{*/
 
 	/*xinput is of the same type, so cast it: */
@@ -400,5 +406,20 @@
 	for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xpentainput->element_values[i] + this->element_values[i];
 }
 /*}}}*/
+void PentaInput2::PointWiseMult(Input2* xinput){/*{{{*/
 
+	/*xinput is of the same type, so cast it: */
+	if(xinput->ObjectEnum()!=PentaInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
+	PentaInput2* xpentainput=xDynamicCast<PentaInput2*>(xinput);
+	if(xpentainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
+
+	/* we need to check that the vector sizes are identical*/
+	if(xpentainput->M!=this->M||xpentainput->N!=this->N) _error_("Operation not permitted because the inputs have different sizes");
+
+	/*Carry out the pointwise operation depending on type:*/
+	for(int i=0;i<this->M*this->N;i++) this->values[i] = xpentainput->values[i] * this->values[i];
+	for(int i=0;i<PentaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = xpentainput->element_values[i] * this->element_values[i];
+}
+/*}}}*/
+
 /*Object functions*/
Index: ../trunk-jpl/src/c/classes/Inputs2/Inputs2.h
===================================================================
--- ../trunk-jpl/src/c/classes/Inputs2/Inputs2.h	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Inputs2/Inputs2.h	(revision 24790)
@@ -53,7 +53,7 @@
 		SegInput2*       GetSegInput(int enum_type);
 		TriaInput2*      GetTriaInput(int enum_type);
 		TriaInput2*      GetTriaInput(int enum_type,IssmDouble time);
-		TriaInput2*      GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time);
+		TriaInput2*      GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time,int averaging_method);
 		PentaInput2*     GetPentaInput(int enum_type);
 		PentaInput2*     GetPentaInput(int enum_type,IssmDouble time);
 		TransientInput2* GetTransientInput(int enum_type);
Index: ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.h
===================================================================
--- ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.h	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Inputs2/PentaInput2.h	(revision 24790)
@@ -36,7 +36,9 @@
 		PentaInput2* GetPentaInput(){return this;};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss);
 		void Scale(IssmDouble scalar);
+		void Pow(IssmDouble scalar);
 		void AXPY(Input2* xinput,IssmDouble scalar);
+		void PointWiseMult(Input2* xinput);
 		void Serve(int numindices,int* indices);
 		void Serve(int row,int numindices);
 		void ServeCollapsed(int row,int state);
Index: ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.h
===================================================================
--- ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.h	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.h	(revision 24790)
@@ -49,7 +49,7 @@
 		void         GetAllTimes(IssmDouble** ptimesteps,int* pnumtimesteps);
 		TriaInput2*  GetTriaInput();
 		TriaInput2*  GetTriaInput(IssmDouble time);
-		TriaInput2*  GetTriaInput(IssmDouble start_time,IssmDouble end_time);
+		TriaInput2*  GetTriaInput(IssmDouble start_time,IssmDouble end_time,int averaging_method);
 		TriaInput2*  GetTriaInput(int offset);
 		PentaInput2* GetPentaInput();
 		PentaInput2* GetPentaInput(IssmDouble time);
@@ -58,7 +58,7 @@
 		IssmDouble   GetTimeByOffset(int offset);
 		int          GetTimeInputOffset(IssmDouble time);
 		void         SetCurrentTimeInput(IssmDouble time);
-		void         SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time);
+		void         SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time,int averaging_method);
 		/*numerics:*/
 
 };
Index: ../trunk-jpl/src/c/classes/Inputs2/SegInput2.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Inputs2/SegInput2.cpp	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Inputs2/SegInput2.cpp	(revision 24790)
@@ -136,8 +136,8 @@
 void SegInput2::SetInput(int interp_in,int row,IssmDouble value_in){/*{{{*/
 
 	_assert_(this);
-	_assert_(row>=0); 
-	_assert_(row<this->M); 
+	_assert_(row>=0);
+	_assert_(row<this->M);
 	_assert_(this->N==1);
 
 	this->values[row] = value_in;
@@ -151,8 +151,8 @@
 		_assert_(this->N==1);
 		for(int i=0;i<numindices;i++){
 			int row = indices[i];
-			_assert_(row>=0); 
-			_assert_(row<this->M); 
+			_assert_(row>=0);
+			_assert_(row<this->M);
 			this->values[row] = values_in[i];
 		}
 	}
@@ -160,8 +160,8 @@
 		_assert_(this->N==1);
 		for(int i=0;i<numindices;i++){
 			int row = indices[i];
-			_assert_(row>=0); 
-			_assert_(row<this->M); 
+			_assert_(row>=0);
+			_assert_(row<this->M);
 			this->values[row] = values_in[i];
 		}
 	}
@@ -169,8 +169,8 @@
 		this->Reset(interp_in);
 		for(int i=0;i<numindices;i++){
 			int row = indices[i];
-			_assert_(row>=0); 
-			_assert_(row<this->M); 
+			_assert_(row>=0);
+			_assert_(row<this->M);
 			this->values[row] = values_in[i];
 		}
 	}
@@ -201,8 +201,8 @@
 
 	for(int i=0;i<numindices;i++){
 		int row = indices[i];
-		_assert_(row>=0); 
-		_assert_(row<this->M); 
+		_assert_(row>=0);
+		_assert_(row<this->M);
 		this->element_values[i] = this->values[row];
 	}
 
@@ -308,17 +308,38 @@
 	for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*this->element_values[i];
 }
 /*}}}*/
+void SegInput2::Pow(IssmDouble alpha){/*{{{*/
+
+	for(int i=0;i<this->M*this->N;i++) this->values[i] = pow(this->values[i],alpha);
+	for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = pow(this->element_values[i],alpha);
+}
+/*}}}*/
 void SegInput2::AXPY(Input2* xinput,IssmDouble alpha){/*{{{*/
 
 	/*xinput is of the same type, so cast it: */
 	if(xinput->ObjectEnum()!=SegInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
-	SegInput2* xtriainput=xDynamicCast<SegInput2*>(xinput);
-	if(xtriainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
+	SegInput2* xseginput=xDynamicCast<SegInput2*>(xinput);
+	if(xseginput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
 
 	/*Carry out the AXPY operation depending on type:*/
-	for(int i=0;i<this->M*this->N;i++) this->values[i] = alpha*xtriainput->values[i] + this->values[i];
-	for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xtriainput->element_values[i] + this->element_values[i];
+	for(int i=0;i<this->M*this->N;i++) this->values[i] = alpha*xseginput->values[i] + this->values[i];
+	for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xseginput->element_values[i] + this->element_values[i];
 }
 /*}}}*/
+void SegInput2::PointWiseMult(Input2* xinput){/*{{{*/
 
+	/*xinput is of the same type, so cast it: */
+	if(xinput->ObjectEnum()!=SegInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
+	SegInput2* xseginput=xDynamicCast<SegInput2*>(xinput);
+	if(xseginput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
+
+	/* we need to check that the vector sizes are identical*/
+	if(xseginput->M!=this->M||xseginput->N!=this->N) _error_("Operation not permitted because the inputs have different sizes");
+
+	/*Carry out the AXPY operation depending on type:*/
+	for(int i=0;i<this->M*this->N;i++) this->values[i] = xseginput->values[i] * this->values[i];
+	for(int i=0;i<SegRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = xseginput->element_values[i] * this->element_values[i];
+}
+/*}}}*/
+
 /*Object functions*/
Index: ../trunk-jpl/src/c/classes/Inputs2/SegInput2.h
===================================================================
--- ../trunk-jpl/src/c/classes/Inputs2/SegInput2.h	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Inputs2/SegInput2.h	(revision 24790)
@@ -34,7 +34,9 @@
 		SegInput2* GetSegInput(){return this;};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss);
 		void Scale(IssmDouble scalar);
+		void Pow(IssmDouble scalar);
 		void AXPY(Input2* xinput,IssmDouble scalar);
+		void PointWiseMult(Input2* xinput);
 		void Serve(int numindices,int* indices);
 		void Serve(int row,int numindices);
 		int  GetResultArraySize(void);
Index: ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.cpp	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Inputs2/TriaInput2.cpp	(revision 24790)
@@ -141,8 +141,8 @@
 void TriaInput2::SetInput(int interp_in,int row,IssmDouble value_in){/*{{{*/
 
 	_assert_(this);
-	_assert_(row>=0); 
-	_assert_(row<this->M); 
+	_assert_(row>=0);
+	_assert_(row<this->M);
 	_assert_(this->N==1);
 
 	this->values[row] = value_in;
@@ -156,8 +156,8 @@
 		_assert_(this->N==1);
 		for(int i=0;i<numindices;i++){
 			int row = indices[i];
-			_assert_(row>=0); 
-			_assert_(row<this->M); 
+			_assert_(row>=0);
+			_assert_(row<this->M);
 			this->values[row] = values_in[i];
 		}
 	}
@@ -165,8 +165,8 @@
 		_assert_(this->N==1);
 		for(int i=0;i<numindices;i++){
 			int row = indices[i];
-			_assert_(row>=0); 
-			_assert_(row<this->M); 
+			_assert_(row>=0);
+			_assert_(row<this->M);
 			this->values[row] = values_in[i];
 		}
 	}
@@ -174,8 +174,8 @@
 		this->Reset(interp_in);
 		for(int i=0;i<numindices;i++){
 			int row = indices[i];
-			_assert_(row>=0); 
-			_assert_(row<this->M); 
+			_assert_(row>=0);
+			_assert_(row<this->M);
 			this->values[row] = values_in[i];
 		}
 	}
@@ -206,8 +206,8 @@
 
 	for(int i=0;i<numindices;i++){
 		int row = indices[i];
-		_assert_(row>=0); 
-		_assert_(row<this->M); 
+		_assert_(row>=0);
+		_assert_(row<this->M);
 		this->element_values[i] = this->values[row];
 	}
 
@@ -363,6 +363,12 @@
 	for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*this->element_values[i];
 }
 /*}}}*/
+void TriaInput2::Pow(IssmDouble alpha){/*{{{*/
+
+	for(int i=0;i<this->M*this->N;i++) this->values[i] = pow(this->values[i],alpha);
+	for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = pow(this->element_values[i],alpha);
+}
+/*}}}*/
 void TriaInput2::AXPY(Input2* xinput,IssmDouble alpha){/*{{{*/
 
 	/*xinput is of the same type, so cast it: */
@@ -375,5 +381,20 @@
 	for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xtriainput->element_values[i] + this->element_values[i];
 }
 /*}}}*/
+void TriaInput2::PointWiseMult(Input2* xinput){/*{{{*/
 
+	/*xinput is of the same type, so cast it: */
+	if(xinput->ObjectEnum()!=TriaInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
+	TriaInput2* xtriainput=xDynamicCast<TriaInput2*>(xinput);
+	if(xtriainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
+
+	/* we need to check that the vector sizes are identical*/
+	if(xtriainput->M!=this->M||xtriainput->N!=this->N) _error_("Operation not permitted because the inputs have different sizes");
+
+	/*Carry out the AXPY operation depending on type:*/
+	for(int i=0;i<this->M*this->N;i++) this->values[i] = xtriainput->values[i] * this->values[i];
+	for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = xtriainput->element_values[i] * this->element_values[i];
+}
+/*}}}*/
+
 /*Object functions*/
Index: ../trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp	(revision 24790)
@@ -343,7 +343,7 @@
 		return input->GetTriaInput();
 	}
 }/*}}}*/
-TriaInput2* Inputs2::GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time){/*{{{*/
+TriaInput2* Inputs2::GetTriaInput(int enum_in,IssmDouble start_time,IssmDouble end_time,int averaging_method){/*{{{*/
 
 	/*Get input id*/
 	int id = EnumToIndex(enum_in);
@@ -353,7 +353,7 @@
 	if(!input) return NULL;
 
 	if(input->ObjectEnum()==TransientInput2Enum){
-		return xDynamicCast<TransientInput2*>(input)->GetTriaInput(start_time,end_time);
+		return xDynamicCast<TransientInput2*>(input)->GetTriaInput(start_time,end_time,averaging_method);
 	}
 	else{
 		_error_("Input "<<EnumToStringx(enum_in)<<" is not an TransientInput2");
Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24790)
@@ -177,7 +177,7 @@
 		void           AddInput2(int input_enum, IssmDouble* values, int interpolation_enum);
 		void           AddControlInput(int input_enum,Inputs2* inputs2,IoModel* iomodel,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id);
 		void           DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,Inputs2* inputs2,IoModel* iomodel,int input_enum);
-		void           CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time);
+		void           CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int averaging_method);
 		void           GetInputAveragesUpToCurrentTime(int input_enum,IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime);
 		IssmDouble     GetArea(void);
 		IssmDouble     GetHorizontalSurfaceArea(void);
@@ -188,7 +188,7 @@
 		int            GetElementType(void);
 		Input2*        GetInput2(int enumtype);
 		Input2*        GetInput2(int enumtype,IssmDouble time);
-		Input2*        GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time);
+		Input2*        GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method);
 		DatasetInput2* GetDatasetInput2(int inputenum);
 		void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		void           GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype);
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24790)
@@ -1007,7 +1007,7 @@
 	this->AddInput2(distanceenum,&ls[0],P1Enum);
 }
 /*}}}*/
-void       Penta::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time){/*{{{*/
+void       Penta::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method){/*{{{*/
 
 	_assert_(end_time>start_time);
 
@@ -1016,9 +1016,7 @@
 	IssmDouble current_values[NUMVERTICES];
 	IssmDouble dt;
 	int        found,start_offset,end_offset;
-	int        averaging_method=0;
 
-
 	/*Get transient input time steps*/
 	int         numtimesteps;
 	IssmDouble *timesteps    = NULL;
Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24790)
@@ -67,7 +67,7 @@
 		void           ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
 		void				CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum);
 		ElementMatrix* CreateBasalMassMatrix(void);
-		void           CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time);
+		void           CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time,int averaging_method);
 		void           ElementResponse(IssmDouble* presponse,int response_enum);
 		void           ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
 		int            FiniteElement(void);
@@ -85,7 +85,7 @@
 		void           GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
 		Input2*        GetInput2(int enumtype);
 		Input2*        GetInput2(int enumtype,IssmDouble time);
-		Input2*        GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
+		Input2*        GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet!");};
 		void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
 		void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
 		DatasetInput2* GetDatasetInput2(int inputenum);
Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Seg.h	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Elements/Seg.h	(revision 24790)
@@ -64,7 +64,7 @@
 		void		   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
 		Input2*     GetInput2(int enumtype);
 		Input2*     GetInput2(int enumtype,IssmDouble time);
-		Input2*     GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
+		Input2*     GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet!");};
 		void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
 		void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
 		void        GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");};
Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 24790)
@@ -68,7 +68,7 @@
 		void		   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
 		Input2*     GetInput2(int enumtype);
 		Input2*     GetInput2(int enumtype,IssmDouble time);
-		Input2*     GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time){_error_("not implemented yet!");};
+		Input2*     GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet!");};
 		void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
 		void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
 		void        GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
Index: ../trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 24790)
@@ -233,7 +233,7 @@
 		virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;
 		virtual void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum)=0;
 		virtual void       CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum){_error_("not implemented yet");};
-		virtual void       CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time){_error_("not implemented yet "<<this->ObjectEnum());};
+		virtual void       CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int averaging_method){_error_("not implemented yet "<<this->ObjectEnum());};
 		virtual void       ElementResponse(IssmDouble* presponse,int response_enum)=0;
 		virtual void       ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
 		virtual int        FiniteElement(void)=0;
@@ -249,7 +249,7 @@
 		virtual DatasetInput2* GetDatasetInput2(int inputenum){_error_("not implemented");};
 		virtual Input2*    GetInput2(int inputenum)=0;
 		virtual Input2*    GetInput2(int inputenum,IssmDouble time)=0;
-		virtual Input2*    GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time)=0;
+		virtual Input2*    GetInput2(int inputenum,IssmDouble start_time,IssmDouble end_time,int averaging_method)=0;
 		virtual void       GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");};
 		virtual void       GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
 		virtual void       GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value)=0;
Index: ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp	(revision 24790)
@@ -289,10 +289,10 @@
 
 }
 /*}}}*/
-TriaInput2* TransientInput2::GetTriaInput(IssmDouble start_time, IssmDouble end_time){/*{{{*/
+TriaInput2* TransientInput2::GetTriaInput(IssmDouble start_time, IssmDouble end_time, int averaging_method){/*{{{*/
 
 	/*Set current time input*/
-	this->SetAverageAsCurrentTimeInput(start_time,end_time);
+	this->SetAverageAsCurrentTimeInput(start_time,end_time,averaging_method);
 	_assert_(this->current_input);
 
 	/*Cast and return*/
@@ -420,13 +420,12 @@
 	}
 
 }/*}}}*/
-void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time){/*{{{*/
+void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time, int averaging_method){/*{{{*/
 
-	IssmDouble  dt;
+	IssmDouble  dt,durinv;
 	IssmPDouble eps=1.0e-6;
 	IssmDouble  dtsum=0;
 	int         found,start_offset,end_offset;
-	int         averaging_method = 0;
 
 	/*go through the timesteps, and grab offset for start and end*/
 	IssmDouble temp = start_time-eps;
@@ -462,9 +461,25 @@
 				}
 				break;
 			case 1: /*Geometric mean*/
-				_error_("Geometric not implemented yet");
+				if(offset==start_offset){
+					this->current_input = stepinput->copy();
+					this->current_input->Scale(dt);
+				}
+				else{
+					this->stepinput->Scale(dt);
+					this->current_input->PointWiseMult(stepinput)
+				}
+				break;
 			case 2: /*Harmonic mean*/
-				_error_("Harmonic not implemented yet");
+				if(offset==start_offset){
+					this->current_input = stepinput->copy();
+					this->current_input->Pow(-1);
+					this->current_input->Scale(dt);
+				}
+				else{
+					this->stepinput->Pow(-1);
+					this->current_input->AXPY(stepinput,dt);
+				}
 			default:
 				_error_("averaging method is not recognised");
 		}
@@ -471,15 +486,19 @@
 		dtsum+=dt;
 		offset+=1;
 	}
-		/*Integration done, now normalize*/
+	_assert_(dtsum>0);
+	durinv=1./dtsum;
+	/*Integration done, now normalize*/
 	switch(averaging_method){
 		case 0: //Arithmetic mean
-			this->current_input->Scale(1/(dtsum));
+			this->current_input->Scale(durinv);
 			break;
 		case 1: /*Geometric mean*/
-			_error_("Geometric not implemented yet");
+			this->current_input->Pow(durinv);
+			break;
 		case 2: /*Harmonic mean*/
-			_error_("Harmonic not implemented yet");
+			this->current_input->Scale(durinv);
+			this->current_input->Pow(-1);
 		default:
 			_error_("averaging method is not recognised");
 	}
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 24789)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 24790)
@@ -5274,12 +5274,12 @@
 	}
 }
 /*}}}*/
-void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs){ /*{{{*/
+void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method){ /*{{{*/
 
 	for(int i=0;i<numoutputs;i++){
 		for(int j=0;j<this->elements->Size();j++){
 			Element* element = xDynamicCast<Element*>(elements->GetObjectByOffset(j));
-			element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time);
+			element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time,averaging_method);
 		}
 	}
 }/*}}}*/
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24789)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24790)
@@ -1912,7 +1912,7 @@
 		return input;
 	}
 }/*}}}*/
-Input2*    Tria::GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time){/*{{{*/
+Input2*    Tria::GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time, int averaging_method){/*{{{*/
 
 	/*Get Input from dataset*/
 	if(this->iscollapsed){
@@ -1919,7 +1919,7 @@
 		_error_("Get Average input not implemented in Penta yet");
 	}
 	else{
-		TriaInput2* input = this->inputs2->GetTriaInput(inputenum,start_time, end_time);
+		TriaInput2* input = this->inputs2->GetTriaInput(inputenum,start_time, end_time, int averaging_method);
 		if(!input) return input;
 
 		this->InputServe(input);
@@ -2108,7 +2108,7 @@
 
 	return datasetinput;
 }/*}}}*/
-void       Tria::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time){/*{{{*/
+void       Tria::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time, int averaging_method){/*{{{*/
 
 	_assert_(end_time>start_time);
 
@@ -2117,9 +2117,7 @@
 	IssmDouble current_values[NUMVERTICES];
 	IssmDouble dt;
 	int        found,start_offset,end_offset;
-	int        averaging_method=0;
 
-
 	/*Get transient input time steps*/
 	int         numtimesteps;
 	IssmDouble *timesteps    = NULL;
@@ -2212,7 +2210,7 @@
 			for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = pow(averaged_values[iv], 1./(end_time - start_time));
 			break;
 		case 2: /*Harmonic mean*/
-			for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = 1./(averaged_values[iv]/(end_time - start_time));
+			for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = (end_time - start_time)/(averaged_values[iv];
 			break;
 		default:
 			_error_("averaging method is not recognised");
Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 24789)
+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 24790)
@@ -335,6 +335,7 @@
 	int        smb_model;
 	int        smbsubstepping;
 	int        hydrologysubstepping;
+	int        smb_averaging;
 	IssmDouble dt,scalar,sediment_storing;
 	IssmDouble water_head,sediment_transmitivity;
 	IssmDouble water_load,runoff_value,transfer;
@@ -358,6 +359,7 @@
 	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
 	basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
 	basalelement->FindParam(&smb_model,SmbEnum);
+	basalelement->FindParam(&smb_averaging,SmbAveragingEnum);
 
 	Input2*	sed_head_input   = basalelement->GetInput2(SedimentHeadSubstepEnum);
 	Input2*	epl_head_input   = basalelement->GetInput2(EplHeadSubstepEnum);
@@ -383,7 +385,7 @@
 			dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum, time); _assert_(dummy_input);
 		}
 		else{
-			dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time); _assert_(dummy_input);
+			dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time,smb_averaging); _assert_(dummy_input);
 		}
 		surface_runoff_input=xDynamicCast<Input2*>(dummy_input); _assert_(surface_runoff_input);
 	}