Index: issm/trunk-jpl/src/m/classes/SMBgemb.py
===================================================================
--- issm/trunk-jpl/src/m/classes/SMBgemb.py	(revision 24213)
+++ issm/trunk-jpl/src/m/classes/SMBgemb.py	(revision 24240)
@@ -31,16 +31,16 @@
 
         #inputs:
-        self.Ta = float('NaN')  #2 m air temperature, in Kelvin
-        self.V = float('NaN')  #wind speed (m / s - 1)
-        self.dswrf = float('NaN')  #downward shortwave radiation flux [W / m^2]
-        self.dlwrf = float('NaN')  #downward longwave radiation flux [W / m^2]
-        self.P = float('NaN')  #precipitation [mm w.e. / m^2]
-        self.eAir = float('NaN')  #screen level vapor pressure [Pa]
-        self.pAir = float('NaN')  #surface pressure [Pa]
-        self.Tmean = float('NaN')  #mean annual temperature [K]
-        self.Vmean = float('NaN')  #mean annual wind velocity [m s - 1]
-        self.C = float('NaN')  #mean annual snow accumulation [kg m - 2 yr - 1]
-        self.Tz = float('NaN')  #height above ground at which temperature (T) was sampled [m]
-        self.Vz = float('NaN')  #height above ground at which wind (V) was sampled [m]
+        self.Ta = float('NaN')       #2 m air temperature, in Kelvin
+        self.V = float('NaN')        #wind speed (m/s-1)
+        self.dswrf = float('NaN')    #downward shortwave radiation flux [W/m^2]
+        self.dlwrf = float('NaN')    #downward longwave radiation flux [W/m^2]
+        self.P = float('NaN')        #precipitation [mm w.e. / m^2]
+        self.eAir = float('NaN')     #screen level vapor pressure [Pa]
+        self.pAir = float('NaN')     #surface pressure [Pa]
+        self.Tmean = float('NaN')    #mean annual temperature [K]
+        self.Vmean = float('NaN')    #mean annual wind velocity [m s-1]
+        self.C = float('NaN')        #mean annual snow accumulation [kg m-2 yr-1]
+        self.Tz = float('NaN')       #height above ground at which temperature (T) was sampled [m]
+        self.Vz = float('NaN')       #height above ground at which wind (V) eas sampled [m]
 
         #optional inputs:
@@ -49,80 +49,63 @@
 
         # Initialization of snow properties
-        self.Dzini = float('NaN')  #cell depth (m)
-        self.Dini = float('NaN')  #snow density (kg m - 3)
-        self.Reini = float('NaN')  #effective grain size (mm)
-        self.Gdnini = float('NaN')  #grain dricity (0 - 1)
-        self.Gspini = float('NaN')  #grain sphericity (0 - 1)
-        self.ECini = float('NaN')  #evaporation/condensation (kg m - 2)
-        self.Wini = float('NaN')  #Water content (kg m - 2)
-        self.Aini = float('NaN')  #albedo (0 - 1)
-        self.Tini = float('NaN')  #snow temperature (K)
+        self.Dzini = float('NaN')    #cell depth (m)
+        self.Dini = float('NaN')     #snow density (kg m-3)
+        self.Reini = float('NaN')    #effective grain size (mm)
+        self.Gdnini = float('NaN')   #grain dricity (0-1)
+        self.Gspini = float('NaN')   #grain sphericity (0-1)
+        self.ECini = float('NaN')    #evaporation/condensation (kg m-2)
+        self.Wini = float('NaN')     #Water content (kg m-2)
+        self.Aini = float('NaN')     #albedo (0-1)
+        self.Tini = float('NaN')     #snow temperature (K)
         self.Sizeini = float('NaN')  #Number of layers
 
         #settings:
-        self.aIdx = float('NaN')  #method for calculating albedo and subsurface absorption (default is 1)
-        # 0: direct input from aValue parameter
-        # 1: effective grain radius [Gardner & Sharp, 2009]
-        # 2: effective grain radius [Brun et al., 2009]
-        # 3: density and cloud amount [Greuell & Konzelmann, 1994]
-        # 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
-        # 5: ingest MODIS mode, direct input from aValue parameter applied to surface ice only
-        self.swIdx = float('NaN')  # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1)
-        self.denIdx = float('NaN')  #densification model to use (default is 2):
-        # 1 = emperical model of Herron and Langway (1980)
-        # 2 = semi - emperical model of Anthern et al. (2010)
-        # 3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010)
-        # 4 = DO NOT USE: emperical model of Li and Zwally (2004)
-        # 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)
-        # 6 = Antarctica semi - emperical model of Ligtenberg et al. (2011)
-        # 7 = Greenland semi - emperical model of Kuipers Munneke et al. (2015)
+        self.aIdx = float('NaN')     #method for calculating albedo and subsurface absorption (default is 1)
+        self.swIdx = float('NaN')    # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1)
+        self.denIdx = float('NaN')   #densification model to use (default is 2):
         self.dsnowIdx = float('NaN')  #model for fresh snow accumulation density (default is 1):
-        # 0 = Original GEMB value, 150 kg / m^3
-        # 1 = Antarctica value of fresh snow density, 350 kg / m^3
-        # 2 = Greenland value of fresh snow density, 315 kg / m^3, Fausto et al. (2008)
-        # 3 = Antarctica model of Kaspers et al. (2004)
-        # 4 = Greenland model of Kuipers Munneke et al. (2015)
-
-        self.zTop = float('NaN')  # depth over which grid length is constant at the top of the snopack (default 10) [m]
-        self.dzTop = float('NaN')  # initial top vertical grid spacing (default .05) [m]
-        self.dzMin = float('NaN')  # initial min vertical allowable grid spacing (default dzTop / 2) [m]
-        self.zY = float('NaN')  # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]
-        self.zMax = float('NaN')  #initial max model depth (default is min(thickness, 250)) [m]
-        self.zMin = float('NaN')  #initial min model depth (default is min(thickness, 130)) [m]
-        self.outputFreq = float('NaN')  #output frequency in days (default is monthly, 30)
+        self.zTop = float('NaN')     # depth over which grid length is constant at the top of the snopack (default 10) [m]
+        self.dzTop = float('NaN')    # initial top vertical grid spacing (default .05) [m]
+        self.dzMin = float('NaN')    # initial min vertical allowable grid spacing (default dzMin/2) [m]
+        self.zY = float('NaN')       # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]
+        self.zMax = float('NaN')     #initial max model depth (default is min(thickness, 500)) [m]
+        self.zMin = float('NaN')     #initial min model depth (default is min(thickness, 30)) [m]
+        self.outputFreq = float('NaN')       #output frequency in days (default is monthly, 30)
 
         #specific albedo parameters:
         #Method 1 and 2:
-        self.aSnow = float('NaN')  # new snow albedo (0.64 - 0.89)
-        self.aIce = float('NaN')  # range 0.27 - 0.58 for old snow
-        #Method 3: Radiation Correction Factors - > only used for met station data and Greuell & Konzelmann, 1994 albedo
+        self.aSnow = float('NaN')    # new snow albedo (0.64 - 0.89)
+        self.aIce = float('NaN')     # range 0.27-0.58 for old snow
+        #Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo
         self.cldFrac = float('NaN')  # average cloud amount
         #Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005)
-        self.t0wet = float('NaN')  # time scale for wet snow (15 - 21.9)
-        self.t0dry = float('NaN')  # warm snow timescale (30)
-        self.K = float('NaN')  # time scale temperature coef. (7)
+        self.t0wet = float('NaN')    # time scale for wet snow (15-21.9)
+        self.t0dry = float('NaN')    # warm snow timescale (30)
+        self.K = float('NaN')        # time scale temperature coef. (7)
         self.adThresh = float('NaN')  # Apply aIdx method to all areas with densities below this value,
         # or else apply direct input value from aValue, allowing albedo to be altered.
-        # Default value is rho water (1023 kg m - 3).
+        # Default value is rho water (1023 kg m-3).
 
         #densities:
-        self.InitDensityScaling = float('NaN')  #initial scaling factor multiplying the density of ice, which describes the density of the snowpack.
+        self.InitDensityScaling = float('NaN')       #initial scaling factor multiplying the density of ice, which describes the density of the snowpack.
 
         #thermo:
         self.ThermoDeltaTScaling = float('NaN')  #scaling factor to multiply the thermal diffusion timestep (delta t)
 
+        self.steps_per_step = 1
         self.requested_outputs = []
 
-    #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM.
-    #dateN: that's the last row of the above fields.
-    #dt:    included in dateN. Not an input.
-    #elev:  this is taken from the ISSM surface itself.
-
-    #}}}
+        #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM.
+        #dateN: that's the last row of the above fields.
+        #dt:    included in dateN. Not an input.
+        #elev:  this is taken from the ISSM surface itself.
+
+        #}}}
+
 
     def __repr__(self):  # {{{
         #string = "   surface forcings parameters:"
-        #string = "  #s\n  #s" % (string, fielddisplay(self, 'mass_balance', 'surface mass balance [m / yr ice eq]'))
-        #string = "  #s\n  #s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
+        #string = "#s\n#s"%(string, fielddisplay(self, 'mass_balance', 'surface mass balance [m/yr ice eq]'))
+        #string = "#s\n#s"%(string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
         string = '   surface forcings for SMB GEMB model :'
         string = "%s\n%s" % (string, fielddisplay(self, 'issmbgradients', 'is smb gradients method activated (0 or 1, default is 0)'))
@@ -137,18 +120,18 @@
         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, 'Ta', '2 m air temperature, in Kelvin'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'V', 'wind speed (m s - 1)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'dlwrf', 'downward shortwave radiation flux [W / m^2]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'dswrf', 'downward longwave radiation flux [W / m^2]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'V', 'wind speed (m s-1)'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'dlwrf', 'downward shortwave radiation flux [W/m^2]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'dswrf', 'downward longwave radiation flux [W/m^2]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'P', 'precipitation [mm w.e. / m^2]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'eAir', 'screen level vapor pressure [Pa]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'pAir', 'surface pressure [Pa]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'Tmean', 'mean annual temperature [K]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'C', 'mean annual snow accumulation [kg m - 2 yr - 1]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'Vmean', 'mean annual wind velocity [m s - 1] (default 10 m / s)'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'C', 'mean annual snow accumulation [kg m-2 yr-1]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'Vmean', 'mean annual temperature [m s-1] (default 10 m/s)'))
         string = "%s\n%s" % (string, fielddisplay(self, 'Tz', 'height above ground at which temperature (T) was sampled [m]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'Vz', 'height above ground at which wind (V) was sampled [m]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'Vz', 'height above ground at which wind (V) eas sampled [m]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'zTop', 'depth over which grid length is constant at the top of the snopack (default 10) [m]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'dzTop', 'initial top vertical grid spacing (default .05) [m] '))
-        string = "%s\n%s" % (string, fielddisplay(self, 'dzMin', 'initial min vertical allowable grid spacing (default dzMin / 2) [m] '))
+        string = "%s\n%s" % (string, fielddisplay(self, 'dzMin', 'initial min vertical allowable grid spacing (default dzMin/2) [m] '))
         string = "%s\n%s" % (string, fielddisplay(self, 'zMax', 'initial max model depth (default is min(thickness, 500)) [m]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'zMin', 'initial min model depth (default is min(thickness, 30)) [m]'))
@@ -167,26 +150,26 @@
         #snow properties init
         string = "%s\n%s" % (string, fielddisplay(self, 'Dzini', 'Initial cell depth when restart [m]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'Dini', 'Initial snow density when restart [kg m - 3]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'Dini', 'Initial snow density when restart [kg m-3]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'Reini', 'Initial grain size when restart [mm]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'Gdnini', 'Initial grain dricity when restart [ - ]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'Gspini', 'Initial grain sphericity when restart [ - ]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'ECini', 'Initial evaporation / condensation when restart [kg m - 2]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'Wini', 'Initial snow water content when restart [kg m - 2]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'Aini', 'Initial albedo when restart [ - ]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'Gdnini', 'Initial grain dricity when restart [-]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'Gspini', 'Initial grain sphericity when restart [-]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'ECini', 'Initial evaporation/condensation when restart [kg m-2]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'Wini', 'Initial snow water content when restart [kg m-2]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'Aini', 'Initial albedo when restart [-]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'Tini', 'Initial snow temperature when restart [K]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'Sizeini', 'Initial number of layers when restart [K]'))
 
         #additional albedo parameters:
-        if isinstance(self.aIdx, list) or isinstance(self.aIdx, np.ndarray) and (self.aIdx == [1, 2] or (1 in self.aIdx and 2 in self.aIdx)):
+        if isinstance(self.aIdx, (list, type(np.array([1, 2])))) and (self.aIdx == [1, 2] or (1 in self.aIdx and 2 in self.aIdx)):
             string = "%s\n%s" % (string, fielddisplay(self, 'aSnow', 'new snow albedo (0.64 - 0.89)'))
-            string = "%s\n%s" % (string, fielddisplay(self, 'aIce', 'albedo of ice (0.27 - 0.58)'))
-        elif elf.aIdx == 0:
+            string = "%s\n%s" % (string, fielddisplay(self, 'aIce', 'albedo of ice (0.27-0.58)'))
+        elif self.aIdx == 0:
             string = "%s\n%s" % (string, fielddisplay(self, 'aValue', 'Albedo forcing at every element.  Used only if aIdx == {0, 5}'))
-        elif elf.aIdx == 5:
+        elif self.aIdx == 5:
             string = "%s\n%s" % (string, fielddisplay(self, 'aValue', 'Albedo forcing at every element.  Used only if aIdx == {0, 5}'))
         elif self.aIdx == 3:
             string = "%s\n%s" % (string, fielddisplay(self, 'cldFrac', 'average cloud amount'))
         elif self.aIdx == 4:
-            string = "%s\n%s" % (string, fielddisplay(self, 't0wet', 'time scale for wet snow (15 - 21.9) [d]'))
+            string = "%s\n%s" % (string, fielddisplay(self, 't0wet', 'time scale for wet snow (15-21.9) [d]'))
             string = "%s\n%s" % (string, fielddisplay(self, 't0dry', 'warm snow timescale (30) [d]'))
             string = "%s\n%s" % (string, fielddisplay(self, 'K', 'time scale temperature coef. (7) [d]'))
@@ -195,16 +178,17 @@
         string = "%s\n%s" % (string, fielddisplay(self, 'denIdx', ['densification model to use (default is 2):',
                                                                    '1 = emperical model of Herron and Langway (1980)',
-                                                                   '2 = semi - emperical model of Anthern et al. (2010)',
+                                                                   '2 = semi-emperical model of Anthern et al. (2010)',
                                                                    '3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010)',
                                                                    '4 = DO NOT USE: emperical model of Li and Zwally (2004)',
                                                                    '5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)',
-                                                                   '6 = Antarctica semi - emperical model of Ligtenberg et al. (2011)',
-                                                                   '7 = Greenland semi - emperical model of Kuipers Munneke et al. (2015)']))
+                                                                   '6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)',
+                                                                   '7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)']))
         string = "%s\n%s" % (string, fielddisplay(self, 'dsnowIdx', ['model for fresh snow accumulation density (default is 1):',
-                                                                     '0 = Original GEMB value, 150 kg / m^3',
-                                                                     '1 = Antarctica value of fresh snow density, 350 kg / m^3',
-                                                                     '2 = Greenland value of fresh snow density, 315 kg / m^3, Fausto et al. (2008)',
+                                                                     '0 = Original GEMB value, 150 kg/m^3',
+                                                                     '1 = Antarctica value of fresh snow density, 350 kg/m^3',
+                                                                     '2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)',
                                                                      '3 = Antarctica model of Kaspers et al. (2004), Make sure to set Vmean accurately',
                                                                      '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, 'requested_outputs', 'additional outputs requested'))
         return string
@@ -247,18 +231,18 @@
         self.denIdx = 2
         self.dsnowIdx = 1
-        self.zTop = 10 * np.ones((mesh.numberofelements, ))
-        self.dzTop = .05 * np.ones((mesh.numberofelements, ))
+        self.zTop = 10 * np.ones((mesh.numberofelements,))
+        self.dzTop = .05 * np.ones((mesh.numberofelements,))
         self.dzMin = self.dzTop / 2
         self.InitDensityScaling = 1.0
         self.ThermoDeltaTScaling = 1 / 11.0
 
-        self.Vmean = 10 * np.ones((mesh.numberofelements, ))
-
-        self.zMax = 250 * np.ones((mesh.numberofelements, ))
-        self.zMin = 130 * np.ones((mesh.numberofelements, ))
-        self.zY = 1.10 * np.ones((mesh.numberofelements, ))
+        self.Vmean = 10 * np.ones((mesh.numberofelements,))
+
+        self.zMax = 250 * np.ones((mesh.numberofelements,))
+        self.zMin = 130 * np.ones((mesh.numberofelements,))
+        self.zY = 1.10 * np.ones((mesh.numberofelements,))
         self.outputFreq = 30
 
-    #additional albedo parameters
+        #additional albedo parameters
         self.aSnow = 0.85
         self.aIce = 0.48
@@ -269,6 +253,6 @@
         self.adThresh = 1023
 
-        self.teValue = np.ones((mesh.numberofelements, ))
-        self.aValue = self.aSnow * np.ones(mesh.numberofelements, )
+        self.teValue = np.ones((mesh.numberofelements,))
+        self.aValue = self.aSnow * np.ones(mesh.numberofelements,)
 
         self.Dzini = 0.05 * np.ones((mesh.numberofelements, 2))
@@ -277,14 +261,15 @@
         self.Gdnini = 0.0 * np.ones((mesh.numberofelements, 2))
         self.Gspini = 0.0 * np.ones((mesh.numberofelements, 2))
-        self.ECini = 0.0 * np.ones((mesh.numberofelements, ))
+        self.ECini = 0.0 * np.ones((mesh.numberofelements,))
         self.Wini = 0.0 * np.ones((mesh.numberofelements, 2))
         self.Aini = self.aSnow * np.ones((mesh.numberofelements, 2))
         self.Tini = 273.15 * np.ones((mesh.numberofelements, 2))
-        #          / !\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh).
-        #         If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp
-        self.Sizeini = 2 * np.ones((mesh.numberofelements, ))
+#       /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh).
+#       If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp
+        self.Sizeini = 2 * np.ones((mesh.numberofelements,))
     #}}}
 
-    def checkconsistency(self, md, solution, analyses):  # {{{
+    def checkconsistency(self, md, solution, analyses):    # {{{
+
         md = checkfield(md, 'fieldname', 'smb.isgraingrowth', 'values', [0, 1])
         md = checkfield(md, 'fieldname', 'smb.isalbedo', 'values', [0, 1])
@@ -297,24 +282,24 @@
         md = checkfield(md, 'fieldname', 'smb.isclimatology', 'values', [0, 1])
 
-        md = checkfield(md, 'fieldname', 'smb.Ta', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100)  # - 100 / 100 celsius min / max value
-        md = checkfield(md, 'fieldname', 'smb.V', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<', 45, 'size', np.shape(self.Ta))  #max 500 km / h
-        md = checkfield(md, 'fieldname', 'smb.dswrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1400, 'size', np.shape(self.Ta))
-        md = checkfield(md, 'fieldname', 'smb.dlwrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, 'size', np.shape(self.Ta))
-        md = checkfield(md, 'fieldname', 'smb.P', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 100, 'size', np.shape(self.Ta))
+        md = checkfield(md, 'fieldname', 'smb.Ta', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100)  #-100/100 celsius min/max value
+        md = checkfield(md, 'fieldname', 'smb.V', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '<', 45, 'size', np.shape(self.Ta))  #max 500 km/h
+        md = checkfield(md, 'fieldname', 'smb.dswrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1400, 'size', np.shape(self.Ta))
+        md = checkfield(md, 'fieldname', 'smb.dlwrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, 'size', np.shape(self.Ta))
+        md = checkfield(md, 'fieldname', 'smb.P', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 100, 'size', np.shape(self.Ta))
         md = checkfield(md, 'fieldname', 'smb.eAir', 'timeseries', 1, 'NaN', 1, 'Inf', 1, 'size', np.shape(self.Ta))
 
         if (self.isclimatology > 0):
-            md = checkfield(md, 'fieldname', 'smb.Ta', 'size', [md.mesh.numberofelements + 1], 'message', 'Ta must have md.mesh.numberofelements + 1 rows in order to force a climatology')
-            md = checkfield(md, 'fieldname', 'smb.V', 'size', [md.mesh.numberofelements + 1], 'message', 'V must have md.mesh.numberofelements + 1 rows in order to force a climatology')
-            md = checkfield(md, 'fieldname', 'smb.dswrf', 'size', [md.mesh.numberofelements + 1], 'message', 'dswrf must have md.mesh.numberofelements + 1 rows in order to force a climatology')
-            md = checkfield(md, 'fieldname', 'smb.dlwrf', 'size', [md.mesh.numberofelements + 1], 'message', 'dlwrf must have md.mesh.numberofelements + 1 rows in order to force a climatology')
-            md = checkfield(md, 'fieldname', 'smb.P', 'size', [md.mesh.numberofelements + 1], 'message', 'P must have md.mesh.numberofelements + 1 rows in order to force a climatology')
-            md = checkfield(md, 'fieldname', 'smb.eAir', 'size', [md.mesh.numberofelements + 1], 'message', 'eAir must have md.mesh.numberofelements + 1 rows in order to force a climatology')
-
-        md = checkfield(md, 'fieldname', 'smb.Tmean', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100)  # - 100 / 100 celsius min / max value
-        md = checkfield(md, 'fieldname', 'smb.C', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>=', 0)
-        md = checkfield(md, 'fieldname', 'smb.Vmean', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>=', 0)
-        md = checkfield(md, 'fieldname', 'smb.Tz', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 5000)
-        md = checkfield(md, 'fieldname', 'smb.Vz', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 5000)
+            md = checkfield(md, 'fieldname', 'smb.Ta', 'size', [md.mesh.numberofelements + 1], 'message', 'Ta must have md.mesh.numberofelements+1 rows in order to force a climatology')
+            md = checkfield(md, 'fieldname', 'smb.V', 'size', [md.mesh.numberofelements + 1], 'message', 'V must have md.mesh.numberofelements+1 rows in order to force a climatology')
+            md = checkfield(md, 'fieldname', 'smb.dswrf', 'size', [md.mesh.numberofelements + 1], 'message', 'dswrf must have md.mesh.numberofelements+1 rows in order to force a climatology')
+            md = checkfield(md, 'fieldname', 'smb.dlwrf', 'size', [md.mesh.numberofelements + 1], 'message', 'dlwrf must have md.mesh.numberofelements+1 rows in order to force a climatology')
+            md = checkfield(md, 'fieldname', 'smb.P', 'size', [md.mesh.numberofelements + 1], 'message', 'P must have md.mesh.numberofelements+1 rows in order to force a climatology')
+            md = checkfield(md, 'fieldname', 'smb.eAir', 'size', [md.mesh.numberofelements + 1], 'message', 'eAir must have md.mesh.numberofelements+1 rows in order to force a climatology')
+
+        md = checkfield(md, 'fieldname', 'smb.Tmean', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100)  #-100/100 celsius min/max value
+        md = checkfield(md, 'fieldname', 'smb.C', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0)
+        md = checkfield(md, 'fieldname', 'smb.Vmean', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0)
+        md = checkfield(md, 'fieldname', 'smb.Tz', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 5000)
+        md = checkfield(md, 'fieldname', 'smb.Vz', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 5000)
 
         md = checkfield(md, 'fieldname', 'smb.teValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
@@ -325,35 +310,36 @@
         md = checkfield(md, 'fieldname', 'smb.dsnowIdx', 'NaN', 1, 'Inf', 1, 'values', [0, 1, 2, 3, 4])
 
-        md = checkfield(md, 'fieldname', 'smb.zTop', 'NaN', 1, 'Inf', 1, '>=', 0)
+        md = checkfield(md, 'fieldname', 'smb.zTop', 'NaN', 1, 'Inf', 1, '> = ', 0)
         md = checkfield(md, 'fieldname', 'smb.dzTop', 'NaN', 1, 'Inf', 1, '>', 0)
         md = checkfield(md, 'fieldname', 'smb.dzMin', 'NaN', 1, 'Inf', 1, '>', 0)
-        md = checkfield(md, 'fieldname', 'smb.zY', 'NaN', 1, 'Inf', 1, '>=', 1)
+        md = checkfield(md, 'fieldname', 'smb.zY', 'NaN', 1, 'Inf', 1, '> = ', 1)
         md = checkfield(md, 'fieldname', 'smb.outputFreq', 'NaN', 1, 'Inf', 1, '>', 0, '<', 10 * 365)  #10 years max
-        md = checkfield(md, 'fieldname', 'smb.InitDensityScaling', 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
-        md = checkfield(md, 'fieldname', 'smb.ThermoDeltaTScaling', 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
+        md = checkfield(md, 'fieldname', 'smb.InitDensityScaling', 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1)
+        md = checkfield(md, 'fieldname', 'smb.ThermoDeltaTScaling', 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1)
         md = checkfield(md, 'fieldname', 'smb.adThresh', 'NaN', 1, 'Inf', 1, '>=', 0)
 
-        if isinstance(self.aIdx, list) or isinstance(self.aIdx, np.ndarray) and (self.aIdx == [1, 2] or (1 in self.aIdx and 2 in self.aIdx)):
-            md = checkfield(md, 'fieldname', 'smb.aSnow', 'NaN', 1, 'Inf', 1, '>=', .64, '<=', .89)
-            md = checkfield(md, 'fieldname', 'smb.aIce', 'NaN', 1, 'Inf', 1, '>=', .27, '<=', .58)
+        if isinstance(self.aIdx, (list, type(np.array([1, 2])))) and (self.aIdx == [1, 2] or (1 in self.aIdx and 2 in self.aIdx)):
+            md = checkfield(md, 'fieldname', 'smb.aSnow', 'NaN', 1, 'Inf', 1, '> = ', .64, '< = ', .89)
+            md = checkfield(md, 'fieldname', 'smb.aIce', 'NaN', 1, 'Inf', 1, '> = ', .27, '< = ', .58)
         elif self.aIdx == 0:
             md = checkfield(md, 'fieldname', 'smb.aValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
         elif self.aIdx == 3:
-            md = checkfield(md, 'fieldname', 'smb.cldFrac', 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
+            md = checkfield(md, 'fieldname', 'smb.cldFrac', 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1)
         elif self.aIdx == 4:
-            md = checkfield(md, 'fieldname', 'smb.t0wet', 'NaN', 1, 'Inf', 1, '>=', 15, '<=', 21.9)
-            md = checkfield(md, 'fieldname', 'smb.t0dry', 'NaN', 1, 'Inf', 1, '>=', 30, '<=', 30)
-            md = checkfield(md, 'fieldname', 'smb.K', 'NaN', 1, 'Inf', 1, '>=', 7, '<=', 7)
+            md = checkfield(md, 'fieldname', 'smb.t0wet', 'NaN', 1, 'Inf', 1, '> = ', 15, '< = ', 21.9)
+            md = checkfield(md, 'fieldname', 'smb.t0dry', 'NaN', 1, 'Inf', 1, '> = ', 30, '< = ', 30)
+            md = checkfield(md, 'fieldname', 'smb.K', 'NaN', 1, 'Inf', 1, '> = ', 7, '< = ', 7)
 
         #check zTop is < local thickness:
         he = np.sum(md.geometry.thickness[md.mesh.elements - 1], axis=1) / np.size(md.mesh.elements, 1)
         if np.any(he < self.zTop):
-            error('SMBgemb consistency check error: zTop should be smaller than local ice thickness')
-
+            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.requested_outputs', 'stringrow', 1)
         return md
     # }}}
 
-    def marshall(self, prefix, md, fid):  # {{{
+    def marshall(self, prefix, md, fid):    # {{{
+
         yts = md.constants.yts
 
@@ -408,5 +394,5 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'teValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
 
-    #snow properties init
+        #snow properties init
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Dzini', 'format', 'DoubleMat', 'mattype', 3)
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Dini', 'format', 'DoubleMat', 'mattype', 3)
@@ -419,6 +405,6 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Tini', 'format', 'DoubleMat', 'mattype', 3)
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Sizeini', 'format', 'IntMat', 'mattype', 2)
-
-    #figure out dt from forcings:
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
+        #figure out dt from forcings:
         time = self.Ta[-1]  #assume all forcings are on the same time step
         dtime = np.diff(time, n=1, axis=0)
@@ -427,9 +413,9 @@
         WriteData(fid, prefix, 'data', dt, 'name', 'md.smb.dt', 'format', 'Double', 'scale', yts)
 
-    # Check if smb_dt goes evenly into transient core time step
+        # Check if smb_dt goes evenly into transient core time step
         if (md.timestepping.time_step % dt >= 1e-10):
-            error('smb_dt/dt =  #f. The number of SMB time steps in one transient core time step has to be an an integer', md.timestepping.time_step / dt)
-
-    #process requested outputs
+            raise IOError('smb_dt/dt = #f. The number of SMB time steps in one transient core time step has to be an an integer', md.timestepping.time_step / dt)
+
+        #process requested outputs
         outputs = self.requested_outputs
         indices = [i for i, x in enumerate(outputs) if x == 'default']
