Index: /issm/trunk-jpl/src/m/classes/timestepping.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/timestepping.m	(revision 26302)
+++ /issm/trunk-jpl/src/m/classes/timestepping.m	(revision 26303)
@@ -23,4 +23,16 @@
 					error('constructor not supported');
 			end
+		end % }}}
+		function disp(self) % {{{
+			disp(sprintf('   timestepping parameters:'));
+
+			unit = 'yr';
+			fielddisplay(self,'start_time',['simulation starting time [' unit ']']);
+			fielddisplay(self,'final_time',['final time to stop the simulation [' unit ']']);
+			fielddisplay(self,'time_step',['length of time steps [' unit ']']);
+			fielddisplay(self,'interp_forcing','interpolate in time between requested forcing values? (0 or 1)');
+			fielddisplay(self,'cycle_forcing','cycle through forcing? (0 or 1)');
+			fielddisplay(self,'coupling_time',['length of coupling time step with ocean model [' unit ']']);
+
 		end % }}}
 		function self = setdefaultparameters(self) % {{{
@@ -51,16 +63,4 @@
 			end
 		end % }}}
-		function disp(self) % {{{
-			disp(sprintf('   timestepping parameters:'));
-
-			unit = 'yr';
-			fielddisplay(self,'start_time',['simulation starting time [' unit ']']);
-			fielddisplay(self,'final_time',['final time to stop the simulation [' unit ']']);
-			fielddisplay(self,'time_step',['length of time steps [' unit ']']);
-			fielddisplay(self,'interp_forcing','interpolate in time between requested forcing values ? (0 or 1)');
-			fielddisplay(self,'cycle_forcing','cycle through forcing ? (0 or 1)');
-			fielddisplay(self,'coupling_time',['length of coupling time step with ocean model  [' unit ']']);
-
-		end % }}}
 		function marshall(self,prefix,md,fid) % {{{
 
Index: /issm/trunk-jpl/src/m/classes/timestepping.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/timestepping.py	(revision 26302)
+++ /issm/trunk-jpl/src/m/classes/timestepping.py	(revision 26303)
@@ -5,41 +5,44 @@
 
 class timestepping(object):
-    """
-    TIMESTEPPING Class definition
+    """TIMESTEPPING Class definition
 
-       Usage:
-          timestepping = timestepping()
+    Usage:
+        timestepping = timestepping()
     """
 
-    def __init__(self):  # {{{
-        self.start_time = 0.
-        self.final_time = 0.
-        self.time_step = 0.
+    def __init__(self, *args): #{{{
+        self.start_time = 0
+        self.final_time = 0
+        self.time_step = 0
         self.interp_forcing = 1
         self.cycle_forcing = 0
-        self.coupling_time = 0.
+        self.coupling_time = 0
 
-    #set defaults
-        self.setdefaultparameters()
-
+        if len(args) == 0:
+            self.setdefaultparameters()
+        else:
+            raise RuntimeError('constructor not supported')
     #}}}
 
-    def __repr__(self):  # {{{
-        string = "   timestepping parameters:"
-        string = "%s\n%s" % (string, fielddisplay(self, "start_time", "simulation starting time [yr]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "final_time", "final time to stop the simulation [yr]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "time_step", "length of time steps [yr]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "interp_forcing", "interpolate in time between requested forcing values ? (0 or 1)"))
-        string = "%s\n%s" % (string, fielddisplay(self, "cycle_forcing", "cycle through forcing ? (0 or 1)"))
-        string = "%s\n%s" % (string, fielddisplay(self, "coupling_time", "length of coupling time steps with ocean model [yr]"))
-        return string
+    def __repr__(self): #{{{
+        s = '   timestepping parameters:\n'
+        unit = 'yr'
+        s += '{}\n'.format(fielddisplay(self, 'start_time', 'simulation starting time [' + unit + ']'))
+        s += '{}\n'.format(fielddisplay(self, 'final_time', 'final time to stop the simulation [' + unit + ']'))
+        s += '{}\n'.format(fielddisplay(self, 'time_step', 'length of time steps [' + unit + ']'))
+        s += '{}\n'.format(fielddisplay(self, 'interp_forcing', 'interpolate in time between requested forcing values? (0 or 1)'))
+        s += '{}\n'.format(fielddisplay(self, 'cycle_forcing', 'cycle through forcing? (0 or 1)'))
+        s += '{}\n'.format(fielddisplay(self, 'coupling_time', 'length of coupling time steps with ocean model [' + unit + ']'))
+        return s
     #}}}
 
-    def setdefaultparameters(self):  # {{{
-        #time between 2 time steps
-        self.time_step = 1. / 2.
-        #final time
-        self.final_time = 10. * self.time_step
-        #should we interpolate forcing between timesteps?
+    def setdefaultparameters(self): #{{{
+        # Time between 2 time steps
+        self.time_step = 1 / 2
+
+        # Final time
+        self.final_time = 10 * self.time_step
+
+        # Should we interpolate forcing between timesteps?
         self.interp_forcing = 1
         self.cycle_forcing = 0
@@ -48,27 +51,26 @@
     #}}}
 
-    def checkconsistency(self, md, solution, analyses):  # {{{
-
+    def checkconsistency(self, md, solution, analyses): #{{{
         md = checkfield(md, 'fieldname', 'timestepping.start_time', 'numel', [1], 'NaN', 1, 'Inf', 1)
         md = checkfield(md, 'fieldname', 'timestepping.final_time', 'numel', [1], 'NaN', 1, 'Inf', 1)
         md = checkfield(md, 'fieldname', 'timestepping.time_step', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
-        if self.final_time - self.start_time < 0:
-            md.checkmessage("timestepping.final_time should be larger than timestepping.start_time")
-            md = checkfield(md, 'fieldname', 'timestepping.coupling_time', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
-        md = checkfield(md, 'fieldname', 'timestepping.interp_forcing', 'numel', [1], 'values', [0, 1])
         md = checkfield(md, 'fieldname', 'timestepping.cycle_forcing', 'numel', [1], 'values', [0, 1])
+        if (self.final_time - self.start_time) < 0:
+            md.checkmessage('timestepping.final_time should be larger than timestepping.start_time')
+        if solution == 'TransientSolution':
+            md = checkfield(md, 'fieldname', 'timestepping.time_step', 'numel', [1], '>', 0, 'NaN', 1, 'Inf', 1)
+            md = checkfield(md, 'fieldname', 'timestepping.time_step', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
 
         return md
-    # }}}
+    #}}}
 
-    def marshall(self, prefix, md, fid):  # {{{
-
-        yts = md.constants.yts
+    def marshall(self, prefix, md, fid): #{{{
+        scale = md.constants.yts
         WriteData(fid, prefix, 'name', 'md.timestepping.type', 'data', 1, 'format', 'Integer')
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'start_time', 'format', 'Double', 'scale', yts)
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'final_time', 'format', 'Double', 'scale', yts)
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'time_step', 'format', 'Double', 'scale', yts)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'start_time', 'format', 'Double', 'scale', scale)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'final_time', 'format', 'Double', 'scale', scale)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'time_step', 'format', 'Double', 'scale', scale)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'interp_forcing', 'format', 'Boolean')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'cycle_forcing', 'format', 'Boolean')
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'coupling_time', 'format', 'Double', 'scale', yts)
-    # }}}
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'coupling_time', 'format', 'Double', 'scale', scale)
+    #}}}
Index: /issm/trunk-jpl/test/NightlyRun/test2002.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2002.m	(revision 26302)
+++ /issm/trunk-jpl/test/NightlyRun/test2002.m	(revision 26303)
@@ -98,5 +98,4 @@
 Seustatic=md.results.TransientSolution.Sealevel;
 Beustatic=md.results.TransientSolution.Bed;
-pause
 
 %eustatic + selfattraction run:
Index: /issm/trunk-jpl/test/NightlyRun/test2002.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2002.py	(revision 26302)
+++ /issm/trunk-jpl/test/NightlyRun/test2002.py	(revision 26303)
@@ -23,7 +23,7 @@
 
 # Geometry for the bed, arbitrary thickness of 100
-md.geometry.bed = np.zeros((md.mesh.numberofvertices, 1))
+md.geometry.bed = np.zeros((md.mesh.numberofvertices, ))
 md.geometry.base = md.geometry.bed
-md.geometry.thickness = 100 * np.ones((md.mesh.numberofvertices, 1))
+md.geometry.thickness = 100 * np.ones((md.mesh.numberofvertices, ))
 md.geometry.surface = md.geometry.bed + md.geometry.thickness
 
@@ -55,9 +55,9 @@
 # Mask: {{{
 mask = gmtmask(md.mesh.lat, md.mesh.long)
-oceanmask = -1 * np.ones((md.mesh.numberofvertices, 1))
+oceanmask = -1 * np.ones((md.mesh.numberofvertices, ))
 pos = np.where(mask == 0)[0]
 oceanmask[pos] = 1
 
-icemask = np.ones((md.mesh.numberofvertices, 1))
+icemask = np.ones((md.mesh.numberofvertices, ))
 icemask[md.mesh.elements[posant]] = -1
 icemask[md.mesh.elements[posgre]] = -1
@@ -74,9 +74,9 @@
 
 # Masstransport
-md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices, 1))
-md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, 1))
-md.initialization.vx = np.zeros((md.mesh.numberofvertices, 1))
-md.initialization.vy = np.zeros((md.mesh.numberofvertices, 1))
-md.initialization.sealevel = np.zeros((md.mesh.numberofvertices, 1))
+md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices, ))
+md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, ))
+md.initialization.vx = np.zeros((md.mesh.numberofvertices, ))
+md.initialization.vy = np.zeros((md.mesh.numberofvertices, ))
+md.initialization.sealevel = np.zeros((md.mesh.numberofvertices, ))
 md.initialization.str = 0
 
Index: /issm/trunk-jpl/test/NightlyRun/test2003.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2003.m	(revision 26302)
+++ /issm/trunk-jpl/test/NightlyRun/test2003.m	(revision 26303)
@@ -12,10 +12,8 @@
 md.geometry.surface=md.geometry.bed+md.geometry.thickness;
 
-
 %parameterize slc solution:
 %solidearth loading:  {{{
 md.masstransport.spcthickness=[md.geometry.thickness;0];
 md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
-
 
 xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
@@ -44,5 +42,5 @@
 md.mask.ocean_levelset=oceanmask;
 
-% use model representation of ocen area (not the true area)
+% use model representation of ocean area (not the true area)
 md.solidearth.settings.ocean_area_scaling = 0;
 
Index: /issm/trunk-jpl/test/NightlyRun/test2003.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2003.py	(revision 26302)
+++ /issm/trunk-jpl/test/NightlyRun/test2003.py	(revision 26303)
@@ -12,90 +12,115 @@
 
 
-#mesh earth:
+# Mesh earth
+#
+# NOTE: In MATLAB, we currently use cached mesh to account for differences in 
+# mesh generated under Linux versus under macOS
+#
 md = model()
-md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 1000.)  #1000 km resolution mesh
+md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) # 700 km resolution mesh
 
-#parameterize solidearth solution:
+# Geometry for the bed, arbitrary thickness of 100
+md.geometry.bed = -1 * np.ones(md.mesh.numberofvertices)
+md.geometry.base = md.geometry.bed
+md.geometry.thickness = 1000 * np.ones(md.mesh.numberofvertices)
+md.geometry.surface = md.geometry.bed + md.geometry.thickness
+
+# Parameterize slc solution:
 #solidearth loading:  {{{
-md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, ))
-md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, ))
-md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, ))
-md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, ))
-md.dsl.sea_water_pressure_change_at_sea_floor=np.zeros((md.mesh.numberofvertices+1, ))
+md.masstransport.spcthickness = np.append(md.geometry.thickness, 0)
+md.smb.mass_balance = np.zeros(md.mesh.numberofvertices)
 
-#antarctica
-#Access every element in lat using the indices in elements
-# - 1 to convert to base 0 indexing, 1 (not 2, in matlab) to sum over rows
-late = md.mesh.lat[md.mesh.elements - 1].sum(axis=1)/ 3
-longe = md.mesh.long[md.mesh.elements - 1].sum(axis=1)/ 3
-pos = np.intersect1d(np.array(np.where(late < -75)[0]), np.array(np.where(longe < 0)[0]))
-md.solidearth.surfaceload.icethicknesschange[pos] = -1
+xe = md.mesh.x[md.mesh.elements - 1].sum(axis=1) / 3
+ye = md.mesh.y[md.mesh.elements - 1].sum(axis=1) / 3
+ze = md.mesh.z[md.mesh.elements - 1].sum(axis=1) / 3
+re = pow((pow(xe, 2) + pow(ye, 2) + pow(ze, 2)), 0.5)
 
-#elastic loading from love numbers:
-md.solidearth.lovenumbers = lovenumbers('maxdeg', 1000)
+late = asind(ze / re)
+longe = atan2d(ye, xe)
+# Greenland
+pos = np.where(np.logical_and.reduce((late > 60, late < 90, longe > -75, longe < -15)))[0]
+md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 1000
+posice = pos
+
+# Elastic loading from love numbers
+md.solidearth.lovenumbers = lovenumbers('maxdeg', 100)
 #}}}
 
-#mask:  {{{
+# Mask: {{{
 mask = gmtmask(md.mesh.lat, md.mesh.long)
 icemask = np.ones(md.mesh.numberofvertices)
+icemask[md.mesh.elements[posice]] = -1
+md.mask.ice_levelset = icemask
+oceanmask = -1 * np.ones(md.mesh.numberofvertices)
 pos = np.where(mask == 0)[0]
-icemask[pos] = -1
-pos = np.where(mask[md.mesh.elements - 1].sum(axis=1) < 3)[0]
-icemask[md.mesh.elements[pos, :] - 1] = -1
-md.mask.ice_levelset = icemask
-md.mask.ocean_levelset = -icemask
+oceanmask[pos] = 1
+md.mask.ocean_levelset = oceanmask
 
-#make sure that the elements that have loads are fully grounded
-pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0]
-md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1
-
-#make sure wherever there is an ice load, that the mask is set to ice:
-#pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # Do we need to do this twice?
-md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = -1
-# }}}
-
-# use model representation of ocen area (not the true area)
+# Use model representation of ocean area (not the true area)
 md.solidearth.settings.ocean_area_scaling = 0
 
-#geometry
-di = md.materials.rho_ice / md.materials.rho_water
-md.geometry.thickness = np.ones(md.mesh.numberofvertices)
-md.geometry.surface = (1 - di) * np.zeros(md.mesh.numberofvertices)
-md.geometry.base = md.geometry.surface - md.geometry.thickness
-md.geometry.bed = md.geometry.base
+# Materials
+md.initialization.temperature = 273.25 * np.ones(md.mesh.numberofvertices)
+md.initialization.sealevel = np.zeros(md.mesh.numberofvertices)
+md.initialization.str = 0
 
-#materials
-md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, ))
-md.materials.rheology_B = paterson(md.initialization.temperature)
-md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, ))
+md.basalforcings.groundedice_melting_rate = np.zeros(md.mesh.numberofvertices)
+md.basalforcings.floatingice_melting_rate = np.zeros(md.mesh.numberofvertices)
+md.initialization.vx = np.zeros(md.mesh.numberofvertices)
+md.initialization.vy = np.zeros(md.mesh.numberofvertices)
 
-#Miscellaneous
+# Miscellaneous
 md.miscellaneous.name = 'test2003'
 
-#Solution parameters
+# Solution parameters
 md.solidearth.settings.reltol = np.nan
 md.solidearth.settings.abstol = 1e-3
-md.solidearth.settings.computesealevelchange = 1
+md.solidearth.settings.sealevelloading = 0
+md.solidearth.settings.grdocean = 0
+md.solidearth.settings.isgrd = 1
+md.solidearth.settings.ocean_area_scaling = 0
+md.solidearth.settings.grdmodel = 1
+md.solidearth.settings.horiz = 1
+md.solidearth.requested_outputs = ['Sealevel', 'Bed', 'BedEast', 'BedNorth']
 
-#eustatic + rigid + elastic run:
-md.solidearth.settings.rigid = 1
+# Physics
+md.transient.issmb = 0
+md.transient.isstressbalance = 0
+md.transient.isthermal = 0
+md.transient.ismasstransport = 1
+md.transient.isslc = 1
+
+md.timestepping.start_time = 0
+md.timestepping.time_step = 1
+md.timestepping.final_time = 1
+
+# Eustatic + selfattraction + elastic run:
+md.solidearth.settings.selfattraction = 1
 md.solidearth.settings.elastic = 1
 md.solidearth.settings.rotation = 0
+md.solidearth.settings.viscous = 0
 md.cluster = generic('name', gethostname(), 'np', 3)
 #md.verbose = verbose('111111111')
-md = solve(md, 'Sealevelrise')
-SnoRotation = md.results.SealevelriseSolution.Sealevel
+md = solve(md, 'Transient')
+SnoRotation = md.results.TransientSolution.Sealevel
+BUnoRotation = md.results.TransientSolution.Bed
+BEnoRotation = md.results.TransientSolution.BedEast
+BNnoRotation = md.results.TransientSolution.BedNorth
 
-#eustatic + rigid + elastic + rotation run:
-md.solidearth.settings.rigid = 1
+# Eustatic + selfattraction + elastic + rotation run
+md.solidearth.settings.selfattraction = 1
 md.solidearth.settings.elastic = 1
 md.solidearth.settings.rotation = 1
+md.solidearth.settings.viscous = 0
 md.cluster = generic('name', gethostname(), 'np', 3)
 #md.verbose = verbose('111111111')
-md = solve(md, 'Sealevelrise')
-SRotation = md.results.SealevelriseSolution.Sealevel
+md = solve(md, 'Transient')
+SRotation = md.results.TransientSolution.Sealevel
+BURotation = md.results.TransientSolution.Bed
+BERotation = md.results.TransientSolution.BedEast
+BNRotation = md.results.TransientSolution.BedNorth
 
-#Fields and tolerances to track changes
-field_names = ['noRotation', 'Rotation']
-field_tolerances = [1e-13, 1e-13]
-field_values = [SnoRotation, SRotation]
+# Fields and tolerances to track changes
+field_names = ['Sealevel', 'Uplift', 'NorthDisplacement', 'EastDisplacement']
+field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13]
+field_values = [SRotation - SnoRotation, BURotation - BUnoRotation, BNRotation - BNnoRotation,BERotation - BEnoRotation]