Index: /issm/trunk-jpl/src/m/classes/SMBcomponents.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBcomponents.m	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBcomponents.m	(revision 24793)
@@ -11,4 +11,5 @@
 		evaporation = NaN;
 		steps_per_step = 1;
+		averaging = 0;
 		requested_outputs     = {};
 	end
@@ -62,4 +63,5 @@
 			end
 			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
+			md = checkfield(md,'fieldname','smb.averaging',,'numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
 			md = checkfield(md,'fieldname','smb.isclimatology','values',[0 1]);
@@ -80,4 +82,8 @@
 			fielddisplay(self,'isclimatology','repeat all forcings when past last forcing time (default false)');
 			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%55s  0: Arithmetic (default)');
+			disp(sprintf('%55s  0: Geometric');
+			disp(sprintf('%55s  0: Harmonic');
 			fielddisplay(self,'requested_outputs','additional outputs requested');
 		end % }}}
@@ -91,4 +97,5 @@
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','evaporation','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 			WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
 			%process requested outputs
 			outputs = self.requested_outputs;
Index: /issm/trunk-jpl/src/m/classes/SMBcomponents.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBcomponents.py	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBcomponents.py	(revision 24793)
@@ -19,6 +19,7 @@
         self.isclimatology = 0
         self.steps_per_step = 1
+        self.averaging = 0
         self.requested_outputs = []
-        #}}}
+    # }}}
 
     def __repr__(self):  # {{{
@@ -29,7 +30,11 @@
         string = "%s\n%s" % (string, fielddisplay(self, 'isclimatology', 'repeat all forcings when past last forcing time (default false)'))
         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
+        string = "%s\n\t\t%s" % (string, '1: Geometric')
+        string = "%s\n\t\t%s" % (string, '2: Harmonic')
         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
         return string
-        #}}}
+    # }}}
 
     def extrude(self, md):  # {{{
@@ -38,9 +43,9 @@
         self.mass_balance = project3d(md, 'vector', self.evaporation, 'type', 'node')
         return self
-    #}}}
+    # }}}
 
     def defaultoutputs(self, md):  # {{{
         return []
-    #}}}
+    # }}}
 
     def initialize(self, md):  # {{{
@@ -58,5 +63,5 @@
 
         return self
-    #}}}
+    # }}}
 
     def checkconsistency(self, md, solution, analyses):  # {{{
@@ -71,4 +76,5 @@
 
         md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
+        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
         md = checkfield(md, 'fieldname', 'masstransport.requested_outputs', 'stringrow', 1)
         md = checkfield(md, 'fieldname', 'smb.isclimatology', 'values', [0, 1])
@@ -85,4 +91,5 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'evaporation', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
     #process requested outputs
         outputs = self.requested_outputs
Index: /issm/trunk-jpl/src/m/classes/SMBd18opdd.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBd18opdd.m	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBd18opdd.m	(revision 24793)
@@ -30,4 +30,5 @@
 		pddfac_ice                = NaN;
 		steps_per_step            = 1;
+		averaging                 = 0;
 		requested_outputs         = {};
 	end
@@ -115,4 +116,5 @@
 			end
 			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
+			md = checkfield(md,'fieldname','smb.averaging',,'numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
 		end % }}}
@@ -147,4 +149,8 @@
 			end
 			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%55s  0: Arithmetic (default)');
+			disp(sprintf('%55s  0: Geometric');
+			disp(sprintf('%55s  0: Harmonic');
 			fielddisplay(self,'requested_outputs','additional outputs requested');
 			% No need to display rlapslgm, Tdiff, ismungsm
@@ -169,4 +175,5 @@
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','sealev','format','DoubleMat','mattype',1,'timeserieslength',2,'yts',md.constants.yts);
 			WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
 
 			if self.isd18opd
Index: /issm/trunk-jpl/src/m/classes/SMBd18opdd.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBd18opdd.py	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBd18opdd.py	(revision 24793)
@@ -37,4 +37,5 @@
         self.pddfac_ice = float('NaN')
         self.steps_per_step = 1
+        self.averaging = 0
 
     #set defaults
@@ -72,4 +73,9 @@
             string = "%s\n%s" % (string, fielddisplay(self, 'pddfac_ice', 'Pdd factor for ice, at each vertex [mm ice equiv/day/degree C]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
+        string = "%s\n\t\t%s" % (string, '1: Geometric')
+        string = "%s\n\t\t%s" % (string, '2: Harmonic')
+
         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
 
@@ -135,4 +141,5 @@
             md = checkfield(md, 'fieldname', 'smb.rlapslgm', '>=', 0, 'numel', [1])
             md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
+            md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
 
             if self.isd18opd:
@@ -181,4 +188,5 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'sealev', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', 2, 'yts', md.constants.yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
 
         if self.isd18opd:
Index: /issm/trunk-jpl/src/m/classes/SMBforcing.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBforcing.py	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBforcing.py	(revision 24793)
@@ -19,4 +19,5 @@
         self.isclimatology = 0
         self.steps_per_step = 1
+        self.averaging = 0
     #}}}
 
@@ -26,4 +27,8 @@
         string = "%s\n%s" % (string, fielddisplay(self, 'isclimatology', 'repeat all forcings when past last forcing time (default false)'))
         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
+        string = "%s\n\t\t%s" % (string, '1: Geometric')
+        string = "%s\n\t\t%s" % (string, '2: Harmonic')
         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
         return string
@@ -56,4 +61,5 @@
 
         md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
+        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
         md = checkfield(md, 'fieldname', 'masstransport.requested_outputs', 'stringrow', 1)
         md = checkfield(md, 'fieldname', 'smb.isclimatology', 'values', [0, 1])
@@ -69,4 +75,5 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'mass_balance', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
 
         #process requested outputs
Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgemb.m	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBgemb.m	(revision 24793)
@@ -112,4 +112,5 @@
 
 		steps_per_step = 1;
+		averaging = 0;
 		requested_outputs      = {};
 
@@ -286,4 +287,5 @@
 			end
 			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
+			md = checkfield(md,'fieldname','smb.averaging',,'numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
 
@@ -379,4 +381,8 @@
 
 			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%55s  0: Arithmetic (default)');
+			disp(sprintf('%55s  0: Geometric');
+			disp(sprintf('%55s  0: Harmonic');
 			fielddisplay(self,'requested_outputs','additional outputs requested');
 
@@ -450,4 +456,5 @@
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Sizeini','format','IntMat','mattype',2);
 			WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
 			%figure out dt from forcings:
 			if (any(md.smb.P(end,:) - md.smb.Ta(end,:) ~= 0) | ...
Index: /issm/trunk-jpl/src/m/classes/SMBgemb.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgemb.py	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBgemb.py	(revision 24793)
@@ -21,5 +21,6 @@
         #solution choices
         #check these:
-        #isgraingrowth
+        #self.isgraingrowth = 0
+        #self.issmbgradients = 0
         #isalbedo
         #isshortwave
@@ -94,4 +95,5 @@
 
         self.steps_per_step = 1
+        self.averaging = 0
         self.requested_outputs = []
 
@@ -191,4 +193,8 @@
                                                                      '4 = Greenland model of Kuipers Munneke et al. (2015)']))
         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
+        string = "%s\n\t\t%s" % (string, '1: Geometric')
+        string = "%s\n\t\t%s" % (string, '2: Harmonic')
         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
         return string
@@ -218,4 +224,5 @@
     def setdefaultparameters(self, mesh, geometry):  # {{{
         self.isgraingrowth = 1
+        self.issmbgradients = 0
         self.isalbedo = 1
         self.isshortwave = 1
@@ -226,5 +233,5 @@
         self.isturbulentflux = 1
         self.isclimatology = 0
-        self.isconstrainsurfaceT = 0;
+        self.isconstrainsurfaceT = 0
 
         self.aIdx = 1
@@ -274,4 +281,5 @@
 
         md = checkfield(md, 'fieldname', 'smb.isgraingrowth', 'values', [0, 1])
+        #md = checkfield(md, 'fieldname', 'smb.issmbgradients', 'values', [0, 1])
         md = checkfield(md, 'fieldname', 'smb.isalbedo', 'values', [0, 1])
         md = checkfield(md, 'fieldname', 'smb.isshortwave', 'values', [0, 1])
@@ -338,4 +346,5 @@
             raise IOError('SMBgemb consistency check error: zTop should be smaller than local ice thickness')
         md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
+        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
         md = checkfield(md, 'fieldname', 'smb.requested_outputs', 'stringrow', 1)
         return md
@@ -349,4 +358,5 @@
 
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isgraingrowth', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'issmbgradients', 'format', 'Boolean')
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isalbedo', 'format', 'Boolean')
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isshortwave', 'format', 'Boolean')
@@ -409,6 +419,7 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Sizeini', 'format', 'IntMat', 'mattype', 2)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
         #figure out dt from forcings:
-        if (np.any(self.P[-1] - self.Ta[-1] != 0) | np.any(self.V[-1] - self.Ta[-1] != 0) | np.any(self.dswrf[-1] - self.Ta[-1] !=0) | np.any(self.dlwrf[-1] - self.Ta[-1] != 0) | np.any(self.eAir[-1] - self.Ta[-1] != 0) | np.any(self.pAir[-1] - self.Ta[-1] != 0)):
+        if (np.any(self.P[-1] - self.Ta[-1] != 0) | np.any(self.V[-1] - self.Ta[-1] != 0) | np.any(self.dswrf[-1] - self.Ta[-1] != 0) | np.any(self.dlwrf[-1] - self.Ta[-1] != 0) | np.any(self.eAir[-1] - self.Ta[-1] != 0) | np.any(self.pAir[-1] - self.Ta[-1] != 0)):
             raise IOError('All GEMB forcings (Ta, P, V, dswrf, dlwrf, eAir, pAir) must have the same time steps in the final row!')
 
Index: /issm/trunk-jpl/src/m/classes/SMBgradients.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgradients.m	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBgradients.m	(revision 24793)
@@ -11,4 +11,5 @@
 		b_neg             = NaN;
 		steps_per_step    = 1;
+		averaging         = 0;
 		requested_outputs = {};
 	end
@@ -49,4 +50,5 @@
 			end
 			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
+			md = checkfield(md,'fieldname','smb.averaging',,'numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
 		end % }}}
@@ -60,4 +62,8 @@
 			fielddisplay(self,'b_neg',' slope of hs - smb regression line for ablation regime required if smb gradients is activated');
 			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%55s  0: Arithmetic (default)');
+			disp(sprintf('%55s  0: Geometric');
+			disp(sprintf('%55s  0: Harmonic');
 			fielddisplay(self,'requested_outputs','additional outputs requested');
 
@@ -74,4 +80,5 @@
 				  'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 			WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
 
 			%process requested outputs
Index: /issm/trunk-jpl/src/m/classes/SMBgradients.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgradients.py	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBgradients.py	(revision 24793)
@@ -18,4 +18,5 @@
         self.b_neg = float('NaN')
         self.steps_per_step = 1
+        self.averaging = 0
         self.requested_outputs = []
     #}}}
@@ -30,4 +31,8 @@
         string = "%s\n%s" % (string, fielddisplay(self, 'b_neg', ' slope of hs - smb regression line for ablation regime required if smb gradients is activated'))
         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
+        string = "%s\n\t\t%s" % (string, '1: Geometric')
+        string = "%s\n\t\t%s" % (string, '2: Harmonic')
         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
 
@@ -57,4 +62,5 @@
 
         md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
+        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
         md = checkfield(md, 'fieldname', 'masstransport.requested_outputs', 'stringrow', 1)
         return md
@@ -70,4 +76,5 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'b_neg', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
 
         #process requested outputs
Index: /issm/trunk-jpl/src/m/classes/SMBgradientscomponents.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgradientscomponents.m	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBgradientscomponents.m	(revision 24793)
@@ -14,4 +14,5 @@
 		runoffgrad        = NaN;
 		steps_per_step    = 1;
+		averaging         = 0;
 		requested_outputs = {};
 	end
@@ -54,4 +55,5 @@
 			end
 			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
+			md = checkfield(md,'fieldname','smb.averaging',,'numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1);
 		end % }}}
@@ -67,4 +69,8 @@
 			fielddisplay(self,'runoffgrad',' Gradient of the variation of the runoff (0 for uniform runoff) m w.e. m-1 y-1 (lpase rate times ddf)');
 			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%55s  0: Arithmetic (default)');
+			disp(sprintf('%55s  0: Geometric');
+			disp(sprintf('%55s  0: Harmonic');
 			fielddisplay(self,'requested_outputs','additional outputs requested');
 
@@ -81,4 +87,5 @@
 				  'format','Double','scale',1./md.constants.yts);
 			WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
 
 
Index: /issm/trunk-jpl/src/m/classes/SMBgradientscomponents.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgradientscomponents.py	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBgradientscomponents.py	(revision 24793)
@@ -21,4 +21,5 @@
         self.runoffgrad = float('NaN')
         self.steps_per_step = 1
+        self.averaging = 0
         self.requested_outputs = ['default']
     #}}}
@@ -33,4 +34,8 @@
         string = "%s\n%s" % (string, fielddisplay(self, 'runoffgrad', ' Gradient of the variation of the runoff (0 for uniform runoff) m w.e. m-1 y-1 (lapse rate times ddf)'))
         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
+        string = "%s\n\t\t%s" % (string, '1: Geometric')
+        string = "%s\n\t\t%s" % (string, '2: Harmonic')
         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
 
@@ -65,4 +70,5 @@
 
         md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
+        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
         md = checkfield(md, 'fieldname', 'masstransport.requested_outputs', 'stringrow', 1)
         return md
@@ -79,4 +85,5 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'runoffgrad', 'format', 'Double', 'scale', 1. / md.constants.yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
 
         #process requested outputs
Index: /issm/trunk-jpl/src/m/classes/SMBgradientsela.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgradientsela.m	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBgradientsela.m	(revision 24793)
@@ -12,4 +12,5 @@
 		b_min             = NaN;
 		steps_per_step    = 1;
+		averaging         = 0;
 		requested_outputs = {};
 	end
@@ -52,4 +53,5 @@
 			end
 			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
+			md = checkfield(md,'fieldname','smb.averaging',,'numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
 		end % }}}
@@ -64,4 +66,8 @@
 			fielddisplay(self,'b_min',' lower cap on smb rate, default: -9999 (no cap) [m ice eq./yr]');
 			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%55s  0: Arithmetic (default)');
+			disp(sprintf('%55s  0: Geometric');
+			disp(sprintf('%55s  0: Harmonic');
 			fielddisplay(self,'requested_outputs','additional outputs requested');
 
@@ -79,4 +85,5 @@
 				  'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 			WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
 
 			%process requested outputs
Index: /issm/trunk-jpl/src/m/classes/SMBgradientsela.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgradientsela.py	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBgradientsela.py	(revision 24793)
@@ -19,4 +19,5 @@
         self.b_min = float('NaN')
         self.steps_per_step = 1
+        self.averaging = 0
         self.requested_outputs = []
         self.setdefaultparameters()
@@ -33,4 +34,9 @@
         string = "%s\n%s" % (string, fielddisplay(self, 'b_min', ' lower cap on smb rate, default: -9999 (no cap) [m ice eq./yr]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
+        string = "%s\n\t\t%s" % (string, '1: Geometric')
+        string = "%s\n\t\t%s" % (string, '2: Harmonic')
+
         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
         return string
@@ -66,4 +72,5 @@
 
         md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
+        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
         md = checkfield(md, 'fieldname', 'smb.requested_outputs', 'stringrow', 1)
         return md
@@ -80,4 +87,5 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'b_min', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
 
         #process requested outputs
Index: /issm/trunk-jpl/src/m/classes/SMBhenning.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBhenning.m	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBhenning.m	(revision 24793)
@@ -8,4 +8,5 @@
 		smbref = NaN;
 		steps_per_step=1;
+		averaging=0;
 		requested_outputs      = {};
 	end
@@ -53,4 +54,5 @@
 			end
 			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
+			md = checkfield(md,'fieldname','smb.averaging',,'numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
 		end % }}}
@@ -59,4 +61,8 @@
 			fielddisplay(self,'smbref','reference smb from which deviation is calculated [m/yr ice eq]');
 			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%55s  0: Arithmetic (default)');
+			disp(sprintf('%55s  0: Geometric');
+			disp(sprintf('%55s  0: Harmonic');
 			fielddisplay(self,'requested_outputs','additional outputs requested');
 		end % }}}
@@ -68,4 +74,5 @@
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','smbref','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 			WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
 
 			%process requested outputs
Index: /issm/trunk-jpl/src/m/classes/SMBmeltcomponents.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBmeltcomponents.m	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBmeltcomponents.m	(revision 24793)
@@ -12,4 +12,5 @@
 		refreeze = NaN;
 		steps_per_step=1;
+		averaging = 0;
 		requested_outputs      = {};
 	end
@@ -68,4 +69,5 @@
 			end
 			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
+			md = checkfield(md,'fieldname','smb.averaging',,'numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
 			md = checkfield(md,'fieldname','smb.isclimatology','values',[0 1]);
@@ -89,4 +91,8 @@
 			fielddisplay(self,'isclimatology','repeat all forcings when past last forcing time (default false)');
 			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%55s  0: Arithmetic (default)');
+			disp(sprintf('%55s  0: Geometric');
+			disp(sprintf('%55s  0: Harmonic');
 			fielddisplay(self,'requested_outputs','additional outputs requested');
 		end % }}}
@@ -101,4 +107,5 @@
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','refreeze','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 			WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
 
 			%process requested outputs
Index: /issm/trunk-jpl/src/m/classes/SMBmeltcomponents.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBmeltcomponents.py	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBmeltcomponents.py	(revision 24793)
@@ -19,4 +19,5 @@
         self.isclimatology = 0
         self.steps_per_step = 1
+        self.averaging = 0
         self.requested_outputs = []
     #}}}
@@ -30,4 +31,8 @@
         string = "%s\n%s" % (string, fielddisplay(self, 'isclimatology', 'repeat all forcings when past last forcing time (default false)'))
         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
+        string = "%s\n\t\t%s" % (string, '1: Geometric')
+        string = "%s\n\t\t%s" % (string, '2: Harmonic')
         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
         return string
@@ -80,4 +85,5 @@
 
         md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
+        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
         md = checkfield(md, 'fieldname', 'masstransport.requested_outputs', 'stringrow', 1)
         md = checkfield(md, 'fieldname', 'smb.isclimatology', 'values', [0, 1])
@@ -94,4 +100,5 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'refreeze', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
 
         #process requested outputs
Index: /issm/trunk-jpl/src/m/classes/SMBpdd.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBpdd.m	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBpdd.m	(revision 24793)
@@ -26,4 +26,5 @@
 		precipitations_lgm        = NaN;
 		steps_per_step            = 1;
+		averaging                 = 0;
 		requested_outputs         = {};
 	end
@@ -107,4 +108,5 @@
 			end
 			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
+			md = checkfield(md,'fieldname','smb.averaging',,'numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
 
@@ -143,4 +145,8 @@
 			end
 			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%55s  0: Arithmetic (default)');
+			disp(sprintf('%55s  0: Geometric');
+			disp(sprintf('%55s  0: Harmonic');
 			fielddisplay(self,'requested_outputs','additional outputs requested');
 		end % }}}
@@ -160,4 +166,5 @@
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','rlapslgm','format','Double');
 			WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
 			if(self.isdelta18o==0 & self.ismungsm==0)
 				%WriteData(fid,prefix,'object',self,'class','smb','fieldname','monthlytemperatures','format','DoubleMat','mattype',1);
Index: /issm/trunk-jpl/src/m/classes/SMBpdd.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBpdd.py	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBpdd.py	(revision 24793)
@@ -35,4 +35,5 @@
         self.precipitations_lgm = float('NaN')
         self.steps_per_step = 1
+        self.averaging = 0
 
     #set defaults
@@ -74,4 +75,9 @@
 
         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
+        string = "%s\n\t\t%s" % (string, '1: Geometric')
+        string = "%s\n\t\t%s" % (string, '2: Harmonic')
+
         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
 
@@ -155,8 +161,9 @@
                 md = checkfield(md, 'fieldname', 'smb.precipitations_lgm', 'size', [md.mesh.numberofvertices, 12], 'NaN', 1, 'Inf', 1)
                 md = checkfield(md, 'fieldname', 'smb.Pfac', 'NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1)
-                md = checkfield(md, 'fieldname', 'smb.Tdiff','NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1)
-                md = checkfield(md, 'fieldname', 'smb.sealev','NaN',1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1)
+                md = checkfield(md, 'fieldname', 'smb.Tdiff', 'NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1)
+                md = checkfield(md, 'fieldname', 'smb.sealev', 'NaN', 1, 'Inf', 1, 'size', [2, np.nan], 'singletimeseries', 1)
 
         md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
+        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
         md = checkfield(md, 'fieldname', 'masstransport.requested_outputs', 'stringrow', 1)
         return md
@@ -177,4 +184,5 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'rlapslgm', 'format', 'Double')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
 
         if (self.isdelta18o == 0 and self.ismungsm == 0):
Index: /issm/trunk-jpl/src/m/classes/SMBpddSicopolis.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBpddSicopolis.py	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBpddSicopolis.py	(revision 24793)
@@ -28,7 +28,31 @@
         self.isfirnwarming = 0
         self.steps_per_step = 1
+        self.averaging = 0
         self.requested_outputs = []
 
         self.setdefaultparameters()
+    # }}}
+
+    def __repr__(self):  # {{{
+        string = '   surface forcings parameters:'
+        string += '\n   SICOPOLIS PDD scheme (Calov & Greve, 2005) :'
+
+        string = "%s\n%s" % (string, fielddisplay(self, 'monthlytemperatures', 'monthly surface temperatures [K]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'precipitation', 'monthly surface precipitation [m/yr water eq]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'temperature_anomaly', 'anomaly to monthly reference temperature (additive [K])'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'precipitation_anomaly', 'anomaly to monthly precipitation (multiplicative, e.g. q = q0*exp(0.070458*DeltaT) after Huybrechts (2002)) [no unit])'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'smb_corr', 'correction of smb after PDD call [m/a]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 's0p', 'should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 's0t', 'should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'rlaps', 'present day lapse rate (default is 7.4 degree/km)'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'desfac', 'desertification elevation factor (default is -log(2.0)/1000)'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'isfirnwarming', 'is firnwarming (Reeh 1991) activated (0 or 1, default is 1)'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
+        string = "%s\n\t\t%s" % (string, '1: Geometric')
+        string = "%s\n\t\t%s" % (string, '2: Harmonic')
+
+        string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested (TemperaturePDD, SmbAccumulation, SmbMelt)'))
     # }}}
 
@@ -98,28 +122,11 @@
             md = checkfield(md, 'fieldname', 'smb.rlaps', '>=', 0, 'numel', 1)
             md = checkfield(md, 'fieldname', 'smb.monthlytemperatures', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 12])
-            md = checkfield(md, 'fieldname', 'smb.precipitation','NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 12])
+            md = checkfield(md, 'fieldname', 'smb.precipitation', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 12])
 
         md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
+        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
         md = checkfield(md, 'fieldname', 'smb.requested_outputs', 'stringrow', 1)
 
         return md
-    # }}}
-
-    def __repr__(self):  # {{{
-        string = '   surface forcings parameters:'
-        string += '\n   SICOPOLIS PDD scheme (Calov & Greve, 2005) :'
-
-        string = "%s\n%s" % (string, fielddisplay(self, 'monthlytemperatures', 'monthly surface temperatures [K]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'precipitation', 'monthly surface precipitation [m/yr water eq]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'temperature_anomaly', 'anomaly to monthly reference temperature (additive [K])'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'precipitation_anomaly', 'anomaly to monthly precipitation (multiplicative, e.g. q = q0*exp(0.070458*DeltaT) after Huybrechts (2002)) [no unit])'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'smb_corr', 'correction of smb after PDD call [m/a]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 's0p', 'should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 's0t', 'should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'rlaps', 'present day lapse rate (default is 7.4 degree/km)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'desfac', 'desertification elevation factor (default is -log(2.0)/1000)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'isfirnwarming', 'is firnwarming (Reeh 1991) activated (0 or 1, default is 1)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested (TemperaturePDD, SmbAccumulation, SmbMelt)'))
     # }}}
 
@@ -142,4 +149,5 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'smb_corr', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
 
         #process requested outputs
Index: /issm/trunk-jpl/src/m/classes/SMBsemic.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBsemic.m	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/SMBsemic.m	(revision 24793)
@@ -20,4 +20,5 @@
 		s0gcm					= NaN;
 		steps_per_step = 1;
+		averaging = 0;
 		requested_outputs = {};
 	end
@@ -80,4 +81,5 @@
 			end
 			md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
+			md = checkfield(md,'fieldname','smb.averaging',,'numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
 
@@ -104,4 +106,8 @@
 			fielddisplay(self,'s0gcm','GCM reference elevation; (default is 0) [m]');
 			fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%55s  0: Arithmetic (default)');
+			disp(sprintf('%55s  0: Geometric');
+			disp(sprintf('%55s  0: Harmonic');
 			fielddisplay(self,'requested_outputs','additional outputs requested');
 		end % }}}
@@ -124,4 +130,5 @@
 			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailytemperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
 			WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer');
+			WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer');
 			%process requested outputs
 			outputs = self.requested_outputs;
Index: /issm/trunk-jpl/src/m/classes/hydrologydc.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/hydrologydc.m	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/hydrologydc.m	(revision 24793)
@@ -12,5 +12,6 @@
 		rel_tol                  = 0;
 		max_iter                 = 0;
-    steps_per_step           = 0;
+		steps_per_step           = 0;
+		averaging                = 0;
 		sedimentlimit_flag       = 0;
 		sedimentlimit            = 0;
@@ -19,8 +20,8 @@
 		leakage_factor           = 0;
 		basal_moulin_input       = NaN;
-    requested_outputs        = {};
+		requested_outputs        = {};
 
 		spcsediment_head         = NaN;
-    mask_thawed_node         = NaN;
+		mask_thawed_node         = NaN;
 		sediment_transmitivity   = NaN;
 		sediment_compressibility = 0;
@@ -39,20 +40,17 @@
 		epl_conductivity         = 0;
 		eplflip_lock             = 0;
-  end
+	end
 	methods
-    % {{{ function self = extrude(self,md)
-		function self = extrude(self,md)
+		function self = extrude(self,md)    % {{{
 			self.spcsediment_head=project3d(md,'vector',self.spcsediment_head,'type','node','layer',1);
 			self.sediment_transmitivity=project3d(md,'vector',self.sediment_transmitivity,'type','node','layer',1);
 			self.basal_moulin_input=project3d(md,'vector',self.basal_moulin_input,'type','node','layer',1);
-      self.mask_thawed_node=project3d(md,'vector',self.mask_thawed_node,'type','node','layer',1);
+			self.mask_thawed_node=project3d(md,'vector',self.mask_thawed_node,'type','node','layer',1);
 			if(self.isefficientlayer==1);
 				self.spcepl_head=project3d(md,'vector',self.spcepl_head,'type','node','layer',1);
-        self.mask_eplactive_node=project3d(md,'vector',self.mask_eplactive_node,'type','node','layer',1);
-			end
-		end
-    % }}}
-		% {{{ function self = hydrologydc(varargin)
-		function self = hydrologydc(varargin)
+				self.mask_eplactive_node=project3d(md,'vector',self.mask_eplactive_node,'type','node','layer',1);
+			end
+		end    % }}}
+		function self = hydrologydc(varargin) % {{{
 			switch nargin
 				case 0
@@ -61,19 +59,17 @@
 					error('constructor not supported');
 			end
-		end
-		% }}}
-    function list = defaultoutputs(self,md) % {{{
-        list = {'SedimentHead','SedimentHeadResidual','EffectivePressure'};
-        if self.isefficientlayer,
-            list=[list,{'EplHead','HydrologydcMaskEplactiveNode','HydrologydcMaskEplactiveElt','EplHeadSlopeX','EplHeadSlopeY','HydrologydcEplThickness'}];
-        end
-        if self.steps_per_step>1,
-            list = [list,'EffectivePressureSubstep','SedimentHeadSubstep'];
-            if self.isefficientlayer,
-                list = [list,'EplHeadSubstep','HydrologydcEplThicknessSubstep'];
-            end
-        end
-    end % }}}
-
+		end% }}}
+		function list = defaultoutputs(self,md) % {{{
+			list = {'SedimentHead','SedimentHeadResidual','EffectivePressure'};
+			if self.isefficientlayer,
+				list=[list,{'EplHead','HydrologydcMaskEplactiveNode','HydrologydcMaskEplactiveElt','EplHeadSlopeX','EplHeadSlopeY','HydrologydcEplThickness'}];
+			end
+			if self.steps_per_step>1,
+				list = [list,'EffectivePressureSubstep','SedimentHeadSubstep'];
+				if self.isefficientlayer,
+					list = [list,'EplHeadSubstep','HydrologydcEplThicknessSubstep'];
+				end
+			end
+		end % }}}
 		function self = initialize(self,md) % {{{
 			self.epl_colapse_thickness = self.sediment_transmitivity/self.epl_conductivity;
@@ -84,8 +80,5 @@
 
 		end % }}}
-		% {{{ function self = setdefaultparameters(self)
-
-		function self = setdefaultparameters(self)
-
+		function self = setdefaultparameters(self)% {{{
 			%Parameters from de Fleurian 2014
 			self.water_compressibility    = 5.04e-10;
@@ -96,4 +89,5 @@
 			self.max_iter                 = 100;
 			self.steps_per_step           = 1;
+			self.averaging                = 0;
 			self.sedimentlimit_flag       = 0;
 			self.sedimentlimit            = 0;
@@ -101,5 +95,5 @@
 			self.unconfined_flag          = 0;
 			self.leakage_factor           = 1.0e-10;
-      self.requested_outputs        = {'default'};
+			self.requested_outputs        = {'default'};
 
 
@@ -117,11 +111,7 @@
 			self.epl_max_thickness        = 5.0;
 			self.eplflip_lock             = 0;
-		end
-
-		% }}}
-		% {{{ function md = checkconsistency(self,md,solution,analyses)
-
-		function md = checkconsistency(self,md,solution,analyses)
-			%Early return
+		end  % }}}
+		function md = checkconsistency(self,md,solution,analyses)% {{{
+		%Early return
 			if ~ismember('HydrologyDCInefficientAnalysis',analyses) & ~ismember('HydrologyDCEfficientAnalysis',analyses),
 				return;
@@ -135,8 +125,9 @@
 			md = checkfield(md,'fieldname','hydrology.max_iter','>',0,'numel',1);
 			md = checkfield(md,'fieldname','hydrology.steps_per_step','>',0,'numel',1);
+			md = checkfield(md,'fieldname','hydrology.averaging',,'numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]);
 			md = checkfield(md,'fieldname','hydrology.transfer_flag','numel',[1],'values',[0 1]);
 			md = checkfield(md,'fieldname','hydrology.unconfined_flag','numel',[1],'values',[0 1]);
-      md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1);
+			md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1);
 			if self.sedimentlimit_flag==1,
 				md = checkfield(md,'fieldname','hydrology.sedimentlimit','>',0,'numel',1);
@@ -152,5 +143,5 @@
 			md = checkfield(md,'fieldname','hydrology.sediment_thickness','>',0,'numel',1);
 			md = checkfield(md,'fieldname','hydrology.sediment_transmitivity','>=',0,'size',[md.mesh.numberofvertices 1]);
-      md = checkfield(md,'fieldname','hydrology.mask_thawed_node','size',[md.mesh.numberofvertices 1],'values',[0 1]);
+			md = checkfield(md,'fieldname','hydrology.mask_thawed_node','size',[md.mesh.numberofvertices 1],'values',[0 1]);
 
 			if self.isefficientlayer==1,
@@ -169,10 +160,6 @@
 				end
 			end
-		end
-
-		% }}}
-		% {{{ function disp(self)
-
-		function disp(self)
+		end% }}}
+		function disp(self)% {{{
 			disp(sprintf('   hydrology Dual Porous Continuum Equivalent parameters:'));
 			disp(sprintf('   - general parameters'));
@@ -184,4 +171,8 @@
 			fielddisplay(self,'max_iter','maximum number of nonlinear iteration');
 			fielddisplay(self,'steps_per_step','number of hydrology steps per timestep');
+			fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
+			disp(sprintf('%55s  0: Arithmetic (default)');
+			disp(sprintf('%55s  0: Geometric');
+			disp(sprintf('%55s  0: Harmonic');
 			fielddisplay(self,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer');
 			disp(sprintf('%55s  0: no limit',' '));
@@ -224,11 +215,6 @@
 				fielddisplay(self,'eplflip_lock','lock the epl activation to avoid fli-floping (default is 0, no stabilization)');
 			end
-
-		end
-
-		% }}}
-		% {{{ function marshall(self,prefix,md,fid)
-
-		function marshall(self,prefix,md,fid)
+		end % }}}
+		function marshall(self,prefix,md,fid)% {{{
 			WriteData(fid,prefix,'name','md.hydrology.model','data',1,'format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','water_compressibility','format','Double');
@@ -239,4 +225,5 @@
 			WriteData(fid,prefix,'object',self,'fieldname','max_iter','format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer');
+			WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','sedimentlimit_flag','format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','transfer_flag','format','Integer');
@@ -255,6 +242,5 @@
 			WriteData(fid,prefix,'object',self,'fieldname','sediment_thickness','format','Double');
 			WriteData(fid,prefix,'object',self,'fieldname','sediment_transmitivity','format','DoubleMat','mattype',1');
-      WriteData(fid,prefix,'object',self,'fieldname','mask_thawed_node','format','DoubleMat','mattype',1);
-
+			WriteData(fid,prefix,'object',self,'fieldname','mask_thawed_node','format','DoubleMat','mattype',1);
 			if self.isefficientlayer==1,
 				WriteData(fid,prefix,'object',self,'fieldname','spcepl_head','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
@@ -269,14 +255,12 @@
 				WriteData(fid,prefix,'object',self,'fieldname','eplflip_lock','format','Integer');
 			end
-      outputs = self.requested_outputs;
-      pos  = find(ismember(outputs,'default'));
-      if ~isempty(pos),
-        outputs(pos) = [];  %remove 'default' from outputs
-        outputs      = [outputs defaultoutputs(self,md)]; %add defaults
-      end
-      WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
-		end
-
-		% }}}
+			outputs = self.requested_outputs;
+			pos  = find(ismember(outputs,'default'));
+			if ~isempty(pos),
+				outputs(pos) = [];  %remove 'default' from outputs
+				outputs      = [outputs defaultoutputs(self,md)]; %add defaults
+			end
+			WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
+		end% }}}
 	end
 end
Index: /issm/trunk-jpl/src/m/classes/hydrologydc.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/hydrologydc.py	(revision 24792)
+++ /issm/trunk-jpl/src/m/classes/hydrologydc.py	(revision 24793)
@@ -22,4 +22,5 @@
         self.max_iter = 0
         self.steps_per_step = 0
+        self.averaging = 0
         self.sedimentlimit_flag = 0
         self.sedimentlimit = 0
@@ -50,5 +51,5 @@
     #set defaults
         self.setdefaultparameters()
-    #}}}
+    # }}}
 
     def __repr__(self):  # {{{
@@ -62,4 +63,8 @@
         string = "%s\n%s" % (string, fielddisplay(self, 'max_iter', 'maximum number of nonlinear iteration'))
         string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of hydrology steps per time step'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
+        string = "%s\n\t\t%s" % (string, '1: Geometric')
+        string = "%s\n\t\t%s" % (string, '2: Harmonic')
         string = "%s\n%s" % (string, fielddisplay(self, 'basal_moulin_input', 'water flux at a given point [m3 s - 1]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
@@ -105,4 +110,5 @@
             string = "%s\n%s" % (string, fielddisplay(self, 'eplflip_lock', 'lock epl activity to avoid flip - floping (default is 0, no stabilization)'))
         return string
+    # }}}
 
     def extrude(self, md):  # {{{
@@ -115,7 +121,7 @@
             self.mask_eplactive_node = project3d(md, 'vector', self.mask_eplactive_node, 'type', 'node', 'layer', 1)
         return self
-    #}}}
-
-    def setdefaultparameters(self):  #{{{
+    # }}}
+
+    def setdefaultparameters(self):  # {{{
         #Parameters from de Fleurian 2014
         self.water_compressibility = 5.04e-10
@@ -126,4 +132,5 @@
         self.max_iter = 100
         self.steps_per_step = 1
+        self.averaging = 0
         self.sedimentlimit_flag = 0
         self.sedimentlimit = 0
@@ -158,5 +165,5 @@
                 list.extend(['EplHeadSubstep', 'HydrologydcEplThicknessSubstep'])
         return list
-    #}}}
+    # }}}
 
     def initialize(self, md):  # {{{
@@ -169,5 +176,5 @@
     # }}}
 
-    def checkconsistency(self, md, solution, analyses):  #{{{
+    def checkconsistency(self, md, solution, analyses):  # {{{
         #Early return
         if 'HydrologyDCInefficientAnalysis' not in analyses and 'HydrologyDCEfficientAnalysis' not in analyses:
@@ -181,4 +188,5 @@
         md = checkfield(md, 'fieldname', 'hydrology.max_iter', '>', 0., 'numel', [1])
         md = checkfield(md, 'fieldname', 'hydrology.steps_per_step', '>=', 1, 'numel', [1])
+        md = checkfield(md, 'fieldname', 'hydrology.averaging', 'numel', [1], 'values', [0, 1, 2])
         md = checkfield(md, 'fieldname', 'hydrology.sedimentlimit_flag', 'numel', [1], 'values', [0, 1, 2, 3])
         md = checkfield(md, 'fieldname', 'hydrology.transfer_flag', 'numel', [1], 'values', [0, 1])
@@ -214,5 +222,5 @@
     # }}}
 
-    def marshall(self, prefix, md, fid):  #{{{
+    def marshall(self, prefix, md, fid):  # {{{
         WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 1, 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'water_compressibility', 'format', 'Double')
@@ -223,4 +231,5 @@
         WriteData(fid, prefix, 'object', self, 'fieldname', 'max_iter', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'sedimentlimit_flag', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'transfer_flag', 'format', 'Integer')