Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 23665)
+++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 23666)
@@ -85,47 +85,75 @@
 
 			if(hydroslices>1){
-				/*define which variable needs to be averaged on the sub-timestep and initialize as needed*/
-				numaveragedinput = 2;
-				std::vector<int> inputtostack	=	{EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum};
-				std::vector<int> stackedinput	=	{EffectivePressureStackedEnum,SedimentHeadStackedEnum};
-				std::vector<int> averagedinput	=	{EffectivePressureEnum,SedimentHeadEnum};
-
 				if (isefficientlayer){
+					/*define which variable needs to be averaged on the sub-timestep and initialize as needed*/
 					numaveragedinput = 4;
-					inputtostack	=	{EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum};
-					stackedinput	=	{EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum};
-					averagedinput	=	{EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum};
-				}
-
-				femmodel->InitMeanOutputx(&stackedinput[0],numaveragedinput);
-
-				while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
-					hydrostep+=1;
-					hydrotime+=hydrodt;
-					/*Setting substep time as global time*/
-					femmodel->parameters->SetParam(hydrotime,TimeEnum);
-					if(VerboseSolution()) _printf0_("sub iteration " << hydrostep << "/" << hydroslices << "  time [yr]: " << setprecision(4) << hydrotime/yts << " (time step: " << hydrodt/yts << ")\n");
-					if(issmb){
-						if(VerboseSolution()) _printf0_("   computing mass balance\n");
-						SmbAnalysis* analysis = new SmbAnalysis();
-						analysis->Core(femmodel);
-						delete analysis;
+					int inputtostack[4]	=	{EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum};
+					int stackedinput[4]	=	{EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum};
+					int averagedinput[4]	=	{EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum};
+					femmodel->InitMeanOutputx(&stackedinput[0],numaveragedinput);
+
+					while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
+						hydrostep+=1;
+						hydrotime+=hydrodt;
+						/*Setting substep time as global time*/
+						femmodel->parameters->SetParam(hydrotime,TimeEnum);
+						if(VerboseSolution()) _printf0_("sub iteration " << hydrostep << "/" << hydroslices << "  time [yr]: " << setprecision(4) << hydrotime/yts << " (time step: " << hydrodt/yts << ")\n");
+						if(issmb){
+							if(VerboseSolution()) _printf0_("   computing mass balance\n");
+							SmbAnalysis* analysis = new SmbAnalysis();
+							analysis->Core(femmodel);
+							delete analysis;
+						}
+						if(VerboseSolution()) _printf0_("   computing water heads\n");
+						/*save preceding timestep*/
+						InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum);
+						if(isefficientlayer){
+							InputDuplicatex(femmodel,EplHeadHydrostepEnum,EplHeadOldEnum);
+							InputDuplicatex(femmodel,HydrologydcEplThicknessHydrostepEnum,HydrologydcEplThicknessOldEnum);
+						}
+						/*Proceed now to heads computations*/
+						solutionsequence_hydro_nonlinear(femmodel);
+						/*If we have a sub-timestep we stack the variables here*/
+						femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],numaveragedinput);
 					}
-					if(VerboseSolution()) _printf0_("   computing water heads\n");
-					/*save preceding timestep*/
-					InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum);
-					if(isefficientlayer){
-						InputDuplicatex(femmodel,EplHeadHydrostepEnum,EplHeadOldEnum);
-						InputDuplicatex(femmodel,HydrologydcEplThicknessHydrostepEnum,HydrologydcEplThicknessOldEnum);
+					/*Reseting to global time*/
+					femmodel->parameters->SetParam(time,TimeEnum);
+					femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],numaveragedinput);
+				}
+				else{
+					/*define which variable needs to be averaged on the sub-timestep and initialize as needed*/
+					numaveragedinput = 2;
+					int inputtostack[2]	=	{EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum};
+					int stackedinput[2]	=	{EffectivePressureStackedEnum,SedimentHeadStackedEnum};
+					int averagedinput[2]	=	{EffectivePressureEnum,SedimentHeadEnum};
+					femmodel->InitMeanOutputx(&stackedinput[0],numaveragedinput);
+					while(hydrotime<time-(yts*DBL_EPSILON)){ //loop on hydro dts
+						hydrostep+=1;
+						hydrotime+=hydrodt;
+						/*Setting substep time as global time*/
+						femmodel->parameters->SetParam(hydrotime,TimeEnum);
+						if(VerboseSolution()) _printf0_("sub iteration " << hydrostep << "/" << hydroslices << "  time [yr]: " << setprecision(4) << hydrotime/yts << " (time step: " << hydrodt/yts << ")\n");
+						if(issmb){
+							if(VerboseSolution()) _printf0_("   computing mass balance\n");
+							SmbAnalysis* analysis = new SmbAnalysis();
+							analysis->Core(femmodel);
+							delete analysis;
+						}
+						if(VerboseSolution()) _printf0_("   computing water heads\n");
+						/*save preceding timestep*/
+						InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum);
+						if(isefficientlayer){
+							InputDuplicatex(femmodel,EplHeadHydrostepEnum,EplHeadOldEnum);
+							InputDuplicatex(femmodel,HydrologydcEplThicknessHydrostepEnum,HydrologydcEplThicknessOldEnum);
+						}
+						/*Proceed now to heads computations*/
+						solutionsequence_hydro_nonlinear(femmodel);
+						/*If we have a sub-timestep we stack the variables here*/
+						femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],numaveragedinput);
 					}
-					/*Proceed now to heads computations*/
-					solutionsequence_hydro_nonlinear(femmodel);
-					/*If we have a sub-timestep we stack the variables here*/
-					femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],numaveragedinput);
-				}
-				/*Reseting to global time*/
-				femmodel->parameters->SetParam(time,TimeEnum);
-				femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],numaveragedinput);
-			}
+					/*Reseting to global time*/
+					femmodel->parameters->SetParam(time,TimeEnum);
+					femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],numaveragedinput);
+				}
 			else{
 				if(issmb){