Index: /issm/trunk-jpl/src/m/classes/dsl.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/dsl.py	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/dsl.py	(revision 25169)
@@ -1,18 +1,19 @@
 import numpy as np
+
+from checkfield import checkfield
 from fielddisplay import fielddisplay
-from checkfield import checkfield
+from project3d import project3d
 from WriteData import WriteData
-from project3d import project3d
 
 
 class dsl(object):
-    """
-    dsl Class definition
+    '''
+    DSL - slass definition
 
-       Usage:
-          dsl = dsl()
-    """
+        Usage:
+            dsl = dsl()
+    '''
 
-    def __init__(self):  # {{{
+    def __init__(self): #{{{
         self.global_average_thermosteric_sea_level_change = 0 #corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)
         self.sea_surface_height_change_above_geoid = float('NaN') #corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)
@@ -21,14 +22,15 @@
     #}}}
 
-    def __repr__(self):  # {{{
-        string = "   dsl parameters:"
-        string = "%s\n%s" % (string, fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'))
-        return string
+    def __repr__(self): #{{{
+        s = "   dsl parameters:"
+        s = "%s\n%s" % (s, fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'))
+
+        return s
     #}}}
 
-    def extrude(self, md):  # {{{
+    def extrude(self, md): #{{{
         self.sea_surface_height_change_above_geoid = project3d(md, 'vector', self.sea_surface_height_change_above_geoid, 'type', 'node')
         self.sea_water_pressure_change_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor, 'type', 'node')
@@ -36,11 +38,11 @@
     #}}}
 
-    def defaultoutputs(self, md):  # {{{
+    def defaultoutputs(self, md): #{{{
         return []
     #}}}
 
-    def checkconsistency(self, md, solution, analyses):  # {{{
+    def checkconsistency(self, md, solution, analyses): #{{{
         # Early retun
-        if not 'SealevelriseAnalysis' in analyses:
+        if 'SealevelriseAnalysis' not in analyses:
             return md
 
@@ -60,5 +62,5 @@
     # }}}
 
-    def marshall(self, prefix, md, fid):    # {{{
+    def marshall(self, prefix, md, fid):   #{{{
         yts = md.constants.yts
         WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer')
Index: /issm/trunk-jpl/src/m/classes/dslmme.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/dslmme.m	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/dslmme.m	(revision 25169)
@@ -15,10 +15,4 @@
 	end
 	methods
-		function self = extrude(self,md) % {{{
-			for i=1:length(self.global_average_thermosteric_sea_level_change),
-				self.sea_surface_height_change_above_geoid{i}=project3d(md,'vector',self.sea_surface_height_change_above_geoid{i},'type','node','layer',1);
-				self.sea_water_pressure_change_at_sea_floor{i}=project3d(md,'vector',self.sea_water_pressure_change_at_sea_floor{i},'type','node','layer',1);
-			end
-		end % }}}
 		function self = dslmme(varargin) % {{{
 			switch nargin
@@ -80,4 +74,10 @@
 
 		end % }}}
+		function self = extrude(self,md) % {{{
+			for i=1:length(self.global_average_thermosteric_sea_level_change),
+				self.sea_surface_height_change_above_geoid{i}=project3d(md,'vector',self.sea_surface_height_change_above_geoid{i},'type','node','layer',1);
+				self.sea_water_pressure_change_at_sea_floor{i}=project3d(md,'vector',self.sea_water_pressure_change_at_sea_floor{i},'type','node','layer',1);
+			end
+		end % }}}
 		function savemodeljs(self,fid,modelname) % {{{
 			
Index: /issm/trunk-jpl/src/m/classes/dslmme.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/dslmme.py	(revision 25169)
+++ /issm/trunk-jpl/src/m/classes/dslmme.py	(revision 25169)
@@ -0,0 +1,80 @@
+import numpy as np
+
+from checkfield import checkfield
+from fielddisplay import fielddisplay
+from project3d import project3d
+from WriteData import WriteData
+
+
+class dslmme(object):
+    '''
+    DSLMME class definition
+
+        Usage:
+            dsl = dslmme() #dynamic sea level class based on a multi-model ensemble of CMIP5 outputs
+    '''
+
+    def __init__(self, *args): #{{{
+        self.modelid                                        = 0 #index into the multi-model ensemble
+        self.global_average_thermosteric_sea_level_change   = [] #corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble.
+        self.sea_surface_height_change_above_geoid          = [] #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble
+        self.sea_water_pressure_change_at_sea_floor         = [] #corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr equivalent, not in Pa/yr!) for each ensemble
+        self.compute_fingerprints                           = 0 #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble
+        
+        nargin = len(args)
+
+        if nargin == 0:
+            self.setdefaultparameters()
+        else:
+            raise Exception('constructor not supported')
+    #}}}
+
+    def __repr__(self): # {{{
+        s = '   gia mme parameters:\n'
+        s += '{}\n'.format(fielddisplay(self, 'modelid', 'index into the multi-model ensemble, determines which field will be used.'))
+        s += '{}\n'.format(fielddisplay(self, 'Ngia', 'geoid (mm/yr).'))
+        s += '{}\n'.format(fielddisplay(self, 'Ugia', 'uplift (mm/yr).'))
+
+        return s
+    #}}}
+
+    def setdefaultparameters(self): # {{{
+        return
+    #}}}
+
+    def checkconsistency(self, md, solution, analyses): # {{{
+        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
+            return md
+
+        for i in range(len(self.global_average_thermosteric_sea_level_change)):
+            md = checkfield(md, 'field', self.global_average_thermosteric_sea_level_change[i], 'NaN', 1, 'Inf', 1)
+            md = checkfield(md, 'field', self.sea_surface_height_change_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
+            md = checkfield(md, 'field', self.sea_water_pressure_change_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
+
+        md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level_change))
+        
+        if self.compute_fingerprints:
+            #check geodetic flag of slr is on
+            if md.slr.geodetic==0,
+                raise Exception('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on')
+
+        return md
+    #}}}
+
+    def marshall(self, prefix, md, fid): #{{{
+        WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 2, 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_fingerprints', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'modelid', 'format', 'Integer')
+        WriteData(fid, prefix, 'name', 'md.dsl.nummodels', 'data', len(self.global_average_thermosteric_sea_level_change), 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level_change', 'format', 'MatArray', 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_change_at_sea_floor', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_change_above_geoid', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
+    #}}}
+
+    def extrude(self, md): #{{{
+        for i in range(len(self.global_average_thermosteric_sea_level_change)):
+            self.sea_surface_height_change_above_geoid[i] = project3d(md, 'vector', self.self.sea_surface_height_change_above_geoid[i], 'type', 'node', 'layer', 1)
+            self.sea_water_pressure_change_at_sea_floor[i] = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor[i], 'type', 'node', 'layer', 1)
+
+        return self
+    #}}}
Index: /issm/trunk-jpl/src/m/classes/geometry.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/geometry.m	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/geometry.m	(revision 25169)
@@ -61,5 +61,5 @@
 				return; 
 			else
-				md = checkfield(md,'fieldname','geometry.surface'  ,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
+				md = checkfield(md,'fieldname','geometry.surface' ,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
 				md = checkfield(md,'fieldname','geometry.base'      ,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
 				md = checkfield(md,'fieldname','geometry.thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
Index: /issm/trunk-jpl/src/m/classes/geometry.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/geometry.py	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/geometry.py	(revision 25169)
@@ -1,40 +1,40 @@
 import numpy as np
+
+from checkfield import checkfield
+from fielddisplay import fielddisplay
 from project3d import project3d
-from fielddisplay import fielddisplay
-from checkfield import checkfield
 from WriteData import WriteData
 
 
 class geometry(object):
-    """
+    '''
     GEOMETRY class definition
 
-       Usage:
-          geometry = geometry()
-    """
+        Usage:
+            geometry = geometry()
+    '''
 
-    def __init__(self):  # {{{
-        self.surface = float('NaN')
-        self.thickness = float('NaN')
-        self.base = float('NaN')
-        self.bed = float('NaN')
-        self.hydrostatic_ratio = float('NaN')
+    def __init__(self): #{{{
+        self.surface = np.nan
+        self.thickness = np.nan
+        self.base = np.nan
+        self.bed = np.nan
+        self.hydrostatic_ratio = np.nan
 
-    #set defaults
+        #set defaults
         self.setdefaultparameters()
-
     #}}}
 
-    def __repr__(self):  # {{{
-
-        string = "   geometry parameters:"
-        string = "%s\n%s" % (string, fielddisplay(self, 'surface', 'ice upper surface elevation [m]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'thickness', 'ice thickness [m]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'base', 'ice base elevation [m]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'bed', 'bed elevation [m]'))
-        return string
+    def __repr__(self): #{{{
+        s = "   geometry parameters:"
+        s = "%s\n%s" % (s, fielddisplay(self, 'surface', 'ice upper surface elevation [m]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'thickness', 'ice thickness [m]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'base', 'ice base elevation [m]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'bed', 'bed elevation [m]'))
+        
+        return s
     #}}}
 
-    def extrude(self, md):  # {{{
+    def extrude(self, md): #{{{
         self.surface = project3d(md, 'vector', self.surface, 'type', 'node')
         self.thickness = project3d(md, 'vector', self.thickness, 'type', 'node')
@@ -45,11 +45,13 @@
     #}}}
 
-    def setdefaultparameters(self):  # {{{
+    def setdefaultparameters(self): #{{{
         return self
     #}}}
 
-    def checkconsistency(self, md, solution, analyses):  # {{{
+    def checkconsistency(self, md, solution, analyses): #{{{
         if (solution == 'TransientSolution' and md.transient.isgia) or (solution == 'GiaSolution'):
-            md = checkfield(md, 'fieldname', 'geometry.thickness', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
+            md = checkfield(md, 'fieldname', 'geometry.thickness', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
+        elif solution == 'SealevelriseSolution':
+            md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
         elif solution == 'LoveSolution':
             return
@@ -72,5 +74,5 @@
     # }}}
 
-    def marshall(self, prefix, md, fid):  # {{{
+    def marshall(self, prefix, md, fid): #{{{
         WriteData(fid, prefix, 'object', self, 'fieldname', 'surface', 'format', 'DoubleMat', 'mattype', 1)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
Index: /issm/trunk-jpl/src/m/classes/giamme.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/giamme.py	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/giamme.py	(revision 25169)
@@ -56,6 +56,6 @@
 
         if isinstance(self.Ngia, list) or self.Ngia.ndim == 1: #list or 1D numpy.ndarray
-            WriteData(fid, prefix, 'data', np.zeros(md.mesh.numberofvertices), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ngia')
-            WriteData(fid, prefix, 'data', np.zeros(md.mesh.numberofvertices), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ugia')
+            WriteData(fid, prefix, 'data', np.zeros((md.mesh.numberofvertices, 1)), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ngia')
+            WriteData(fid, prefix, 'data', np.zeros((md.mesh.numberofvertices, 1)), 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.gia.Ugia')
             WriteData(fid, prefix, 'data', 1, 'name', 'md.gia.modelid', 'format', 'Double')
             WriteData(fid, prefix, 'name', 'md.gia.nummodels', 'data', 1, 'format', 'Integer')
Index: /issm/trunk-jpl/src/m/classes/matdamageice.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/matdamageice.m	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/matdamageice.m	(revision 25169)
@@ -120,5 +120,5 @@
 			md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'});
 			md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
-            
+			
 			if ismember('GiaAnalysis',analyses),
 				md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1);
Index: /issm/trunk-jpl/src/m/classes/matdamageice.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/matdamageice.py	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/matdamageice.py	(revision 25169)
@@ -130,4 +130,14 @@
         md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
         md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
+
+        if 'GiaAnalysis' in analyses:
+            md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1)
+            md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1)
+            md = checkfield(md,'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1)
+            md = checkfield(md,'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1)
+
+        if 'SealevelriseAnalysis' in analyses:
+                md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1)
+
         return md
     # }}}
Index: /issm/trunk-jpl/src/m/classes/matice.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/matice.m	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/matice.m	(revision 25169)
@@ -68,5 +68,5 @@
 
 			%water viscosity (N.s/m^2)
-			self.mu_water=0.001787;  
+			self.mu_water=0.001787;
 
 			%ice heat capacity cp (J/kg/K)
@@ -178,5 +178,4 @@
 			WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2);
 			WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String');
-
 			WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_shear_modulus','format','Double');
 			WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3);
Index: /issm/trunk-jpl/src/m/classes/matice.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/matice.py	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/matice.py	(revision 25169)
@@ -1,2 +1,4 @@
+import numpy as np
+
 from fielddisplay import fielddisplay
 from project3d import project3d
@@ -6,12 +8,12 @@
 
 class matice(object):
-    """
+    '''
     MATICE class definition
 
-       Usage:
-          matice = matice()
-    """
+    Usage:
+            matice = matice()
+    '''
 
-    def __init__(self):  # {{{
+    def __init__(self): #{{{
         self.rho_ice = 0.
         self.rho_water = 0.
@@ -27,9 +29,9 @@
         self.mixed_layer_capacity = 0.
         self.thermal_exchange_velocity = 0.
-        self.rheology_B = float('NaN')
-        self.rheology_n = float('NaN')
+        self.rheology_B = np.nan
+        self.rheology_n = np.nan
         self.rheology_law = ''
 
-        #giaivins:
+        #giaivins
         self.lithosphere_shear_modulus = 0.
         self.lithosphere_density = 0.
@@ -37,37 +39,37 @@
         self.mantle_density = 0.
 
-        #SLR
-        self.earth_density = 5512
+        #slr
+        self.earth_density = 0
         self.setdefaultparameters()
     #}}}
 
-    def __repr__(self):  # {{{
-        string = "   Materials:"
+    def __repr__(self): #{{{
+        s = "   Materials:"
+        s = "%s\n%s" % (s, fielddisplay(self, "rho_ice", "ice density [kg/m^3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "rho_water", "water density [kg/m^3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "rho_freshwater", "fresh water density [kg/m^3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "mu_water", "water viscosity [Ns/m^2]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "heatcapacity", "heat capacity [J/kg/K]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "thermalconductivity", "ice thermal conductivity [W/m/K]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "temperateiceconductivity", "temperate ice thermal conductivity [W/m/K]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "effectiveconductivity_averaging", "computation of effectiveconductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)"))
+        s = "%s\n%s" % (s, fielddisplay(self, "meltingpoint", "melting point of ice at 1atm in K"))
+        s = "%s\n%s" % (s, fielddisplay(self, "latentheat", "latent heat of fusion [J/m^3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "beta", "rate of change of melting point with pressure [K/Pa]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "mixed_layer_capacity", "mixed layer capacity [W/kg/K]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "thermal_exchange_velocity", "thermal exchange velocity [m/s]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
+        s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', or 'NyeH2O'"))
+        s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]"))
+        s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]"))
 
-        string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg / m^3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "rho_water", "water density [kg / m^3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "rho_freshwater", "fresh water density [kg / m^3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "mu_water", "water viscosity [N s / m^2]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "heatcapacity", "heat capacity [J / kg / K]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "thermalconductivity", "ice thermal conductivity [W / m / K]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "temperateiceconductivity", "temperate ice thermal conductivity [W / m / K]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "effectiveconductivity_averaging", "computation of effectiveconductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)"))
-        string = "%s\n%s" % (string, fielddisplay(self, "meltingpoint", "melting point of ice at 1atm in K"))
-        string = "%s\n%s" % (string, fielddisplay(self, "latentheat", "latent heat of fusion [J / m^3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "beta", "rate of change of melting point with pressure [K / Pa]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "mixed_layer_capacity", "mixed layer capacity [W / kg / K]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "thermal_exchange_velocity", "thermal exchange velocity [m / s]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1 / n)]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
-        string = "%s\n%s" % (string, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', or 'NyeH2O'"))
-        string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_density", "Lithosphere density [g / cm^ - 3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "mantle_density", "Mantle density [g / cm^ - 3]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "earth_density", "Mantle density [kg / m^ - 3]"))
-        return string
+        return s
     #}}}
 
-    def extrude(self, md):  # {{{
+    def extrude(self, md): #{{{
         self.rheology_B = project3d(md, 'vector', self.rheology_B, 'type', 'node')
         self.rheology_n = project3d(md, 'vector', self.rheology_n, 'type', 'element')
@@ -75,30 +77,30 @@
     #}}}
 
-    def setdefaultparameters(self):  # {{{
-        #ice density (kg / m^3)
+    def setdefaultparameters(self): #{{{
+        #ice density (kg/m^3)
         self.rho_ice = 917.
-        #ocean water density (kg / m^3)
+        #ocean water density (kg/m^3)
         self.rho_water = 1023.
-        #fresh water density (kg / m^3)
+        #fresh water density (kg/m^3)
         self.rho_freshwater = 1000.
-        #water viscosity (N.s / m^2)
+        #water viscosity (N.s/m^2)
         self.mu_water = 0.001787
-        #ice heat capacity cp (J / kg / K)
+        #ice heat capacity cp (J/kg/K)
         self.heatcapacity = 2093.
-        #ice latent heat of fusion L (J / kg)
+        #ice latent heat of fusion L (J/kg)
         self.latentheat = 3.34e5
-        #ice thermal conductivity (W / m / K)
+        #ice thermal conductivity (W/m/K)
         self.thermalconductivity = 2.4
+        #temperate ice thermal conductivity (W/m/K)
+        self.temperateiceconductivity = 0.24
         #computation of effective conductivity
         self.effectiveconductivity_averaging = 1
-        #temperate ice thermal conductivity (W / m / K)
-        self.temperateiceconductivity = 0.24
         #the melting point of ice at 1 atmosphere of pressure in K
         self.meltingpoint = 273.15
-        #rate of change of melting point with pressure (K / Pa)
+        #rate of change of melting point with pressure (K/Pa)
         self.beta = 9.8e-8
-        #mixed layer (ice-water interface) heat capacity (J / kg / K)
+        #mixed layer (ice-water interface) heat capacity (J/kg/K)
         self.mixed_layer_capacity = 3974.
-        #thermal exchange velocity (ice-water interface) (m / s)
+        #thermal exchange velocity (ice-water interface) (m/s)
         self.thermal_exchange_velocity = 1.00e-4
         #Rheology law: what is the temperature dependence of B with T
@@ -107,32 +109,37 @@
 
         # GIA:
-        self.lithosphere_shear_modulus = 6.7e10  # (Pa)
-        self.lithosphere_density = 3.32  # (g / cm^ - 3)
+        self.lithosphere_shear_modulus = 6.7e10 # (Pa)
+        self.lithosphere_density = 3.32  # (g/cm^-3)
         self.mantle_shear_modulus = 1.45e11  # (Pa)
-        self.mantle_density = 3.34  # (g / cm^ - 3)
+        self.mantle_density = 3.34  # (g/cm^-3)
 
-        #SLR
-        self.earth_density = 5512  # average density of the Earth, (kg / m^3)
-        return self
+        # SLR
+        self.earth_density = 5512  # average density of the Earth, (kg/m^3)
     #}}}
 
-    def checkconsistency(self, md, solution, analyses):  # {{{
-        md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
-        md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0)
-        md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0)
-        md = checkfield(md, 'fieldname', 'materials.mu_water', '>', 0)
-        md = checkfield(md, 'fieldname', 'materials.rheology_B', '>', 0, 'timeseries', 1, 'NaN', 1, 'Inf', 1)
-        md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements])
-        md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O'])
-        md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
-        md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', [1])
-        md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', [1])
-        md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', [1])
-        md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
-        md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
+    def checkconsistency(self, md, solution, analyses): #{{{
+        if solution != 'SealevelriseSolution':
+            md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
+            md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0)
+            md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0)
+            md = checkfield(md, 'fieldname', 'materials.mu_water', '>', 0)
+            md = checkfield(md, 'fieldname', 'materials.rheology_B', '>', 0, 'timeseries', 1, 'NaN', 1, 'Inf', 1)
+            md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements])
+            md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O'])
+            md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
+
+        if 'GiaAnalysis' in analyses:
+            md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', [1])
+            md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', [1])
+            md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', [1])
+            md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
+
+        if 'SealevelriseAnalysis' in analyses:
+            md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
+
         return md
-    # }}}
+    #}}}
 
-    def marshall(self, prefix, md, fid):  # {{{
+    def marshall(self, prefix, md, fid): #{{{
         WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 3, 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')
@@ -151,6 +158,5 @@
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_B', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
-        WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
-
+        WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 's')
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10.**3.)
@@ -158,4 +164,3 @@
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10.**3.)
         WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
-
-    # }}}
+    #}}}
Index: /issm/trunk-jpl/src/m/classes/nodalvalue.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/nodalvalue.py	(revision 25169)
+++ /issm/trunk-jpl/src/m/classes/nodalvalue.py	(revision 25169)
@@ -0,0 +1,69 @@
+import numpy as np
+
+from checkfield import checkfield
+from fielddisplay import fielddisplay
+from WriteData import WriteData
+
+
+class dslmme(object):
+    '''
+    NODALVALUE class definition
+
+        Usage:
+            nodalvalue=nodalvalue()
+            nodalvalue=nodalvalue(
+                'name', 'SealevelriseSNodalValue',
+                'definitionstring', 'Outputdefinition1',
+                'model_string', 'SealevelriseS',
+                'node', 1
+                )
+    ''' 
+
+    def __init__(self, *args): #{{{
+        self.name               = ''
+        self.definitionstring   = '' #string that identifies this output definition uniquely, from 'Outputdefinition[1-10]'
+        self.model_string       = '' #string for field that is being retrieved
+        self.node               = np.nan #for which node are we retrieving the value?
+        
+        #use provided options to change fields
+        options = pairoptions(*args)
+
+        #get name
+        self.name               = options.getfieldvalue(options, 'name', '')
+        self.definitionstring   = options.getfieldvalue(options, 'definitionstring', '')
+        self.model_string       = options.getfieldvalue(options, 'model_string', '')
+        self.node               = options.getfieldvalue(options, 'node', '')
+    #}}}
+
+    def __repr__(self): # {{{
+        s = '   Nodalvalue:\n'
+        s += '{}\n'.format(fielddisplay(self, 'name', 'identifier for this nodalvalue response'))
+        s += '{}\n'.format(fielddisplay(self, 'definitionstring', 'string that identifies this output definition uniquely, from \'Outputdefinition[1-10]\''))
+        s += '{}\n'.format(fielddisplay(self, 'model_string', 'string for field that is being retrieved'))
+        s += '{}\n'.format(fielddisplay(self, 'node', 'vertex index at which we retrieve the value'))
+
+        return s
+    #}}}
+
+    def setdefaultparameters(self): # {{{
+        return
+    #}}}
+
+    def checkconsistency(self, md, solution, analyses): # {{{
+        if not isinstance(self.name, basestring):
+            raise Exception('nodalvalue error message: \'name\' field should be a string!')
+        OutputdefinitionStringArray = []
+        for i in range(100):
+            OutputdefinitionStringArray[i] = 'Outputdefinition{}'.format(i)
+        md = checkfield(md, 'fieldname', 'self.definitionstring', 'field', self.definitionstring, 'values', OutputdefinitionStringArray)
+        md = checkfield(md, 'fieldname', 'self.node', 'field', self.node, 'values', range(md.mesh.numberofvertices))
+
+        return md
+    #}}}
+
+    def marshall(self, prefix, md, fid): #{{{
+        WriteData(fid, prefix, 'data', self.name, 'name', 'md.nodalvalue.name', 'format', 'String')
+        WriteData(fid, prefix, 'data', self.definitionstring, 'name', 'md.nodalvalue.definitionenum', 'format', 'String')
+        WriteData(fid, prefix, 'data', self.model_string, 'name', 'md.nodalvalue.model_enum', 'format', 'String')
+        WriteData(fid, prefix, 'data', self.node, 'name', 'md.nodalvalue.node', 'format', 'Integer')
+    #}}}
Index: /issm/trunk-jpl/src/m/classes/slr.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/slr.py	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/slr.py	(revision 25169)
@@ -50,32 +50,32 @@
 
     def __repr__(self):  # {{{
-        string = '   slr parameters:'
-        string = "%s\n%s" % (string, fielddisplay(self, 'deltathickness', 'thickness change: ice height equivalent [m]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'sealevel', 'current sea level (prior to computation) [m]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'spcthickness', 'thickness constraints (NaN means no constraint) [m]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'reltol', 'sea level rise relative convergence criterion, (NaN: not applied)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'abstol', 'sea level rise absolute convergence criterion, (default, NaN: not applied)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'love_h', 'load Love number for radial displacement'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'love_k', 'load Love number for gravitational potential perturbation'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'love_l', 'load Love number for horizontal displaements'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'tide_love_k', 'tidal load Love number (degree 2)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'tide_love_h', 'tidal load Love number (degree 2)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'fluid_love', 'secular fluid Love number'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'equatorial_moi', 'mean equatorial moment of inertia [kg m^2]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'polar_moi', 'polar moment of inertia [kg m^2]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'angular_velocity', 'mean rotational velocity of earth [per second]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'ocean_area_scaling', 'correction for model representation of ocean area [default: No correction]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'hydro_rate', 'rate of hydrological expansion [mm / yr]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'geodetic', 'compute geodetic SLR? (in addition to steric?) default 0'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'geodetic_run_frequency', 'how many time steps we skip before we run SLR solver during transient (default: 1)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'rigid', 'rigid earth graviational potential perturbation'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'rotation', 'earth rotational potential perturbation'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green''s functions'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'transitions', 'indices into parts of the mesh that will be icecaps'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
-
-        return string
+        s = '   slr parameters:'
+        s = "%s\n%s" % (s, fielddisplay(self, 'deltathickness', 'thickness change: ice height equivalent [m]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'sealevel', 'current sea level (prior to computation) [m]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'spcthickness', 'thickness constraints (NaN means no constraint) [m]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'reltol', 'sea level rise relative convergence criterion, (NaN: not applied)'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'abstol', 'sea level rise absolute convergence criterion, (default, NaN: not applied)'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'maxiter', 'maximum number of nonlinear iterations'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'love_h', 'load Love number for radial displacement'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'love_k', 'load Love number for gravitational potential perturbation'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'love_l', 'load Love number for horizontal displaements'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'tide_love_k', 'tidal load Love number (degree 2)'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'tide_love_h', 'tidal load Love number (degree 2)'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'fluid_love', 'secular fluid Love number'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'equatorial_moi', 'mean equatorial moment of inertia [kg m^2]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'polar_moi', 'polar moment of inertia [kg m^2]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'angular_velocity', 'mean rotational velocity of earth [per second]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'ocean_area_scaling', 'correction for model representation of ocean area [default: No correction]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'hydro_rate', 'rate of hydrological expansion [mm / yr]'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'geodetic', 'compute geodetic SLR? (in addition to steric?) default 0'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'geodetic_run_frequency', 'how many time steps we skip before we run SLR solver during transient (default: 1)'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'rigid', 'rigid earth graviational potential perturbation'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'rotation', 'earth rotational potential perturbation'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green''s functions'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'transitions', 'indices into parts of the mesh that will be icecaps'))
+        s = "%s\n%s" % (s, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
+
+        return s
     # }}}
 
@@ -143,5 +143,5 @@
         md = checkfield(md, 'fieldname', 'slr.hydro_rate', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
         md = checkfield(md, 'fieldname', 'slr.degacc', 'size', [1, 1], '>=', 1e-10)
-        md = checkfield(md, 'fieldname', 'slr.requested_outputs', 'stringrow', 1)
+        md = checkfield(md, 'fieldname', 'slr.requested_outputs', 'srow', 1)
         md = checkfield(md, 'fieldname', 'slr.horiz', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
 
@@ -202,4 +202,4 @@
             outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
             outputs = outputscopy
-        WriteData(fid, prefix, 'data', outputs, 'name', 'md.slr.requested_outputs', 'format', 'StringArray')
-    # }}}
+        WriteData(fid, prefix, 'data', outputs, 'name', 'md.slr.requested_outputs', 'format', 'sArray')
+    # }}}
Index: /issm/trunk-jpl/src/m/classes/solidearth.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/solidearth.py	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/solidearth.py	(revision 25169)
@@ -24,5 +24,5 @@
         self.settings = solidearthsettings()
         self.surfaceload = surfaceload()
-        self.love = lovenumbers()
+        self.lovenumbers = lovenumbers()
         self.rotational = rotational()
         self.planetradius = planetradius('earth')
@@ -65,5 +65,5 @@
         self.settings.checkconsistency(md, solution, analyses)
         self.surfaceload.checkconsistency(md, solution, analyses)
-        self.love.checkconsistency(md, solution, analyses)
+        self.lovenumbers.checkconsistency(md, solution, analyses)
         self.rotational.checkconsistency(md, solution, analyses)
 
@@ -82,5 +82,5 @@
         self.settings.marshall(prefix, md, fid)
         self.surfaceload.marshall(prefix, md, fid)
-        self.love.marshall(prefix, md, fid)
+        self.lovenumbers.marshall(prefix, md, fid)
         self.rotational.marshall(prefix, md, fid)
 
Index: /issm/trunk-jpl/src/m/classes/toolkits.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/toolkits.m	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/toolkits.m	(revision 25169)
@@ -5,169 +5,175 @@
 
 classdef toolkits < dynamicprops
-    properties (SetAccess=public) 
-		 DefaultAnalysis  = struct();
-		 RecoveryAnalysis = struct();
-		 %The other properties are dynamic
-	 end
-	 methods (Static)
-		 function self = loadobj(self) % {{{
-			 % This function is directly called by matlab when a model object is
-			 % loaded. Update old properties here
+	properties (SetAccess=public) 
+		DefaultAnalysis  = struct();
+		RecoveryAnalysis = struct();
+		%The other properties are dynamic
+	end
+	methods (Static)
+		function self = loadobj(self) % {{{
+			% This function is directly called by matlab when a model object is
+			% loaded. Update old properties here
 
-			 if isempty(fieldnames(self.RecoveryAnalysis));
-				 disp('WARNING: updating toolkits (RecoveryAnalysis not set)');
-				 self.RecoveryAnalysis  = self.DefaultAnalysis;
-			 end
-		 end% }}}
-	 end
-	 methods
-		 function self = toolkits(varargin) % {{{
-			 switch nargin
-				 case 0
-					 self=setdefaultparameters(self);
-				 case 1
-					 self=structtoobj(self,varargin{1});
-				 otherwise
-					 error('constructor not supported');
-				 end
-			 end % }}}
-		 function self = addoptions(self,analysis,varargin) % {{{
-		 % Usage example:
-		 %    md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis',FSoptions());
-		 %    md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis');
+			if isempty(fieldnames(self.RecoveryAnalysis));
+				disp('WARNING: updating toolkits (RecoveryAnalysis not set)');
+				self.RecoveryAnalysis  = self.DefaultAnalysis;
+			end
+		end% }}}
+	end
+	methods
+		function self = toolkits(varargin) % {{{
+			switch nargin
+				case 0
+					self=setdefaultparameters(self);
+				case 1
+					self=structtoobj(self,varargin{1});
+				otherwise
+					error('constructor not supported');
+				end
+			end % }}}
+		function self = addoptions(self,analysis,varargin) % {{{
+		%ADDOPTIONS - add analysis to md.toolkits.analysis
+		%
+		%   Optional third parameter adds toolkits options to analysis.
+		%
+		%   Usage:
+		%      md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis',FSoptions());
+		%      md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis');
 
-			 %Create dynamic property if property does not exist yet
-			 if ~ismember(analysis,properties(self)),
-				 self.addprop(analysis);
-			 end
+			%Create dynamic property if property does not exist yet
+			if ~ismember(analysis,properties(self)),
+				self.addprop(analysis);
+			end
 
-			 %Add toolkits options to analysis
-			 if nargin==3, self.(analysis) = varargin{1}; end
-		 end
-		 %}}}
-		 function self = setdefaultparameters(self) % {{{
+			%Add toolkits options to analysis
+			if nargin==3,
+				self.(analysis) = varargin{1};
+			end
+		end
+		%}}}
+		function self = setdefaultparameters(self) % {{{
 
-			 %default toolkits: 
-			 if IssmConfig('_HAVE_PETSC_'),
-				 %MUMPS is the default toolkits
-				 if IssmConfig('_HAVE_MUMPS_'),
-					 self.DefaultAnalysis           = mumpsoptions();
-				 else
-					 self.DefaultAnalysis           = iluasmoptions(); 
-				 end
-			 else
-				 if IssmConfig('_HAVE_MUMPS_'),
-					 self.DefaultAnalysis           = issmmumpssolver(); 
-				 elseif IssmConfig('_HAVE_GSL_'),
-					 self.DefaultAnalysis           = issmgslsolver(); 
-				 else 
-					 disp('WARNING: Need at least Mumps or Gsl to define an issm solver type, no default solver assigned');
-				 end
-			 end
+			%default toolkits: 
+			if IssmConfig('_HAVE_PETSC_'),
+				%MUMPS is the default toolkits
+				if IssmConfig('_HAVE_MUMPS_'),
+					self.DefaultAnalysis           = mumpsoptions();
+				else
+					self.DefaultAnalysis           = iluasmoptions(); 
+				end
+			else
+				if IssmConfig('_HAVE_MUMPS_'),
+					self.DefaultAnalysis           = issmmumpssolver(); 
+				elseif IssmConfig('_HAVE_GSL_'),
+					self.DefaultAnalysis           = issmgslsolver(); 
+				else 
+					disp('WARNING: Need at least Mumps or Gsl to define an issm solver type, no default solver assigned');
+				end
+			end
 
-			 %Use same solver for Recovery mode
-			 self.RecoveryAnalysis = self.DefaultAnalysis;
+			%Use same solver for Recovery mode
+			self.RecoveryAnalysis = self.DefaultAnalysis;
 
 
-		 end % }}}
-		 function disp(self) % {{{
-			 analyses=properties(self);
-			 disp(sprintf('List of toolkits options per analysis:\n'));
-			 for i=1:numel(analyses),
-				 analysis=analyses{i};
-				 disp([analysis ':']);
-				 disp(self.(analysis));
-			 end
-		 end % }}}
-		 function md = checkconsistency(self,md,solution,analyses) % {{{
-			 analyses=properties(self);
-			 for i=1:numel(analyses),
-				 switch analyses{i}
-					 case 'DefaultAnalysis'
-					 case 'RecoveryAnalysis'
-					 case 'StressbalanceAnalysis'
-					 case 'GLheightadvectionAnalysis'
-					 case 'MasstransportAnalysis'
-					 case 'ThermalAnalysis'
-					 case 'EnthalpyAnalysis'
-					 case 'AdjointBalancethicknessAnalysis'
-					 case 'BalancethicknessAnalysis'
-					 case 'Balancethickness2Analysis'
-					 case 'BalancethicknessSoftAnalysis'
-					 case 'BalancevelocityAnalysis'
-					 case 'DamageEvolutionAnalysis'
-					 case 'LoveAnalysis'
-					 case 'EsaAnalysis'
-					 case 'SealevelriseAnalysis'
-					 otherwise
+		end % }}}
+		function disp(self) % {{{
+			analyses=properties(self);
+			disp(sprintf('List of toolkits options per analysis:\n'));
+			for i=1:numel(analyses),
+				analysis=analyses{i};
+				disp([analysis ':']);
+				disp(self.(analysis));
+			end
+		end % }}}
+		function md = checkconsistency(self,md,solution,analyses) % {{{
+			analyses=properties(self);
+			for i=1:numel(analyses),
+				switch analyses{i}
+					case 'DefaultAnalysis'
+					case 'RecoveryAnalysis'
+					case 'StressbalanceAnalysis'
+					case 'GLheightadvectionAnalysis'
+					case 'MasstransportAnalysis'
+					case 'ThermalAnalysis'
+					case 'EnthalpyAnalysis'
+					case 'AdjointBalancethicknessAnalysis'
+					case 'BalancethicknessAnalysis'
+					case 'Balancethickness2Analysis'
+					case 'BalancethicknessSoftAnalysis'
+					case 'BalancevelocityAnalysis'
+					case 'DamageEvolutionAnalysis'
+					case 'LoveAnalysis'
+					case 'EsaAnalysis'
+					case 'SealevelriseAnalysis'
+					otherwise
 						md = checkmessage(md,['md.toolkits.' analyses{i} ' not supported yet']);
-				 end
-				 if isempty(fieldnames(self.(analyses{i})))
-					 md = checkmessage(md,['md.toolkits.' analyses{i} ' is empty']);
-				 end
-			 end
-		 end % }}}
-		 function ToolkitsFile(toolkits,filename) % {{{
-			 %TOOLKITSFILE - build toolkits file
-			 %
-			 %   Build a Petsc compatible options file, from the toolkits model field  + return options string. 
-			 %   This file will also be used when the toolkit used is 'issm' instead of 'petsc'
-			 %
-			 %   Usage:     ToolkitsFile(toolkits,filename);
+				end
+				if isempty(fieldnames(self.(analyses{i})))
+					md = checkmessage(md,['md.toolkits.' analyses{i} ' is empty']);
+				end
+			end
+		end % }}}
+		function ToolkitsFile(toolkits,filename) % {{{
+			%TOOLKITSFILE - build toolkits file
+			%
+			%   Build a Petsc compatible options file, from the toolkits model field  + return options string. 
+			%   This file will also be used when the toolkit used is 'issm' instead of 'petsc'
+			%
+			%   Usage:     ToolkitsFile(toolkits,filename);
 
-			 %open file for writing
-			 fid=fopen(filename,'w');
-			 if fid==-1,
-				 error(['ToolkitsFile error: could not open ' filename ' for writing']);
-			 end
+			%open file for writing
+			fid=fopen(filename,'w');
+			if fid==-1,
+				error(['ToolkitsFile error: could not open ' filename ' for writing']);
+			end
 
-			 %write header
-			 fprintf(fid,'%s%s%s\n','%Toolkits options file: ',filename,' written from Matlab toolkits array');
+			%write header
+			fprintf(fid,'%s%s%s\n','%Toolkits options file: ',filename,' written from Matlab toolkits array');
 
-			 %start writing options
-			 analyses=properties(toolkits);
-			 for i=1:numel(analyses),
-				 analysis=analyses{i};
-				 options=toolkits.(analysis);
+			%start writing options
+			analyses=properties(toolkits);
+			for i=1:numel(analyses),
+				analysis=analyses{i};
+				options=toolkits.(analysis);
 
-				 %first write analysis:
-				 fprintf(fid,'\n+%s\n',analysis); %append a + to recognize it's an analysis enum
+				%first write analysis:
+				fprintf(fid,'\n+%s\n',analysis); %append a + to recognize it's an analysis enum
 
-				 %now, write options
-				 optionslist=fieldnames(options);
-				 for j=1:numel(optionslist),
-					 optionname=optionslist{j};
-					 optionvalue=options.(optionname);
+				%now, write options
+				optionslist=fieldnames(options);
+				for j=1:numel(optionslist),
+					optionname=optionslist{j};
+					optionvalue=options.(optionname);
 
-					 if isempty(optionvalue),
-						 %this option has only one argument
-						 fprintf(fid,'-%s\n',optionname);
-					 else
-						 %option with value. value can be string or scalar
-						 if isnumeric(optionvalue),
-							 fprintf(fid,'-%s %g\n',optionname,optionvalue);
-						 elseif ischar(optionvalue),
-							 fprintf(fid,'-%s %s\n',optionname,optionvalue);
-						 else
-							 error(['ToolkitsFile error: option ' optionname ' is not well formatted']);
-						 end
-					 end
-				 end
-			 end
+					if isempty(optionvalue),
+						%this option has only one argument
+						fprintf(fid,'-%s\n',optionname);
+					else
+						%option with value. value can be string or scalar
+						if isnumeric(optionvalue),
+							fprintf(fid,'-%s %g\n',optionname,optionvalue);
+						elseif ischar(optionvalue),
+							fprintf(fid,'-%s %s\n',optionname,optionvalue);
+						else
+							error(['ToolkitsFile error: option ' optionname ' is not well formatted']);
+						end
+					end
+				end
+			end
 
-			 fclose(fid);
-		 end %}}}
-		 function savemodeljs(self,fid,modelname) % {{{
+			fclose(fid);
+		end %}}}
+		function savemodeljs(self,fid,modelname) % {{{
 
-			 analyses=properties(self);
-			 for i=1:numel(analyses),
-				 if isempty(fieldnames(self.(analyses{i})))
-					 error(['md.toolkits.' analyses{i} ' is empty']);
-				 else
-					 writejsstruct(fid,[modelname '.toolkits.' analyses{i}],self.(analyses{i}));
-				 end
-			 end
+			analyses=properties(self);
+			for i=1:numel(analyses),
+				if isempty(fieldnames(self.(analyses{i})))
+					error(['md.toolkits.' analyses{i} ' is empty']);
+				else
+					writejsstruct(fid,[modelname '.toolkits.' analyses{i}],self.(analyses{i}));
+				end
+			end
 
-		 end % }}}
-	 end
+		end % }}}
+	end
  end
Index: /issm/trunk-jpl/src/m/classes/toolkits.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/toolkits.py	(revision 25168)
+++ /issm/trunk-jpl/src/m/classes/toolkits.py	(revision 25169)
@@ -8,12 +8,12 @@
 
 class toolkits(object):
-    """
+    '''
     TOOLKITS class definition
 
-       Usage:
-          self = toolkits()
-    """
+        Usage:
+            self = toolkits()
+    '''
 
-    def __init__(self):  # {{{
+    def __init__(self): #{{{
         #default toolkits
         if IssmConfig('_HAVE_PETSC_')[0]:
@@ -34,16 +34,16 @@
         self.RecoveryAnalysis = self.DefaultAnalysis
 
-    #The other properties are dynamic
-    # }}}
+        #The other properties are dynamic
+    #}}}
 
-    def __repr__(self):  # {{{
+    def __repr__(self): #{{{
         s = "List of toolkits options per analysis:\n\n"
         for analysis in list(vars(self).keys()):
-            s += "%s\n" % fielddisplay(self, analysis, '')
+            s += "{}\n".format(fielddisplay(self, analysis, ''))
 
-            return s
-    # }}}
+        return s
+    #}}}
 
-    def addoptions(self, analysis, *args):  # {{{
+    def addoptions(self, analysis, *args): #{{{
         # Usage example:
         #    md.toolkits = addoptions(md.toolkits, 'StressbalanceAnalysis', FSoptions())
@@ -59,26 +59,32 @@
 
         return self
-    # }}}
+    #}}}
 
-    def checkconsistency(self, md, solution, analyses):  # {{{
+    def checkconsistency(self, md, solution, analyses): #{{{
+        # TODO
+        # - Implement something closer to a switch as in 
+        # src/m/classes/toolkits.m?
+        #
         for analysis in list(vars(self).keys()):
             if not getattr(self, analysis):
-                md.checkmessage("md.toolkits.%s is empty" % analysis)
+                md.checkmessage("md.toolkits.{} is empty".format(analysis))
 
         return md
-    # }}}
+    #}}}
 
-    def ToolkitsFile(self, filename):  # {{{
-        """
+    def ToolkitsFile(self, filename): #{{{
+        '''
         TOOLKITSFILE - build toolkits file
 
-           Build a Petsc compatible options file, from the toolkits model field + return options string
-           This file will also be used when the toolkit used is 'issm' instead of 'petsc'
+            Build a Petsc compatible options file, from the toolkits model 
+            field + return options string.
+            This file will also be used when the toolkit used is 'issm' instead 
+            of 'petsc'.s
 
+            Usage:
+                ToolkitsFile(toolkits, filename)
+        '''
 
-           Usage:     ToolkitsFile(toolkits, filename)
-        """
-
-    #open file for writing
+        #open file for writing
         try:
             fid = open(filename, 'w')
@@ -86,14 +92,14 @@
             raise IOError("ToolkitsFile error: could not open {}' for writing due to".format(filename), e)
 
-    #write header
+        #write header
         fid.write("%s%s%s\n" % ('%Toolkits options file: ', filename, ' written from Python toolkits array'))
 
-    #start writing options
+        #start writing options
         for analysis in list(vars(self).keys()):
             options = getattr(self, analysis)
 
-    #first write analysis:
+            #first write analysis:
             fid.write("\n+{}\n".format(analysis))  #append a + to recognize it's an analysis enum
-    #now, write options
+            #now, write options
             for optionname, optionvalue in list(options.items()):
 
@@ -111,3 +117,3 @@
 
         fid.close()
-    # }}}
+    #}}}
Index: /issm/trunk-jpl/src/m/consistency/checkfield.py
===================================================================
--- /issm/trunk-jpl/src/m/consistency/checkfield.py	(revision 25168)
+++ /issm/trunk-jpl/src/m/consistency/checkfield.py	(revision 25169)
@@ -8,30 +8,33 @@
 
 def checkfield(md, *args):
-    """
+    '''
     CHECKFIELD - check field consistency
 
-       Used to check model consistency.,
-       Requires:
-       'field' or 'fieldname' option. If 'fieldname' is provided, it will retrieve it from the model md. (md.(fieldname))
-             If 'field' is provided, it will assume the argument following 'field' is a numeric array.
-
-       Available options:
- - NaN: 1 if check that there is no NaN
- - size: [lines cols], NaN for non checked dimensions, or 'universal' for any input type (nodal, element, time series, etc)
- -> :  greater than provided value
- ->= : greater or equal to provided value
- - < :  smallerthan provided value
- - <=: smaller or equal to provided value
- - < vec:  smallerthan provided values on each vertex
- - timeseries: 1 if check time series consistency (size and time)
- - values: cell of strings or vector of acceptable values
- - numel: list of acceptable number of elements
- - cell: 1 if check that is cell
- - empty: 1 if check that non empty
- - message: overloaded error message
-
-       Usage:
-          md = checkfield(md, fieldname, options)
-    """
+        Used to check model consistency
+        
+        Requires:
+            'field' or 'fieldname' option. If 'fieldname' is provided, it will retrieve it from the model md. (md.(fieldname))
+            If 'field' is provided, it will assume the argument following 
+            'field' is a numeric array.
+
+        Available options:
+        - NaN: 1 if check that there is no NaN
+        - size: [lines cols], NaN for non checked dimensions, or 'universal' 
+        for any input type (nodal, element, time series, etc)
+        - > :  greater than provided value
+        - >= : greater or equal to provided value
+        - < :  smallerthan provided value
+        - <=: smaller or equal to provided value
+        - < vec:  smallerthan provided values on each vertex
+        - timeseries: 1 if check time series consistency (size and time)
+        - values: list of acceptable values
+        - numel: list of acceptable number of elements
+        - cell: 1 if check that is cell
+        - empty: 1 if check that non empty
+        - message: overloaded error message
+
+        Usage:
+            md = checkfield(md, fieldname, options)
+    '''
 
     #get options
@@ -143,5 +146,5 @@
     if options.getfieldvalue('cell', 0):
         if not isinstance(field, (tuple, list, dict)):
-            md = md.checkmessage(options.getfieldvalue('message', "field '{}' should be a cell".format(fieldname)))
+            md = md.checkmessage(options.getfieldvalue('message', "field '{}' should be a tuple, list, or dict".format(fieldname)))
 
     #check values
Index: /issm/trunk-jpl/src/m/solve/solve.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/solve.py	(revision 25168)
+++ /issm/trunk-jpl/src/m/solve/solve.py	(revision 25169)
@@ -2,11 +2,11 @@
 import os
 import shutil
+
+from ismodelselfconsistent import ismodelselfconsistent
+from loadresultsfromcluster import loadresultsfromcluster
+from marshall import marshall
 from pairoptions import pairoptions
-from ismodelselfconsistent import ismodelselfconsistent
-from marshall import marshall
+from preqmu import *
 from waitonlock import waitonlock
-from loadresultsfromcluster import loadresultsfromcluster
-from preqmu import *
-#from MatlabFuncs import *
 
 
@@ -15,33 +15,35 @@
     SOLVE - apply solution sequence for this model
 
-       Usage:
-          md = solve(md, solutionstring, varargin)
-          where varargin is a list of paired arguments of string OR enums
+        Usage:
+            md = solve(md, solutionstring, varargin)
+            where varargin is a list of paired arguments of string OR enums
 
         solution types available comprise:
-         - 'Stressbalance'    or 'sb'
-         - 'Masstransport'    or 'mt'
-         - 'Thermal'          or 'th'
-         - 'Steadystate'      or 'ss'
-         - 'Transient'        or 'tr'
-         - 'Balancethickness' or 'mc'
-         - 'Balancevelocity'  or 'bv'
-         - 'BedSlope'         or 'bsl'
-         - 'SurfaceSlope'     or 'ssl'
-         - 'Hydrology'        or 'hy'
-         - 'DamageEvolution'  or 'da'
-         - 'Gia'              or 'gia'
-         - 'Esa'          or 'esa'
-         - 'Sealevelrise'     or 'slr'
-         - 'Love'             or 'lv'
+            - 'Stressbalance'      or 'sb'
+            - 'Masstransport'      or 'mt'
+            - 'Thermal'            or 'th'
+            - 'Steadystate'        or 'ss'
+            - 'Transient'          or 'tr'
+            - 'Balancethickness'   or 'mc'
+            - 'Balancevelocity'    or 'bv'
+            - 'BedSlope'           or 'bsl'
+            - 'SurfaceSlope'       or 'ssl'
+            - 'Hydrology'          or 'hy'
+            - 'DamageEvolution'    or 'da'
+            - 'Gia'                or 'gia'
+            - 'Esa'                or 'esa'
+            - 'Sealevelrise'       or 'slr'
+            - 'Love'               or 'lv'
 
-       extra options:
- - loadonly : does not solve. only load results
-         - checkconsistency : 'yes' or 'no' (default is 'yes'), ensures checks on consistency of model
-         - restart: 'directory name (relative to the execution directory) where the restart file is located.
+        extra options:
+            - loadonly          : does not solve. only load results
+            - checkconsistency  : 'yes' or 'no' (default is 'yes'), ensures 
+                                  checks on consistency of model
+            - restart           : 'directory name (relative to the execution 
+                                  directory) where the restart file is located
 
-       Examples:
-          md = solve(md, 'Stressbalance')
-         md = solve(md, 'sb')
+        Examples:
+            md = solve(md, 'Stressbalance')
+            md = solve(md, 'sb')
     '''
 
Index: /issm/trunk-jpl/test/NightlyRun/test2002.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2002.m	(revision 25168)
+++ /issm/trunk-jpl/test/NightlyRun/test2002.m	(revision 25169)
@@ -65,5 +65,7 @@
 
 %eustatic run:
-md.solidearth.settings.rigid=0; md.solidearth.settings.elastic=0;md.solidearth.settings.rotation=0;
+md.solidearth.settings.rigid=0;
+md.solidearth.settings.elastic=0;
+md.solidearth.settings.rotation=0;
 md=solve(md,'Sealevelrise');
 Seustatic=md.results.SealevelriseSolution.Sealevel;
Index: /issm/trunk-jpl/test/NightlyRun/test2002.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2002.py	(revision 25168)
+++ /issm/trunk-jpl/test/NightlyRun/test2002.py	(revision 25169)
@@ -39,5 +39,5 @@
 #mask:  {{{
 mask = gmtmask(md.mesh.lat, md.mesh.long)
-icemask = np.ones(md.mesh.numberofvertices)
+icemask = np.ones((md.mesh.numberofvertices, 1))
 pos = np.where(mask == 0)[0]
 icemask[pos] = -1
@@ -48,6 +48,6 @@
 
 #make sure that the ice level set is all inclusive:
-md.mask.land_levelset = np.zeros((md.mesh.numberofvertices))
-md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices))
+md.mask.land_levelset = np.zeros((md.mesh.numberofvertices, 1))
+md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, 1))
 
 #make sure that the elements that have loads are fully grounded