Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23543)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23544)
@@ -4821,13 +4821,13 @@
 }
 /*}}}*/
-void FemModel::InitTransientOutputx(int* input_enum,int numoutputs){ /*{{{*/
+void FemModel::InitTransientOutputx(int* stackedinput_enum,int numoutputs){ /*{{{*/
 
   for(int i=0;i<numoutputs;i++){
-		if(input_enum[i]<0){
+		if(stackedinput_enum[i]<0){
 			_error_("Can't deal with non enum fields for result Stack");
 		}
 		else{
 			for(int j=0;j<elements->Size();j++){
-				TransientInput* transient_input = new TransientInput(input_enum[i]);
+				TransientInput* transient_input = new TransientInput(stackedinput_enum[i]);
 				/*Intermediaries*/
 				Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
@@ -4850,11 +4850,10 @@
 				Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
 				TransientInput* stacking_input=NULL;
-				Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input);
-				Input* input_to_stack=element->GetInput(input_enum[i]); _assert_(input_to_stack);
+				Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input); //this is the enum stack
 				stacking_input=dynamic_cast<TransientInput*>(input);
 
 				int  numvertices = element->GetNumberOfVertices();
 				IssmDouble* N=xNew<IssmDouble>(numvertices);
-				element->GetInputListOnVertices(&N[0],input_enum[i]);
+				element->GetInputListOnVertices(&N[0],input_enum[i]);	//this is the enum to stack
 				switch(element->ObjectEnum()){
 				case TriaEnum:
@@ -4875,8 +4874,8 @@
 }
 /*}}}*/
-void FemModel::AverageTransientOutputx(int* input_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs){ /*{{{*/
+void FemModel::AverageTransientOutputx(int* stackedinput_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs){ /*{{{*/
 
   for(int i=0;i<numoutputs;i++){
-		if(input_enum[i]<0){
+		if(stackedinput_enum[i]<0){
 			_error_("Can't deal with non enum fields for result Stack");
 		}
@@ -4888,5 +4887,5 @@
 				IssmDouble* time_averaged=NULL;
 
-				Input* input=element->inputs->GetInput(input_enum[i]); _assert_(input);
+				Input* input=element->inputs->GetInput(stackedinput_enum[i]); _assert_(input);
 				TransientInput* stacking_input=NULL;
 				stacking_input=dynamic_cast<TransientInput*>(input);
@@ -4895,4 +4894,93 @@
 				element->AddInput(averagedinput_enum[i],&time_averaged[0],P1Enum);
 				xDelete<IssmDouble>(time_averaged);
+			}
+		}
+	}
+}
+/*}}}*/
+void FemModel::InitMeanOutputx(int* stackedinput_enum,int numoutputs){ /*{{{*/
+
+  for(int i=0;i<numoutputs;i++){
+		if(stackedinput_enum[i]<0){
+			_error_("Can't deal with non enum fields for result Stack");
+		}
+		else{
+			for(int j=0;j<elements->Size();j++){
+				/*Intermediaries*/
+				Element*   element            =xDynamicCast<Element*>(elements->GetObjectByOffset(j));
+				int        numvertices        =element->GetNumberOfVertices();
+				IssmDouble zeros[numvertices] ={0.0};
+				switch(element->ObjectEnum()){
+				case TriaEnum:
+					element->inputs->AddInput(new TriaInput(stackedinput_enum[i],&zeros[0],P1Enum));
+					break;
+				case PentaEnum:
+					element->inputs->AddInput(new PentaInput(stackedinput_enum[i],&zeros[0],P1Enum));
+					break;
+				case TetraEnum:
+					element->inputs->AddInput(new TetraInput(stackedinput_enum[i],&zeros[0],P1Enum));
+					break;
+				default: _error_("Not implemented yet");
+				}
+			}
+		}
+	}
+}
+/*}}}*/
+void FemModel::SumOutputx(int* input_enum,int* stackedinput_enum,int numoutputs){ /*{{{*/
+
+	//First get sub-timestep
+	IssmDouble hydrodt;
+	this->parameters->FindParam(&hydrodt,HydrologydtEnum);
+
+  for(int i=0;i<numoutputs;i++){
+		if(input_enum[i]<0){
+			_error_("Can't deal with non enum fields for result Stack");
+		}
+		else{
+			for(int j=0;j<elements->Size();j++){
+				/*Intermediaries*/
+				Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
+				int  numvertices = element->GetNumberOfVertices();
+				IssmDouble* values_to_add=xNew<IssmDouble>(numvertices);
+				IssmDouble* existing_values=xNew<IssmDouble>(numvertices);
+				element->GetInputListOnVertices(&values_to_add[0],input_enum[i]); //those are the values to add
+				element->GetInputListOnVertices(&existing_values[0],stackedinput_enum[i]); //those are the values to add
+				for(int k=0;k<numvertices;k++){
+					existing_values[k]+=values_to_add[k]*hydrodt;
+				}
+				element->AddInput(stackedinput_enum[i],&existing_values[0],P1Enum);
+				xDelete<IssmDouble>(existing_values);
+				xDelete<IssmDouble>(values_to_add);
+			}
+		}
+	}
+}
+/*}}}*/
+void FemModel::AverageSumOutputx(int* stackedinput_enum,int* averagedinput_enum,int numoutputs){ /*{{{*/
+
+	//First get timestep
+	IssmDouble maindt;
+	this->parameters->FindParam(&maindt,TimesteppingTimeStepEnum);
+  for(int i=0;i<numoutputs;i++){
+		if(stackedinput_enum[i]<0){
+			_error_("Can't deal with non enum fields for result Stack");
+		}
+		else{
+			for(int j=0;j<elements->Size();j++){
+				/*Intermediaries*/
+				Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
+				int  numvertices = element->GetNumberOfVertices();
+				IssmDouble* time_averaged=xNew<IssmDouble>(numvertices);
+				IssmDouble* existing_values=xNew<IssmDouble>(numvertices);
+				element->GetInputListOnVertices(&existing_values[0],stackedinput_enum[i]); //those are the values to add
+
+				for(int k=0;k<numvertices;k++){
+					time_averaged[k]=existing_values[k]/maindt;
+				}
+
+				element->AddInput(averagedinput_enum[i],&time_averaged[0],P1Enum);
+				xDelete<IssmDouble>(time_averaged);
+				xDelete<IssmDouble>(existing_values);
 			}
 		}
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 23543)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 23544)
@@ -141,6 +141,6 @@
 		#endif
 		#ifdef _HAVE_ESA_
-		void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy); 
-		void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); 
+		void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy);
+		void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
 		#endif
 		#ifdef _HAVE_SEALEVELRISE_
@@ -148,5 +148,5 @@
 		void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
 		void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
-		void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); 
+		void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
 		IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);
 		#endif
@@ -159,5 +159,8 @@
 		void InitTransientOutputx(int* input_enum, int numoutputs);
 		void StackTransientOutputx(int* input_enum,int* stackedinput_enum, IssmDouble hydrotime, int numoutputs);
-		void AverageTransientOutputx(int* input_enum,int* averagedinput_enum,IssmDouble hydrotime,int numoutputs);
+		void AverageTransientOutputx(int* stackedinput_enum,int* averagedinput_enum,IssmDouble init_time,int numoutputs);
+		void InitMeanOutputx(int* stackedinput_enum,int numoutputs);
+		void SumOutputx(int* input_enum,int* stackedinput_enum,int numoutputs);
+		void AverageSumOutputx(int* stackedinput_enum,int* averagedinput_enum,int numoutputs);
 		void UpdateConstraintsx(void);
 		int  UpdateVertexPositionsx(void);
Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 23543)
+++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 23544)
@@ -75,5 +75,6 @@
 					int stackedinput[4]={EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum};
 					int averagedinput[4]={EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum};
-					femmodel->InitTransientOutputx(&stackedinput[0],4);
+					//femmodel->InitTransientOutputx(&stackedinput[0],4);
+					femmodel->InitMeanOutputx(&stackedinput[0],4);
 					while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
 						hydrostep+=1;
@@ -88,7 +89,9 @@
 						solutionsequence_hydro_nonlinear(femmodel);
 						/*If we have a sub-timestep we stack the variables here*/
-						femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4);
+						//femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4);
+						femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],4);
 					}
-					femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4);
+					//femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4);
+					femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],4);
 				}
 				else{
@@ -96,5 +99,6 @@
 					int stackedinput[2]={EffectivePressureStackedEnum,SedimentHeadStackedEnum};
 					int averagedinput[2]={EffectivePressureEnum,SedimentHeadEnum};
-					femmodel->InitTransientOutputx(&stackedinput[0],2);
+					//femmodel->InitTransientOutputx(&stackedinput[0],2);
+					femmodel->InitMeanOutputx(&stackedinput[0],2);
 					while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
 						hydrotime+=hydrodt;
@@ -104,7 +108,9 @@
 						solutionsequence_hydro_nonlinear(femmodel);
 						/*If we have a sub-timestep we stack the variables here*/
-						femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2);
+						//femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2);
+						femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],2);
 					}
-					femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2);
+					//femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2);
+					femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],2);
 				}
 			}