Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24350)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24351)
@@ -246,4 +246,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       ElementResponse(IssmDouble* presponse,int response_enum)=0;
 		virtual void       ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24350)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24351)
@@ -1929,7 +1929,46 @@
 		if(!input) return input;
 
+		this->InputServe(input);
+		return input;
+	}
+	else{
+		TriaInput2* input = this->inputs2->GetTriaInput(inputenum);
+		if(!input) return input;
+
+		this->InputServe(input);
+		return input;
+	}
+}/*}}}*/
+Input2*    Tria::GetInput2(int inputenum,IssmDouble time){/*{{{*/
+
+	/*Get Input from dataset*/
+	if(this->iscollapsed){
+		PentaInput2* input = this->inputs2->GetPentaInput(inputenum,time);
+		if(!input) return input;
+
+		this->InputServe(input);
+		return input;
+	}
+	else{
+		TriaInput2* input = this->inputs2->GetTriaInput(inputenum,time);
+		if(!input) return input;
+
+		this->InputServe(input);
+		return input;
+	}
+}/*}}}*/
+void       Tria::InputServe(Input2* input_in){/*{{{*/
+
+	/*Return NULL pointer if input is NULL*/
+	if(!input_in) return;
+
+	/*Get Input from dataset*/
+	if(this->iscollapsed){
+		_assert_(input_in->ObjectEnum()==PentaInput2Enum);
+		PentaInput2* input = xDynamicCast<PentaInput2*>(input_in);
+
 		/*Intermediaries*/
 		int numindices;
-		int indices[7];
+		int indices[3];
 
 		/*Check interpolation*/
@@ -1955,78 +1994,9 @@
 		/*Flag as collapsed for later use*/
 		input->SetServeCollapsed(true);
-
-		return input;
+		return;
 	}
 	else{
-		TriaInput2* input = this->inputs2->GetTriaInput(inputenum);
-		if(!input) return input;
-
-		/*Intermediaries*/
-		int numindices;
-		int indices[7];
-
-		/*Check interpolation*/
-		int interpolation = input->GetInterpolation();
-		switch(interpolation){
-			case P0Enum:
-				numindices = 1;
-				indices[0] = this->lid;
-				input->Serve(numindices,&indices[0]);
-				break;
-			case P1Enum:
-				numindices = 3;
-				for(int i=0;i<3;i++) indices[i] = vertices[i]->lid;
-				input->Serve(numindices,&indices[0]);
-				break;
-			case P1DGEnum:
-				numindices = 3;
-				input->Serve(this->lid,numindices);
-				break;
-			default:
-				input->Serve(this->lid,this->GetNumberOfNodes(interpolation));
-		}
-
-		return input;
-	}
-}/*}}}*/
-Input2*    Tria::GetInput2(int inputenum,IssmDouble time){/*{{{*/
-
-	/*Get Input from dataset*/
-	if(this->iscollapsed){
-		PentaInput2* input = this->inputs2->GetPentaInput(inputenum,time);
-		if(!input) return input;
-
-		/*Intermediaries*/
-		int numindices;
-		int indices[3];
-
-		/*Check interpolation*/
-		int interpolation = input->GetInterpolation();
-		switch(interpolation){
-			case P0Enum:
-				numindices = 1;
-				indices[0] = this->lid;
-				input->Serve(numindices,&indices[0]);
-				break;
-			case P1Enum:
-				numindices = 3;
-				for(int i=0;i<3;i++) indices[i] = vertices[i]->lid;
-				input->Serve(numindices,&indices[0]);
-				break;
-			case P1DGEnum:
-			case P1bubbleEnum:
-				input->ServeCollapsed(this->lid,this->iscollapsed);
-				break;
-			default: _error_("interpolation "<<EnumToStringx(interpolation)<<" not supported");
-		}
-
-		/*Flag as collapsed for later use*/
-		input->SetServeCollapsed(true);
-
-		return input;
-	}
-	else{
-		TriaInput2* input = this->inputs2->GetTriaInput(inputenum,time);
-		if(!input) return input;
+		_assert_(input_in->ObjectEnum()==TriaInput2Enum);
+		TriaInput2* input = xDynamicCast<TriaInput2*>(input_in);
 
 		/*Intermediaries*/
@@ -2053,6 +2023,5 @@
 			default: _error_("interpolation "<<EnumToStringx(interpolation)<<" not supported");
 		}
-
-		return input;
+		return;
 	}
 }/*}}}*/
@@ -2129,4 +2098,125 @@
 	return datasetinput;
 }/*}}}*/
+void       Tria::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time){/*{{{*/
+
+	_assert_(end_time>start_time);
+
+	/*Intermediaries*/
+	IssmDouble averaged_values[NUMVERTICES];
+	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;
+	TransientInput2* transient_input  = this->inputs2->GetTransientInput(transientinput_enum);
+	transient_input->GetAllTimes(&timesteps,&numtimesteps);
+
+	/*go through the timesteps, and grab offset for start and end*/
+	found=binary_search(&start_offset,start_time,timesteps,numtimesteps);
+	if(!found) _error_("Input not found (is TransientInput sorted ?)");
+	found=binary_search(&end_offset,end_time,timesteps,numtimesteps);
+	if(!found) _error_("Input not found (is TransientInput sorted ?)");
+
+	Gauss* gauss=this->NewGauss();
+
+	/*stack the input for each timestep in the slice*/
+	int offset = start_offset;
+	while(offset <= end_offset ){
+
+		if(offset==-1){
+			/*get values for the first time: */
+			_assert_(start_time<timesteps[0]);
+			TriaInput2* input = transient_input->GetTriaInput(0);
+			_assert_(input->GetInterpolation()==P1Enum);
+			this->InputServe(input);
+			for(int iv=0;iv<NUMVERTICES;iv++){
+				gauss->GaussVertex(iv);
+				input->GetInputValue(&current_values[iv],gauss);
+			}
+		}
+		else{
+			TriaInput2* input = transient_input->GetTriaInput(offset);
+			_assert_(input->GetInterpolation()==P1Enum);
+			this->InputServe(input);
+			for(int iv=0;iv<NUMVERTICES;iv++){
+				gauss->GaussVertex(iv);
+				input->GetInputValue(&current_values[iv],gauss);
+			}
+		}
+
+		/*Step between current offset and next*/
+		if(offset==-1){
+			dt = timesteps[0] - start_time;
+		}
+		else if(offset = numtimesteps-1){
+			dt = end_time - timesteps[offset];
+		}
+		else{
+			dt = timesteps[offset+1] - timesteps[offset]; _assert_(dt>0.);
+		}
+
+		switch(averaging_method){
+			case 0: /*Arithmetic mean*/
+				if(offset==start_offset){
+					for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
+				}
+				else{
+					for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] += dt*current_values[iv];
+				}
+				break;
+			case 1: /*Geometric mean*/
+				if(offset==start_offset){
+					for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
+				}
+				else{
+					for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] += dt*current_values[iv];
+				}
+				break;
+			case 2: /*Harmonic mean*/
+				if(offset==start_offset){
+					for(int iv=0;iv<NUMVERTICES;iv++){
+						_assert_(current_values[iv]>1.e-50);
+						averaged_values[iv]  = dt*1./current_values[iv];
+					}
+				}
+				else{
+					for(int iv=0;iv<NUMVERTICES;iv++){
+						_assert_(current_values[iv]>1.e-50);
+						averaged_values[iv]  = dt*1./current_values[iv];
+					}
+				}
+				break;
+			default:
+				_error_("averaging method is not recognised");
+		}
+
+		offset+=1;
+	}
+
+	/*Integration done, now normalize*/
+	switch(averaging_method){
+		case 0: //Arithmetic mean
+			for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] =  averaged_values[iv]/(end_time - start_time);
+			break;
+		case 1: /*Geometric mean*/
+			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));
+			break;
+		default:
+			_error_("averaging method is not recognised");
+	}
+
+	this->AddInput2(averagedinput_enum,&averaged_values[0],P1Enum);
+
+	/*Cleanup*/
+	delete gauss;
+	xDelete<IssmDouble>(timesteps);
+}
+/*}}}*/
 void       Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24350)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24351)
@@ -176,4 +176,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);
 		IssmDouble     GetArea(void);
 		IssmDouble     GetHorizontalSurfaceArea(void);
@@ -221,4 +222,5 @@
 		void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
 		void           SetTemporaryElementType(int element_type_in){_error_("not implemented yet");};
+		void           InputServe(Input2* input_in);
 		Seg*	         SpawnSeg(int index1,int index2);
 		IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 24350)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 24351)
@@ -19,4 +19,5 @@
 #include "./Inputs2/DatasetInput2.h"
 #include "./Inputs2/ElementInput2.h"
+#include "./Inputs2/TransientInput2.h"
 
 #if _HAVE_CODIPACK_
@@ -5154,26 +5155,23 @@
 		else{
 			for(int j=0;j<elements->Size();j++){
-				/*Intermediaries*/
-				_error_("TODO");
+
+				/*Get the right transient input*/
 				Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
-				Input* input=element->inputs->GetInput(transientinput_enum[i]); _assert_(input); //this is the enum stack
-				TransientInput* stacking_input=xDynamicCast<TransientInput*>(input);
-
-				int  numvertices = element->GetNumberOfVertices();
-				IssmDouble* N=xNew<IssmDouble>(numvertices);
-				element->GetInputListOnVertices(&N[0],input_enum[i]);   //this is the enum to stack
+				TransientInput2* transientinput = this->inputs2->GetTransientInput(transientinput_enum[i]); _assert_(transientinput);
+
+				/*Get values and lid list*/
+				const int   numvertices = element->GetNumberOfVertices();
+				IssmDouble* values=xNew<IssmDouble>(numvertices);
+				int        *vertexlids = xNew<int>(numvertices);
+				element->GetInputListOnVertices(&values[0],input_enum[i]);   //this is the enum to stack
+				element->GetVerticesLidList(vertexlids);
+
 				switch(element->ObjectEnum()){
-				case TriaEnum:
-					stacking_input->AddTimeInput(new TriaInput(transientinput_enum[i],&N[0],P1Enum),subtime);
-					break;
-				case PentaEnum:
-					stacking_input->AddTimeInput(new PentaInput(transientinput_enum[i],&N[0],P1Enum),subtime);
-					break;
-				case TetraEnum:
-					stacking_input->AddTimeInput(new TetraInput(transientinput_enum[i],&N[0],P1Enum),subtime);
-					break;
-				default: _error_("Not implemented yet");
+					case TriaEnum:  transientinput->AddTriaTimeInput( subtime,numvertices,vertexlids,values,P1Enum); break;
+					case PentaEnum: transientinput->AddPentaTimeInput(subtime,numvertices,vertexlids,values,P1Enum); break;
+					default: _error_("Not implemented yet");
 				}
-				xDelete<IssmDouble>(N);
+				xDelete<IssmDouble>(values);
+				xDelete<int>(vertexlids);
 			}
 		}
@@ -5184,31 +5182,10 @@
 
 	for(int i=0;i<numoutputs;i++){
-		if(transientinput_enum[i]<0){
-			_error_("Can't deal with non enum fields for result Stack");
-		}
-		else{
-			for(int j=0;j<this->elements->Size();j++){
-				Element*    element       = xDynamicCast<Element*>(elements->GetObjectByOffset(j));
-				int         numnodes      = element->GetNumberOfNodes();
-				IssmDouble* time_averaged = xNew<IssmDouble>(numnodes);
-				Gauss*      gauss         = element->NewGauss();
-
-				Input*      input         = element->GetInput(transientinput_enum[i]); _assert_(input);
-				TransientInput* transient_input=xDynamicCast<TransientInput*>(input);
-
-				for(int iv=0;iv<numnodes;iv++){
-					gauss->GaussNode(element->FiniteElement(),iv);
-					transient_input->GetInputAverageOverTimeSlice(&time_averaged[iv],gauss,init_time,end_time);
-				}
-
-				element->AddInput2(averagedinput_enum[i],&time_averaged[0],element->GetElementType());
-				xDelete<IssmDouble>(time_averaged);
-				delete gauss;
-				xDelete<IssmDouble>(time_averaged);
-			}
-		}
-	}
-}
-/*}}}*/
+		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);
+		}
+	}
+}/*}}}*/
 #ifdef _HAVE_JAVASCRIPT_
 FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/
Index: /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp	(revision 24350)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp	(revision 24351)
@@ -97,6 +97,6 @@
 	for(i=0;i<this->numtimesteps;i++){
 		_printf_("   time: " << this->timesteps[i]<<"  ");
-		//((Input*)this->inputs->GetObjectByOffset(i))->Echo();
-		_error_("not implemented");
+		if(this->inputs[i]) this->inputs[i]->Echo();
+		else                _printf_(" NOT SET! \n");
 	}
 }
@@ -136,22 +136,34 @@
 
 /*Intermediary*/
-void TransientInput2::AddTimeInput(Input2* input,IssmDouble time){/*{{{*/
-
-	/*insert values at time step: */
-	if (this->numtimesteps>0 && time<=this->timesteps[this->numtimesteps-1]) _error_("timestep values must increase sequentially");
-
-	//copy timesteps, add the new time, delete previous timesteps, and add the new input: inputs->AddObject(input);
-	IssmDouble* old_timesteps=NULL;
-
+void TransientInput2::AddTriaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/
+
+	/*Check whether this is the last time step that we have*/
+	if(this->numtimesteps){
+		if(this->timesteps[this->numtimesteps-1]>time-1.e-5 && this->timesteps[this->numtimesteps-1]<time+1.e-5){
+			this->AddTriaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in);
+			return;
+		}
+	}
+
+	/*This is a new time step! we need to add it to the list*/
+	if(this->numtimesteps>0 && time<this->timesteps[this->numtimesteps-1]) _error_("timestep values must increase sequentially");
+
+	IssmDouble *old_timesteps = NULL;
+	Input2    **old_inputs    = NULL;
 	if (this->numtimesteps > 0){
 		old_timesteps=xNew<IssmDouble>(this->numtimesteps);
 		xMemCpy(old_timesteps,this->timesteps,this->numtimesteps);
-		xDelete(this->timesteps);
+		xDelete<IssmDouble>(this->timesteps);
+		old_inputs=xNew<Input2*>(this->numtimesteps);
+		xMemCpy(old_inputs,this->inputs,this->numtimesteps);
+		xDelete<Input2*>(this->inputs);
 	}
 
 	this->numtimesteps=this->numtimesteps+1;
 	this->timesteps=xNew<IssmDouble>(this->numtimesteps);
+	this->inputs   = xNew<Input2*>(this->numtimesteps);
 
 	if (this->numtimesteps > 1){
+		xMemCpy(this->inputs,old_inputs,this->numtimesteps-1);
 		xMemCpy(this->timesteps,old_timesteps,this->numtimesteps-1);
 		xDelete(old_timesteps);
@@ -159,7 +171,13 @@
 
 	/*go ahead and plug: */
-	this->timesteps[this->numtimesteps-1]=time;
-	//inputs->AddObject(input);
-	_error_("not implemented...");
+	this->timesteps[this->numtimesteps-1] = time;
+	this->inputs[this->numtimesteps-1]    = NULL;
+	this->AddTriaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in);
+
+}
+/*}}}*/
+void TransientInput2::AddPentaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/
+
+	_error_("not implemented yet, look at TransientInput2::AddTriaTimeInput");
 
 }
@@ -201,4 +219,16 @@
 }
 /*}}}*/
+void TransientInput2::GetAllTimes(IssmDouble** ptimesteps,int* pnumtimesteps){/*{{{*/
+
+	if(ptimesteps){
+		*ptimesteps=xNew<IssmDouble>(this->numtimesteps);
+		xMemCpy(*ptimesteps,this->timesteps,this->numtimesteps);
+	}
+	if(pnumtimesteps){
+		*pnumtimesteps = this->numtimesteps;
+	}
+
+}
+/*}}}*/
 TriaInput2* TransientInput2::GetTriaInput(){/*{{{*/
 
@@ -223,4 +253,19 @@
 }
 /*}}}*/
+TriaInput2* TransientInput2::GetTriaInput(int offset){/*{{{*/
+
+	/*Check offset*/
+	if(offset<0 || offset>this->numtimesteps-1){
+		_error_("Cannot return input for offset "<<offset);
+	}
+	Input2* input = this->inputs[offset];
+
+	/*Cast and return*/
+	_assert_(input);
+	if(input->ObjectEnum()!=TriaInput2Enum) _error_("Cannot return a TriaInput2");
+	return xDynamicCast<TriaInput2*>(input);
+
+}
+/*}}}*/
 PentaInput2* TransientInput2::GetPentaInput(){/*{{{*/
 
@@ -241,4 +286,19 @@
 	}
 	return xDynamicCast<PentaInput2*>(this->current_input);
+
+}
+/*}}}*/
+PentaInput2* TransientInput2::GetPentaInput(int offset){/*{{{*/
+
+
+	/*Check offset*/
+	if(offset<0 || offset>this->numtimesteps-1){
+		_error_("Cannot return input for offset "<<offset);
+	}
+	Input2* input = this->inputs[offset];
+
+	/*Cast and return*/
+	if(input->ObjectEnum()!=PentaInput2Enum) _error_("Cannot return a PentaInput2");
+	return xDynamicCast<PentaInput2*>(input);
 
 }
@@ -293,5 +353,5 @@
 
 		/*If already processed return*/
-		if(this->current_step>this_step-1.e-5 && this->current_step<this_step-1.e-5) return;
+		if(this->current_step>this_step-1.e-5 && this->current_step<this_step+1.e-5) return;
 
 		/*Prepare input*/
Index: /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h	(revision 24350)
+++ /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.h	(revision 24351)
@@ -31,5 +31,7 @@
 		TransientInput2(int in_enum_type,int nbe,int nbv,IssmDouble* times,int N);
 		~TransientInput2();
-		void AddTimeInput(Input2* input,IssmDouble time);
+		void AddTimeInput(Input2* input,IssmDouble time); /*FIXME: remove!*/
+		void AddTriaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in);
+		void AddPentaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in);
 		void AddTriaTimeInput(int step,int numindices,int* indices,IssmDouble* values_in,int interp_in);
 		void AddPentaTimeInput(int step,int numindices,int* indices,IssmDouble* values_in,int interp_in);
@@ -45,8 +47,11 @@
 		/*}}}*/
 		/*TransientInput2 management:*/
+		void         GetAllTimes(IssmDouble** ptimesteps,int* pnumtimesteps);
 		TriaInput2*  GetTriaInput();
 		TriaInput2*  GetTriaInput(IssmDouble time);
+		TriaInput2*  GetTriaInput(int offset);
 		PentaInput2* GetPentaInput();
 		PentaInput2* GetPentaInput(IssmDouble time);
+		PentaInput2* GetPentaInput(int offset);
 		Input2*      GetTimeInput(IssmDouble time){_error_("This should not happen!");};
 		IssmDouble   GetTimeByOffset(int offset);