Index: /issm/trunk-jpl/src/m/classes/SMBgradients.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgradients.m	(revision 27969)
+++ /issm/trunk-jpl/src/m/classes/SMBgradients.m	(revision 27970)
@@ -48,6 +48,6 @@
 				md = checkfield(md,'fieldname','smb.smbref','timeseries',1,'NaN',1,'Inf',1);
 				if max(max(abs(md.smb.smbref(1:end-1,:))))<1
-              disp('!!! Warning: SMBgradients now expects smbref to be in m/yr ice eq. instead of mm/yr water eq.');
-            end
+					disp('!!! Warning: SMBgradients now expects smbref to be in m/yr ice eq. instead of mm/yr water eq.');
+				end
 				md = checkfield(md,'fieldname','smb.b_pos','timeseries',1,'NaN',1,'Inf',1);
 				md = checkfield(md,'fieldname','smb.b_neg','timeseries',1,'NaN',1,'Inf',1);
@@ -61,8 +61,8 @@
 
 			disp(sprintf('\n   SMB gradients parameters:'));
-			fielddisplay(self,'href',' reference elevation from which deviation is used to calculate SMB adjustment in smb gradients method [m]');
-			fielddisplay(self,'smbref',' reference smb from which deviation is calculated in smb gradients method [m/yr ice equiv]');
-			fielddisplay(self,'b_pos',' slope of hs - smb regression line for accumulation regime required if smb gradients is activated');
-			fielddisplay(self,'b_neg',' slope of hs - smb regression line for ablation regime required if smb gradients is activated');
+			fielddisplay(self,'href','reference elevation from which deviation is used to calculate SMB adjustment in smb gradients method [m]');
+			fielddisplay(self,'smbref','reference smb from which deviation is calculated in smb gradients method [m/yr ice equiv]');
+			fielddisplay(self,'b_pos','slope of hs - smb regression line for accumulation regime required if smb gradients is activated');
+			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');
@@ -90,5 +90,5 @@
 			pos  = find(ismember(outputs,'default'));
 			if ~isempty(pos),
-				outputs(pos) = [];                         %remove 'default' from outputs
+				outputs(pos) = [];                                %remove 'default' from outputs
 				outputs      = [outputs defaultoutputs(self,md)]; %add defaults
 			end
Index: /issm/trunk-jpl/src/m/classes/SMBgradients.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBgradients.py	(revision 27969)
+++ /issm/trunk-jpl/src/m/classes/SMBgradients.py	(revision 27970)
@@ -1,2 +1,3 @@
+import numpy as np
 from fielddisplay import fielddisplay
 from checkfield import checkfield
@@ -5,16 +6,15 @@
 
 class SMBgradients(object):
-    """
-    SMBgradients Class definition
+    """SMBgradients Class definition
 
-       Usage:
-          SMBgradients = SMBgradients();
+    Usage:
+        SMBgradients = SMBgradients();
     """
 
     def __init__(self):  # {{{
-        self.href = float('NaN')
-        self.smbref = float('NaN')
-        self.b_pos = float('NaN')
-        self.b_neg = float('NaN')
+        self.href = np.nan
+        self.smbref = np.nan
+        self.b_pos = np.nan
+        self.b_neg = np.nan
         self.steps_per_step = 1
         self.averaging = 0
@@ -26,23 +26,22 @@
 
     def __repr__(self):  # {{{
-        string = "   surface forcings parameters:"
+        s = '   surface forcings parameters:\n'
+        s += '{}\n'.format(fielddisplay(self, 'issmbgradients', 'is smb gradients method activated (0 or 1, default is 0)'))
+        s += '{}\n'.format(fielddisplay(self, 'href', 'reference elevation from which deviation is used to calculate SMB adjustment in smb gradients method'))
+        s += '{}\n'.format(fielddisplay(self, 'smbref', 'reference smb from which deviation is calculated in smb gradients method [m/yr ice equiv]'))
+        s += '{}\n'.format(fielddisplay(self, 'b_pos', 'slope of hs - smb regression line for accumulation regime required if smb gradients is activated'))
+        s += '{}\n'.format(fielddisplay(self, 'b_neg', 'slope of hs - smb regression line for ablation regime required if smb gradients is activated'))
+        s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
+        s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
+        s += '\t\t{}\n'.format('0: Arithmetic (default)')
+        s += '\t\t{}\n'.format('1: Geometric')
+        s += '\t\t{}\n'.format('2: Harmonic')
+        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
 
-        string = "%s\n%s" % (string, fielddisplay(self, 'issmbgradients', 'is smb gradients method activated (0 or 1, default is 0)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'href', ' reference elevation from which deviation is used to calculate SMB adjustment in smb gradients method'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'smbref', ' reference smb from which deviation is calculated in smb gradients method'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'b_pos', ' slope of hs - smb regression line for accumulation regime required if smb gradients is activated'))
-        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'))
-
-        return string
+        return s
     # }}}
 
     def extrude(self, md):  # {{{
-        #Nothing for now
+        # Nothing for now
         return self
     # }}}
@@ -59,12 +58,14 @@
 
     def initialize(self, md):  # {{{
-        #Nothing for now
+        # Nothing for now
         return self
     # }}}
 
-    def checkconsistency(self, md, solution, analyses):    # {{{
+    def checkconsistency(self, md, solution, analyses):  # {{{
         if 'MasstransportAnalysis' in analyses:
             md = checkfield(md, 'fieldname', 'smb.href', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
             md = checkfield(md, 'fieldname', 'smb.smbref', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
+            if np.max(np.max(np.abs(md.smb.smbref[0:-1,]))) < 1:
+                print('!!! Warning: SMBgradients now expects smbref to be in m/yr ice eq. instead of mm/yr water eq.')
             md = checkfield(md, 'fieldname', 'smb.b_pos', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
             md = checkfield(md, 'fieldname', 'smb.b_neg', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
@@ -76,5 +77,5 @@
     # }}}
 
-    def marshall(self, prefix, md, fid):    # {{{
+    def marshall(self, prefix, md, fid):  # {{{
         yts = md.constants.yts
 
@@ -87,5 +88,5 @@
         WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
 
-        #process requested outputs
+        # Process requested outputs
         outputs = self.requested_outputs
         indices = [i for i, x in enumerate(outputs) if x == 'default']
@@ -94,4 +95,3 @@
             outputs = outputscopy
         WriteData(fid, prefix, 'data', outputs, 'name', 'md.smb.requested_outputs', 'format', 'StringArray')
-
     # }}}
Index: /issm/trunk-jpl/test/NightlyRun/test328.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test328.py	(revision 27969)
+++ /issm/trunk-jpl/test/NightlyRun/test328.py	(revision 27970)
@@ -14,9 +14,9 @@
 md = setflowequation(md, 'SSA', 'all')
 md.smb = SMBgradients()
-md.smb.b_pos = -100. + 0.00005 * md.mesh.x - 0.0001 * md.mesh.y
-md.smb.b_neg = 250. + 0.000051 * md.mesh.x - 0.00011 * md.mesh.y
+md.smb.b_pos = (-100. + 0.00005 * md.mesh.x - 0.0001 * md.mesh.y) / 1000. * md.materials.rho_freshwater / md.materials.rho_ice
+md.smb.b_neg = (250. + 0.000051 * md.mesh.x - 0.00011 * md.mesh.y) / 1000. * md.materials.rho_freshwater / md.materials.rho_ice
 md.transient.requested_outputs = ['default', 'TotalSmb']
 md.smb.href = copy.deepcopy(md.geometry.surface)
-md.smb.smbref = 1000. - 0.001 * md.mesh.x - 0.005 * md.mesh.y
+md.smb.smbref= (1000. - 0.001 * md.mesh.x - 0.005 * md.mesh.y) / 1000. * md.materials.rho_freshwater / md.materials.rho_ice
 md.cluster = generic('name', gethostname(), 'np', 3)
 md = solve(md, 'Transient')
Index: /issm/trunk-jpl/test/NightlyRun/test329.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test329.py	(revision 27969)
+++ /issm/trunk-jpl/test/NightlyRun/test329.py	(revision 27970)
@@ -15,8 +15,8 @@
 md = setflowequation(md, 'HO', 'all')
 md.smb = SMBgradients()
-md.smb.b_pos = -100. + 0.00005 * md.mesh.x - 0.0001 * md.mesh.y
-md.smb.b_neg = 250. + 0.000051 * md.mesh.x - 0.00011 * md.mesh.y
+md.smb.b_pos = (-100. + 0.00005 * md.mesh.x - 0.0001 * md.mesh.y) / 1000. * md.materials.rho_freshwater / md.materials.rho_ice
+md.smb.b_neg = (250. + 0.000051 * md.mesh.x - 0.00011 * md.mesh.y) / 1000. * md.materials.rho_freshwater / md.materials.rho_ice
 md.smb.href = copy.deepcopy(md.geometry.surface)
-md.smb.smbref = 1000. - 0.001 * md.mesh.x - 0.005 * md.mesh.y
+md.smb.smbref = (1000. - 0.001 * md.mesh.x - 0.005 * md.mesh.y) / 1000. * md.materials.rho_freshwater / md.materials.rho_ice
 md.transient.requested_outputs = ['default', 'TotalSmb']
 md.cluster = generic('name', gethostname(), 'np', 3)