Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 24789)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 24790)
@@ -336,4 +336,5 @@
 	int        smbsubstepping;
 	int        hydrologysubstepping;
+	int        smb_averaging;
 	IssmDouble dt,scalar,sediment_storing;
 	IssmDouble water_head,sediment_transmitivity;
@@ -359,4 +360,5 @@
 	basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
 	basalelement->FindParam(&smb_model,SmbEnum);
+	basalelement->FindParam(&smb_averaging,SmbAveragingEnum);
 
 	Input2*	sed_head_input   = basalelement->GetInput2(SedimentHeadSubstepEnum);
@@ -384,5 +386,5 @@
 		}
 		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);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24790)
@@ -234,5 +234,5 @@
 		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;
@@ -250,5 +250,5 @@
 		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");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24790)
@@ -1008,5 +1008,5 @@
 }
 /*}}}*/
-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);
@@ -1017,6 +1017,4 @@
 	IssmDouble dt;
 	int        found,start_offset,end_offset;
-	int        averaging_method=0;
-
 
 	/*Get transient input time steps*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24790)
@@ -68,5 +68,5 @@
 		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);
@@ -86,5 +86,5 @@
 		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);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 24790)
@@ -65,5 +65,5 @@
 		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);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 24790)
@@ -69,5 +69,5 @@
 		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);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24790)
@@ -1913,5 +1913,5 @@
 	}
 }/*}}}*/
-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*/
@@ -1920,5 +1920,5 @@
 	}
 	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;
 
@@ -2109,5 +2109,5 @@
 	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);
@@ -2118,6 +2118,4 @@
 	IssmDouble dt;
 	int        found,start_offset,end_offset;
-	int        averaging_method=0;
-
 
 	/*Get transient input time steps*/
@@ -2213,5 +2211,5 @@
 			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:
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24790)
@@ -178,5 +178,5 @@
 		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);
@@ -189,5 +189,5 @@
 		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);
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 24790)
@@ -5275,10 +5275,10 @@
 }
 /*}}}*/
-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: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 24790)
@@ -176,5 +176,5 @@
 		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: /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.cpp	(revision 24790)
@@ -344,5 +344,5 @@
 	}
 }/*}}}*/
-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*/
@@ -354,5 +354,5 @@
 
 	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{
Index: /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.h	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/Inputs2.h	(revision 24790)
@@ -54,5 +54,5 @@
 		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);
Index: /issm/trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/PentaInput2.cpp	(revision 24790)
@@ -155,6 +155,6 @@
 
 	_assert_(this);
-	_assert_(row>=0); 
-	_assert_(row<this->M); 
+	_assert_(row>=0);
+	_assert_(row<this->M);
 	_assert_(this->N==1);
 
@@ -170,6 +170,6 @@
 		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];
 		}
@@ -179,6 +179,6 @@
 		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];
 		}
@@ -213,6 +213,6 @@
 	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];
 	}
@@ -389,4 +389,10 @@
 }
 /*}}}*/
+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){/*{{{*/
 
@@ -401,4 +407,19 @@
 }
 /*}}}*/
+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: /issm/trunk-jpl/src/c/classes/Inputs2/PentaInput2.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/PentaInput2.h	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/PentaInput2.h	(revision 24790)
@@ -37,5 +37,7 @@
 		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);
Index: /issm/trunk-jpl/src/c/classes/Inputs2/SegInput2.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/SegInput2.cpp	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/SegInput2.cpp	(revision 24790)
@@ -137,6 +137,6 @@
 
 	_assert_(this);
-	_assert_(row>=0); 
-	_assert_(row<this->M); 
+	_assert_(row>=0);
+	_assert_(row<this->M);
 	_assert_(this->N==1);
 
@@ -152,6 +152,6 @@
 		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];
 		}
@@ -161,6 +161,6 @@
 		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];
 		}
@@ -170,6 +170,6 @@
 		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];
 		}
@@ -202,6 +202,6 @@
 	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];
 	}
@@ -309,14 +309,35 @@
 }
 /*}}}*/
+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];
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Inputs2/SegInput2.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/SegInput2.h	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/SegInput2.h	(revision 24790)
@@ -35,5 +35,7 @@
 		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);
Index: /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp	(revision 24790)
@@ -290,8 +290,8 @@
 }
 /*}}}*/
-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);
 
@@ -421,11 +421,10 @@
 
 }/*}}}*/
-void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time){/*{{{*/
-
-	IssmDouble  dt;
+void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time, int averaging_method){/*{{{*/
+
+	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*/
@@ -463,7 +462,23 @@
 				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");
@@ -472,13 +487,17 @@
 		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: /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h	(revision 24790)
@@ -50,5 +50,5 @@
 		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();
@@ -59,5 +59,5 @@
 		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: /issm/trunk-jpl/src/c/classes/Inputs2/TriaInput2.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/TriaInput2.cpp	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/TriaInput2.cpp	(revision 24790)
@@ -142,6 +142,6 @@
 
 	_assert_(this);
-	_assert_(row>=0); 
-	_assert_(row<this->M); 
+	_assert_(row>=0);
+	_assert_(row<this->M);
 	_assert_(this->N==1);
 
@@ -157,6 +157,6 @@
 		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];
 		}
@@ -166,6 +166,6 @@
 		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];
 		}
@@ -175,6 +175,6 @@
 		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];
 		}
@@ -207,6 +207,6 @@
 	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];
 	}
@@ -364,4 +364,10 @@
 }
 /*}}}*/
+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){/*{{{*/
 
@@ -376,4 +382,19 @@
 }
 /*}}}*/
+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: /issm/trunk-jpl/src/c/classes/Inputs2/TriaInput2.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/TriaInput2.h	(revision 24789)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/TriaInput2.h	(revision 24790)
@@ -38,5 +38,7 @@
 		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);
Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 24789)
+++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 24790)
@@ -63,5 +63,5 @@
 
 			if(dtslices>1){
-				int        substep, numaveragedinput;
+				int        substep, numaveragedinput, hydro_averaging;
 				IssmDouble global_time, subtime, yts;
 				IssmDouble dt, subdt;
@@ -70,4 +70,5 @@
             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
@@ -117,5 +118,5 @@
 				}
 				/*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*/
Index: /issm/trunk-jpl/src/c/cores/smb_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/smb_core.cpp	(revision 24789)
+++ /issm/trunk-jpl/src/c/cores/smb_core.cpp	(revision 24790)
@@ -55,5 +55,5 @@
 	/*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;
@@ -62,4 +62,5 @@
 		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
@@ -83,5 +84,5 @@
 		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);