Index: /issm/trunk-jpl/src/m/classes/friction.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/friction.m	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/friction.m	(revision 25499)
@@ -79,11 +79,11 @@
 		end % }}}
 		function marshall(self,prefix,md,fid) % {{{
-			yts=md.constants.yts;
-
 			WriteData(fid,prefix,'name','md.friction.law','data',1,'format','Integer');
 			if(size(self.coefficient,1)==md.mesh.numberofvertices | size(self.coefficient,1)==md.mesh.numberofvertices+1),
-				mattype=1; tsl = md.mesh.numberofvertices;
+				mattype=1;
+				tsl = md.mesh.numberofvertices;
 			else
-				mattype=2; tsl = md.mesh.numberofelements;
+				mattype=2;
+				tsl = md.mesh.numberofelements;
 			end
 			WriteData(fid,prefix,'object',self,'fieldname','coefficient','format','DoubleMat','mattype',mattype,'timeserieslength',tsl+1,'yts',md.constants.yts);
Index: /issm/trunk-jpl/src/m/classes/friction.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/friction.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/friction.py	(revision 25499)
@@ -1,2 +1,4 @@
+import numpy as np
+
 from checkfield import checkfield
 from fielddisplay import fielddisplay
@@ -6,34 +8,38 @@
 
 class friction(object):
-    '''
-    FRICTION class definition
+    """FRICTION class definition
 
-        Usage:
-            friction = friction()
-    '''
+    Usage:
+        friction = friction()
+    """
 
     def __init__(self):  # {{{
-        self.coefficient = float('NaN')
-        self.p = float('NaN')
-        self.q = float('NaN')
+        self.coefficient = np.nan
+        self.p = np.nan
+        self.q = np.nan
         self.coupling = 0
-        self.effective_pressure = float('NaN')
+        self.effective_pressure = np.nan
         self.effective_pressure_limit = 0
+
         #set defaults
         self.setdefaultparameters()
-        #self.requested_outputs = []
     #}}}
 
     def __repr__(self):  # {{{
-        string = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b, \n(effective stress Neff = rho_ice * g * thickness + rho_water * g * base, r = q / p and s = 1 / p)"
+        s = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b, \n(effective stress Neff = rho_ice * g * thickness + rho_water * g * base, r = q / p and s = 1 / p)\n"
+        s = "{}\n".format(fielddisplay(self, "coefficient", "friction coefficient [SI]"))
+        s = "{}\n".format(fielddisplay(self, "p", "p exponent"))
+        s = "{}\n".format(fielddisplay(self, "q", "q exponent"))
+        s = "{}\n".format(fielddisplay(self, 'coupling', 'Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)'))
+        s = "{}\n".format(fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
+        s = "{}\n".format(fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
+        
+        return s
+    #}}}
 
-        string = "%s\n%s" % (string, fielddisplay(self, "coefficient", "friction coefficient [SI]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "p", "p exponent"))
-        string = "%s\n%s" % (string, fielddisplay(self, "q", "q exponent"))
-        string = "%s\n%s" % (string, fielddisplay(self, 'coupling', 'Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
-        #string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
-        return string
+    def setdefaultparameters(self):  # {{{
+        self.effective_pressure_limit = 0
+
+        return self
     #}}}
 
@@ -50,10 +56,4 @@
     #}}}
 
-    def setdefaultparameters(self):  # {{{
-        #self.requested_outputs = ['default']
-        self.effective_pressure_limit = 0
-        return self
-    #}}}
-
     def defaultoutputs(self, md):  # {{{
         list = []
@@ -64,4 +64,6 @@
         #Early return
         if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
+            return md
+        if solution == 'TransientSoluition' and not md.transient.isstressbalance and not md.transient.isthermal:
             return md
 
@@ -74,6 +76,6 @@
             md = checkfield(md, 'fieldname', 'friction.effective_pressure', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
         elif self.coupling > 4:
-            raise ValueError('md.friction.coupling larger than 4, not supported yet')
-        #md = checkfield(md, 'fieldname', 'friction.requested_outputs', 'stringrow', 1)
+            raise ValueError('not supported yet')
+
         return md
     # }}}
@@ -81,5 +83,13 @@
     def marshall(self, prefix, md, fid):  # {{{
         WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 1, 'format', 'Integer')
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', 1)
+
+        if type(self.coefficient) in [np.ndarray] and (self.coefficient.shape[0] == md.mesh.numberofvertices or self.coefficient.shape[0] == (md.mesh.numberofvertices + 1)):
+            mattype = 1
+            tsl = md.mesh.numberofvertices
+        else:
+            mattype = 2
+            tsl = md.mesh.numberofelements
+
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', mattype, 'timeserieslength', tsl + 1, 'yts', md.constants.yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'p', 'format', 'DoubleMat', 'mattype', 2)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'q', 'format', 'DoubleMat', 'mattype', 2)
@@ -88,15 +98,7 @@
         if self.coupling == 3:
             WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
-        if self.coupling == 4:
-            WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1)
+        elif self.coupling == 4:
+            WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
         elif self.coupling > 4:
-            raise ValueError('md.friction.coupling larger than 4, not supported yet')
-
-    # #process requested outputs
-    #     outputs = self.requested_outputs
-    #     indices = [i for i, x in enumerate(outputs) if x == 'default']
-    #     if len(indices) > 0:
-    #         outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
-    #         outputs = outputscopy
-    #     WriteData(fid, prefix, 'data', outputs, 'name', 'md.friction.requested_outputs', 'format', 'StringArray')
+            raise ValueError('not supported yet')
     # }}}
Index: /issm/trunk-jpl/src/m/classes/frictioncoulomb.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/frictioncoulomb.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/frictioncoulomb.py	(revision 25499)
@@ -1,38 +1,47 @@
+import numpy as np
+
+from checkfield import checkfield
 from fielddisplay import fielddisplay
 from project3d import project3d
-from checkfield import checkfield
 from WriteData import WriteData
 
 
 class frictioncoulomb(object):
-    """
-    FRICTIONCOULOMB class definition
+    """FRICTIONCOULOMB class definition
 
     Usage:
-    frictioncoulomb = frictioncoulomb()
+        frictioncoulomb = frictioncoulomb()
     """
 
     def __init__(self):  # {{{
-        self.coefficient = float('NaN')
-        self.coefficientcoulomb = float('NaN')
-        self.p = float('NaN')
-        self.q = float('NaN')
+        self.coefficient = np.nan
+        self.coefficientcoulomb = np.nan
+        self.p = np.nan
+        self.q = np.nan
         self.coupling = 0
-        self.effective_pressure = float('NaN')
+        self.effective_pressure = np.nan
         self.effective_pressure_limit = 0
-    #set defaults
+
+        #set defaults
         self.setdefaultparameters()
     #}}}
 
     def __repr__(self):  # {{{
-        string = "Basal shear stress parameters: Sigma_b = min(coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b, \n coefficientcoulomb^2 * rho_i * g * (h - h_f)), (effective stress Neff = rho_ice * g * thickness + rho_water * g * bed, r = q / p and s = 1 / p)."
-        string = "%s\n%s" % (string, fielddisplay(self, "coefficient", "power law (Weertman) friction coefficient [SI]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "coefficientcoulomb", "Coulomb friction coefficient [SI]"))
-        string = "%s\n%s" % (string, fielddisplay(self, "p", "p exponent"))
-        string = "%s\n%s" % (string, fielddisplay(self, "q", "q exponent"))
-        string = "%s\n%s" % (string, fielddisplay(self, 'coupling', 'Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure)  and 2 for coupled(not implemented yet)'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
-        return string
+        s = "Basal shear stress parameters: Sigma_b = min(coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b, \n coefficientcoulomb^2 * rho_i * g * (h - h_f)), (effective stress Neff = rho_ice * g * thickness + rho_water * g * bed, r = q / p and s = 1 / p).\n"
+        s = "{}\n".format(fielddisplay(self, "coefficient", "power law (Weertman) friction coefficient [SI]"))
+        s = "{}\n".format(fielddisplay(self, "coefficientcoulomb", "Coulomb friction coefficient [SI]"))
+        s = "{}\n".format(fielddisplay(self, "p", "p exponent"))
+        s = "{}\n".format(fielddisplay(self, "q", "q exponent"))
+        s = "{}\n".format(fielddisplay(self, 'coupling', 'Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure)  and 2 for coupled(not implemented yet)'))
+        s = "{}\n".format(fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
+        s = "{}\n".format(fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
+        
+        return s
+    #}}}
+
+    def setdefaultparameters(self):  # {{{
+        self.effective_pressure_limit = 0
+
+        return self
     #}}}
 
@@ -45,15 +54,10 @@
             self.effective_pressure = project3d(md, 'vector', self.effective_pressure, 'type', 'node', 'layer', 1)
         elif self.coupling == 2:
-            raise ValueError('coupling not supported yet')
+            raise ValueError('not implemented yet')
         elif self.coupling > 2:
-            raise ValueError('md.friction.coupling larger than 2, not supported yet')
+            raise ValueError('not supported yet')
+
         return self
     #}}}
-
-    def setdefaultparameters(self):  # {{{
-        self.effective_pressure_limit = 0
-        return self
-    #}}}
-
     def checkconsistency(self, md, solution, analyses):  # {{{
         #Early return
@@ -65,11 +69,13 @@
         md = checkfield(md, 'fieldname', 'friction.q', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofelements])
         md = checkfield(md, 'fieldname', 'friction.p', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofelements])
+        md = checkfield(md, 'fieldname', 'friction.coupling', 'numel', [1], 'values', [0, 1, 2])
         md = checkfield(md, 'fieldname', 'friction.effective_pressure_limit', 'numel', [1], '>=', 0)
         if self.coupling == 1:
             md = checkfield(md, 'fieldname', 'friction.effective_pressure', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
         elif self.coupling == 2:
-            raise ValueError('coupling not supported yet')
+            raise ValueError('not implemented yet')
         elif self.coupling > 2:
-            raise ValueError('md.friction.coupling larger than 2, not supported yet')
+            raise ValueError('not supported yet')
+
         return md
     # }}}
@@ -86,6 +92,6 @@
             WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
         elif self.coupling == 2:
-            raise ValueError('coupling not supported yet')
+            raise ValueError('not implemented yet')
         elif self.coupling > 2:
-            raise ValueError('md.friction.coupling larger than 2, not supported yet')
+            raise ValueError('not supported yet')
     # }}}
Index: /issm/trunk-jpl/src/m/classes/frictionshakti.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/frictionshakti.m	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/frictionshakti.m	(revision 25499)
@@ -36,6 +36,4 @@
 		end % }}}
 		function marshall(self,prefix,md,fid) % {{{
-			yts=md.constants.yts;
-
 			WriteData(fid,prefix,'name','md.friction.law','data',8,'format','Integer');
 			WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
Index: /issm/trunk-jpl/src/m/classes/frictionshakti.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/frictionshakti.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/frictionshakti.py	(revision 25499)
@@ -1,26 +1,32 @@
+import numpy as np
+
+from checkfield import checkfield
 from fielddisplay import fielddisplay
 from project3d import project3d
-from checkfield import checkfield
 from WriteData import WriteData
 
 
 class frictionshakti(object):
-    """
-    FRICTIONSHAKTI class definition
+    """FRICTIONSHAKTI class definition
 
     Usage:
-    friction = frictionshakti()
+        friction = frictionshakti()
     """
 
     def __init__(self, md):  # {{{
-        self.coefficient = md.friction.coefficient
-    #set defaults
+        self.coefficient = np.nan
+        #set defaults
         self.setdefaultparameters()
     #}}}
 
     def __repr__(self):  # {{{
-        string = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff = rho_ice * g * thickness + rho_water * g * (head - b))"
-        string = "%s\n%s" % (string, fielddisplay(self, "coefficient", "friction coefficient [SI]"))
-        return string
+        s = "Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n(effective stress Neff = rho_ice * g * thickness + rho_water * g * (head - b))\n"
+        s += "{}\n".format(fielddisplay(self, "coefficient", "friction coefficient [SI]"))
+
+        return s
+    #}}}
+
+    def setdefaultparameters(self):  # {{{
+        return self
     #}}}
 
@@ -30,10 +36,6 @@
     #}}}
 
-    def setdefaultparameters(self):  # {{{
-        return self
-    #}}}
-
     def checkconsistency(self, md, solution, analyses):  # {{{
-        #Early return
+        # Early return
         if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
             return md
Index: /issm/trunk-jpl/src/m/classes/frictionwaterlayer.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/frictionwaterlayer.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/frictionwaterlayer.py	(revision 25499)
@@ -1,21 +1,38 @@
+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 frictionwaterlayer(object):
-    """
-    frictionwaterlayer class definition
+    """FRICTIONWATERLAYER class definition
 
     Usage:
         friction = frictionwaterlayer(md)
     """
+
     def __init__(self, md):  # {{{
-        self.coefficient = md.friction.coefficient
-        self.f = float('NaN')
-        self.p = md.friction.p
-        self.q = md.friction.q
-        self.water_layer = float('NaN')
+        self.coefficient = np.nan
+        self.f = np.nan
+        self.p = np.nan
+        self.q = np.nan
+        self.water_layer = np.nan
+    #}}}
+
+    def __repr__(self):  # {{{
+        s = 'Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b * 1 / f(T)\n(effective stress Neff = rho_ice * g * thickness + rho_water * g * (bed + water_layer), r = q / p and s = 1 / p)\n'
+        s = "{}\n".format(fielddisplay(self, 'coefficient', 'frictiontemp coefficient [SI]'))
+        s = "{}\n".format(fielddisplay(self, 'f', 'f variable for effective pressure'))
+        s = "{}\n".format(fielddisplay(self, 'p', 'p exponent'))
+        s = "{}\n".format(fielddisplay(self, 'q', 'q exponent'))
+        s = "{}\n".format(fielddisplay(self, 'water_layer', 'water thickness at the base of the ice (m)'))
+
+        return s
+    #}}}
+
+    def setdefaultparameters(self):  # {{{
+        return self
     #}}}
 
@@ -37,25 +54,14 @@
         self.q = project3d(md, 'vector', self.q, 'type', 'element')
         self.water_layer = project3d(md, 'vector', self.water_layer, 'type', 'node', 'layer', 1)
+
         return self
     # }}}
 
-    def __repr__(self):  # {{{
-        string = 'Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b * 1 / f(T)\n(effective stress Neff = rho_ice * g * thickness + rho_water * g * (bed + water_layer), r = q / p and s = 1 / p)'
-        string = "%s\n%s" % (string, fielddisplay(self, 'coefficient', 'frictiontemp coefficient [SI]'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'f', 'f variable for effective pressure'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'p', 'p exponent'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'q', 'q exponent'))
-        string = "%s\n%s" % (string, fielddisplay(self, 'water_layer', 'water thickness at the base of the ice (m)'))
-
-        return string
-    #}}}
-
     def marshall(self, prefix, md, fid):  #{{{
-        yts = md.constants.yts
         WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 5, 'format', 'Integer')
-        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
+        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
         WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'f', 'format', 'Double')
         WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'p', 'format', 'DoubleMat', 'mattype', 2)
         WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'q', 'format', 'DoubleMat', 'mattype', 2)
-        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'water_layer', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
+        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'water_layer', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
     #}}}
Index: /issm/trunk-jpl/src/m/classes/mesh3dsurface.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh3dsurface.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/mesh3dsurface.py	(revision 25499)
@@ -207,23 +207,23 @@
 
                 #first line:
-                contour                 = OrderedStruct()
-                contour.x               = [self.long[el[0]], self.long[el[1]]]
-                contour.y               = [self.lat[el[0]], self.lat[el[1]]]
-                contour.Geometry        = 'Line' 
-                contours[count]       = contour
+                contour             = OrderedStruct()
+                contour.x           = [self.long[el[0]], self.long[el[1]]]
+                contour.y           = [self.lat[el[0]], self.lat[el[1]]]
+                contour.Geometry    = 'Line' 
+                contours[count]     = contour
 
                 #second line:
-                contour                 = OrderedStruct()
-                contour.x               = [self.long[el[1]], self.long[el[2]]]
-                contour.y               = [self.lat[el[1]], self.lat[el[2]]]
-                contour.Geometry        = 'Line'
-                contours[count + 1]   = contour
+                contour             = OrderedStruct()
+                contour.x           = [self.long[el[1]], self.long[el[2]]]
+                contour.y           = [self.lat[el[1]], self.lat[el[2]]]
+                contour.Geometry    = 'Line'
+                contours[count + 1] = contour
 
                 #second line:
-                contour                 = OrderedStruct()
-                contour.x               = [self.long[el[2]], self.long[el[0]]]
-                contour.y               = [self.lat[el[2]], self.lat[el[0]]]
-                contour.Geometry        = 'Line'
-                contours[count + 2]   = contour
+                contour             = OrderedStruct()
+                contour.x           = [self.long[el[2]], self.long[el[0]]]
+                contour.y           = [self.lat[el[2]], self.lat[el[0]]]
+                contour.Geometry    = 'Line'
+                contours[count + 2] = contour
 
                 #increase count:
Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 25499)
@@ -214,5 +214,5 @@
 			%
 			%   This routine collapses a 3d model into a 2d model
-			%   and collapses all the fileds of the 3d model by
+			%   and collapses all the fields of the 3d model by
 			%   taking their depth-averaged values
 			%
@@ -261,29 +261,66 @@
 
 			%observations
-			if ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end;
-			if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end;
-			if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end;
-			if ~isnan(md.inversion.thickness_obs), md.inversion.thickness_obs=project2d(md,md.inversion.thickness_obs,md.mesh.numberoflayers); end;
-			if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end;
-			if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end;
-			if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end;
+			if ~isnan(md.inversion.vx_obs),
+				md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers);
+			end
+			if ~isnan(md.inversion.vy_obs),
+				md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers);
+			end
+			if ~isnan(md.inversion.vel_obs),
+				md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers);
+			end
+			if ~isnan(md.inversion.thickness_obs),
+				md.inversion.thickness_obs=project2d(md,md.inversion.thickness_obs,md.mesh.numberoflayers);
+			end
+			if ~isnan(md.inversion.cost_functions_coefficients),
+				md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers);
+			end
+			if numel(md.inversion.min_parameters)>1,
+				md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers);
+			end
+			if numel(md.inversion.max_parameters)>1,
+				md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers);
+			end
 			if isa(md.smb,'SMBforcing') & ~isnan(md.smb.mass_balance),
 				md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers); 
 			elseif isa(md.smb,'SMBhenning') & ~isnan(md.smb.smbref),
 				md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers);
-			end;
+			end
 
 			%results
-			if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end;
-			if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end;
-			if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end;
-			if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end;
-			if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end;
-			if ~isnan(md.initialization.pressure),md.initialization.pressure=project2d(md,md.initialization.pressure,1);end;
-			if ~isnan(md.initialization.sediment_head),md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);end;
-			if ~isnan(md.initialization.epl_head),md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);end;
-			if ~isnan(md.initialization.epl_thickness),md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);end;
-			if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project2d(md,md.initialization.waterfraction,1);end;
-			if ~isnan(md.initialization.watercolumn),md.initialization.watercolumn=project2d(md,md.initialization.watercolumn,1);end;
+			if ~isnan(md.initialization.vx),
+				md.initialization.vx=DepthAverage(md,md.initialization.vx);
+			end
+			if ~isnan(md.initialization.vy),
+				md.initialization.vy=DepthAverage(md,md.initialization.vy);
+			end
+			if ~isnan(md.initialization.vz),
+				md.initialization.vz=DepthAverage(md,md.initialization.vz);
+			end
+			if ~isnan(md.initialization.vel),
+				md.initialization.vel=DepthAverage(md,md.initialization.vel);
+			end
+			if ~isnan(md.initialization.temperature),
+				md.initialization.temperature=DepthAverage(md,md.initialization.temperature);
+			end
+			if ~isnan(md.initialization.pressure),
+				md.initialization.pressure=project2d(md,md.initialization.pressure,1);
+			end
+			if ~isnan(md.initialization.sediment_head),
+				md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);
+			end
+			if ~isnan(md.initialization.epl_head),
+				md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);
+			end
+			if ~isnan(md.initialization.epl_thickness),
+				md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);
+			end
+			if ~isnan(md.initialization.waterfraction),
+				md.initialization.waterfraction=project2d(md,md.initialization.waterfraction,1);
+			end
+			if ~isnan(md.initialization.watercolumn),
+				md.initialization.watercolumn=project2d(md,md.initialization.watercolumn,1);
+			end
+
 			%giaivins
 			if isa(md.gia,'giaivins'),
@@ -299,5 +336,5 @@
 				md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1);
 				md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1);
-			end	
+			end
 
 			%boundary conditions
@@ -307,7 +344,13 @@
 			md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers);
 			md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers);
-			if numel(md.masstransport.spcthickness)>1 md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers); end
-			if numel(md.damage.spcdamage)>1, md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers); end
-			if numel(md.levelset.spclevelset)>1, md.levelset.spclevelset=project2d(md,md.levelset.spclevelset,md.mesh.numberoflayers); end
+			if numel(md.masstransport.spcthickness)>1,
+				md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers);
+			end
+			if numel(md.damage.spcdamage)>1,
+				md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers);
+			end
+			if numel(md.levelset.spclevelset)>1,
+				md.levelset.spclevelset=project2d(md,md.levelset.spclevelset,md.mesh.numberoflayers);
+			end
 			md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
 
@@ -362,5 +405,4 @@
 				md.geometry.bed=project2d(md,md.geometry.bed,1);
 			end
-
 			if ~isnan(md.mask.ocean_levelset),
 				md.mask.ocean_levelset=project2d(md,md.mask.ocean_levelset,1);
@@ -371,6 +413,10 @@
 
 			%lat long
-			if numel(md.mesh.lat) ==md.mesh.numberofvertices,  md.mesh.lat=project2d(md,md.mesh.lat,1); end
-			if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end
+			if numel(md.mesh.lat)==md.mesh.numberofvertices,
+				md.mesh.lat=project2d(md,md.mesh.lat,1);
+			end
+			if numel(md.mesh.long)==md.mesh.numberofvertices,
+				md.mesh.long=project2d(md,md.mesh.long,1);
+			end
 
 			%outputdefinitions
@@ -388,5 +434,5 @@
 			end
 
-			%Initialize with the 2d mesh
+			%Initialize the 2d mesh
 			mesh=mesh2d();
 			mesh.x=md.mesh.x2d;
@@ -395,10 +441,20 @@
 			mesh.numberofelements=md.mesh.numberofelements2d;
 			mesh.elements=md.mesh.elements2d;
-			if numel(md.mesh.lat) ==md.mesh.numberofvertices,  mesh.lat=project2d(md,md.mesh.lat,1); end
-			if numel(md.mesh.long)==md.mesh.numberofvertices,	mesh.long=project2d(md,md.mesh.long,1); end
+			if numel(md.mesh.lat)==md.mesh.numberofvertices,
+				mesh.lat=project2d(md,md.mesh.lat,1);
+			end
+			if numel(md.mesh.long)==md.mesh.numberofvertices,
+				mesh.long=project2d(md,md.mesh.long,1);
+			end
 			mesh.epsg=md.mesh.epsg;
-			if numel(md.mesh.scale_factor)==md.mesh.numberofvertices,	mesh.scale_factor=project2d(md,md.mesh.scale_factor,1); end
-			if ~isnan(md.mesh.vertexonboundary), mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); end
-			if ~isnan(md.mesh.elementconnectivity), mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); end
+			if numel(md.mesh.scale_factor)==md.mesh.numberofvertices,
+				mesh.scale_factor=project2d(md,md.mesh.scale_factor,1);
+			end
+			if ~isnan(md.mesh.vertexonboundary),
+				mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1);
+			end
+			if ~isnan(md.mesh.elementconnectivity),
+				mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1);
+			end
 			md.mesh=mesh;
 			md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
@@ -622,5 +678,6 @@
 				md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
 				md2.mesh.segments=contourenvelope(md2.mesh);
-				md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
+				md2.mesh.vertexonboundary=zeros(numberofvertices2,1);
+				md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
 			else
 				%First do the connectivity for the contourenvelope in 2d
@@ -628,5 +685,6 @@
 				md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity);
 				segments=contourenvelope(md2.mesh);
-				md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(segments(:,1:2))=1;
+				md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1);
+				md2.mesh.vertexonboundary(segments(:,1:2))=1;
 				md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1);
 				%Then do it for 3d as usual
Index: /issm/trunk-jpl/src/m/classes/model.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/model.py	(revision 25499)
@@ -412,5 +412,5 @@
             md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements, md2.mesh.numberofvertices)
             md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements, md2.mesh.vertexconnectivity)
-            md2.mesh.segments = contourenvelope(md2)
+            md2.mesh.segments = contourenvelope(md2.mesh)
             md2.mesh.vertexonboundary = np.zeros(numberofvertices2, bool)
             md2.mesh.vertexonboundary[md2.mesh.segments[:, 0:2] - 1] = True
@@ -419,5 +419,5 @@
             md2.mesh.vertexconnectivity = NodeConnectivity(md2.mesh.elements2d, md2.mesh.numberofvertices2d)
             md2.mesh.elementconnectivity = ElementConnectivity(md2.mesh.elements2d, md2.mesh.vertexconnectivity)
-            segments = contourenvelope(md2)
+            msegments = contourenvelope(md2.mesh)
             md2.mesh.vertexonboundary = np.zeros(int(numberofvertices2 / md2.mesh.numberoflayers), bool)
             md2.mesh.vertexonboundary[segments[:, 0:2] - 1] = True
@@ -696,39 +696,49 @@
 
         This routine collapses a 3d model into a 2d model and collapses all
-        the fileds of the 3d model by taking their depth - averaged values
+        the fields of the 3d model by taking their depth-averaged values
 
         Usage:
             md = collapse(md)
+
+        See also: EXTRUDE, MODELEXTRACT
         """
 
-        #Check that the model is really a 3d model
-        if md.mesh.domaintype().lower() != '3d':
-            raise Exception("only a 3D model can be collapsed")
-
-        #dealing with the friction law
-        #drag is limited to nodes that are on the bedrock.
-        if hasattr(md.friction, 'coefficient'):
+        # Check that the model is really a 3d model
+        if md.mesh.elementtype() != 'Penta':
+            raise Exception("collapse error message: only a 3d mesh can be collapsed")
+
+        # Start with changing all the fields from the 3d mesh
+
+        # Dealing with the friction law
+        # Drag is limited to nodes that are on the bedrock.
+        if md.friction.__class__.__name__ == 'friction':
             md.friction.coefficient = project2d(md, md.friction.coefficient, 1)
-
-        #p and q (same deal, except for element that are on the bedrock: )
-        if hasattr(md.friction, 'p'):
             md.friction.p = project2d(md, md.friction.p, 1)
-        if hasattr(md.friction, 'q'):
             md.friction.q = project2d(md, md.friction.q, 1)
-
-        if hasattr(md.friction, 'coefficientcoulomb'):
+        elif md.friction.__class__.__name__ == 'frictioncoulomb':
+            md.friction.coefficient = project2d(md, md.friction.coefficient, 1)
             md.friction.coefficientcoulomb = project2d(md, md.friction.coefficientcoulomb, 1)
-        if hasattr(md.friction, 'C'):
+            md.friction.p = project2d(md, md.friction.p, 1)
+            md.friction.q = project2d(md, md.friction.q, 1)
+        elif md.friction.__class__.__name__ == 'frictionhydro':
+            md.friction.q = project2d(md, md.friction.q, 1)
             md.friction.C = project2d(md, md.friction.C, 1)
-        if hasattr(md.friction, 'As'):
             md.friction.As = project2d(md, md.friction.As, 1)
-        if hasattr(md.friction, 'effective_pressure') and not np.isnan(md.friction.effective_pressure).all():
             md.friction.effective_pressure = project2d(md, md.friction.effective_pressure, 1)
-        if hasattr(md.friction, 'water_layer'):
+        elif md.friction.__class__.__name__ == 'frictionwaterlayer':
+            md.friction.coefficient = project2d(md, md.friction.coefficient, 1)
+            md.friction.p = project2d(md, md.friction.p, 1)
+            md.friction.q = project2d(md, md.friction.q, 1)
             md.friction.water_layer = project2d(md, md.friction.water_layer, 1)
-        if hasattr(md.friction, 'm'):
+        elif md.friction.__class__.__name__ == 'frictionweertman':
+            md.friction.C = project2d(md, md.friction.C, 1)
             md.friction.m = project2d(md, md.friction.m, 1)
-
-        #observations
+        elif md.friction.__class__.__name__ == 'frictionweertmantemp':
+            md.friction.C = project2d(md, md.friction.C, 1)
+            md.friction.m = project2d(md, md.friction.m, 1)
+        else:
+            print('friction type not supported')
+
+        # Observations
         if not np.isnan(md.inversion.vx_obs).all():
             md.inversion.vx_obs = project2d(md, md.inversion.vx_obs, md.mesh.numberoflayers)
@@ -741,15 +751,14 @@
         if not np.isnan(md.inversion.cost_functions_coefficients).all():
             md.inversion.cost_functions_coefficients = project2d(md, md.inversion.cost_functions_coefficients, md.mesh.numberoflayers)
-        if isinstance(md.inversion.min_parameters, np.ndarray):
-            if md.inversion.min_parameters.size > 1:
-                md.inversion.min_parameters = project2d(md, md.inversion.min_parameters, md.mesh.numberoflayers)
-        if isinstance(md.inversion.max_parameters, np.ndarray):
-            if md.inversion.max_parameters.size > 1:
-                md.inversion.max_parameters = project2d(md, md.inversion.max_parameters, md.mesh.numberoflayers)
-        if hasattr(md.smb, 'mass_balance'):
-            if not np.isnan(md.smb.mass_balance).all():
+        if isinstance(md.inversion.min_parameters, np.ndarray) and md.inversion.min_parameters.size > 1:
+            md.inversion.min_parameters = project2d(md, md.inversion.min_parameters, md.mesh.numberoflayers)
+        if isinstance(md.inversion.max_parameters, np.ndarray) and md.inversion.max_parameters.size > 1:
+            md.inversion.max_parameters = project2d(md, md.inversion.max_parameters, md.mesh.numberoflayers)
+        if md.smb.__class__.__name__ == 'SMBforcing' and not np.isnan(md.smb.mass_balance).all():
                 md.smb.mass_balance = project2d(md, md.smb.mass_balance, md.mesh.numberoflayers)
-
-        #results
+        elif md.smb.__class__.__name__ == 'SMBhenning' and not np.isnan(md.smb.smbref).all():
+                md.smb.smbref = project2d(md, md.smb.smbref, md.mesh.numberoflayers)
+
+        # Results
         if not np.isnan(md.initialization.vx).all():
             md.initialization.vx = DepthAverage(md, md.initialization.vx)
@@ -770,6 +779,10 @@
         if not np.isnan(md.initialization.epl_thickness).all():
             md.initialization.epl_thickness = project2d(md, md.initialization.epl_thickness, 1)
-
-        #giaivins
+        if not np.isnan(md.initialization.waterfraction).all():
+            md.initialization.waterfraction = project2d(md, md.initialization.waterfraction, 1)
+        if not np.isnan(md.initialization.watercolumn).all():
+            md.initialization.watercolumn = project2d(md, md.initialization.watercolumn, 1)
+
+        # giaivins
         if md.gia.__class__.__name__ == 'giaivins':
             if not np.isnan(md.gia.mantle_viscosity).all():
@@ -778,5 +791,5 @@
                 md.gia.lithosphere_thickness = project2d(md, md.gia.lithosphere_thickness, 1)
 
-        #elementstype
+        # elementstype
         if not np.isnan(md.flowequation.element_equation).all():
             md.flowequation.element_equation = project2d(md, md.flowequation.element_equation, 1)
@@ -786,17 +799,5 @@
             md.flowequation.borderFS = project2d(md, md.flowequation.borderFS, 1)
 
-        # Hydrologydc variables
-        hydrofields = md.hydrology.__dict__.keys()
-        for field in hydrofields:
-            try:
-                isvector = np.logical_or(np.shape(md.hydrology.__dict__[field])[0] == md.mesh.numberofelements,
-                                         np.shape(md.hydrology.__dict__[field])[0] == md.mesh.numberofvertices)
-            except IndexError:
-                isvector = False
-            #we colpase only fields that are vertices or element based
-            if isvector:
-                md.hydrology.__dict__[field] = project2d(md, md.hydrology.__dict__[field], 1)
-
-        #boundary conditions
+        # Boundary conditions
         md.stressbalance.spcvx = project2d(md, md.stressbalance.spcvx, md.mesh.numberoflayers)
         md.stressbalance.spcvy = project2d(md, md.stressbalance.spcvy, md.mesh.numberoflayers)
@@ -804,4 +805,10 @@
         md.stressbalance.referential = project2d(md, md.stressbalance.referential, md.mesh.numberoflayers)
         md.stressbalance.loadingforce = project2d(md, md.stressbalance.loadingforce, md.mesh.numberoflayers)
+        # TODO: 
+        # - Check if md.mesh.numberoflayershould really be offset by 1.
+        # - Find out why md.masstransport.spcthickness is not offset, but the 
+        #   other fields are.
+        # - If offset is required, figure out if it can be abstarcted away to 
+        #   another part of the API.
         if np.size(md.masstransport.spcthickness) > 1:
             md.masstransport.spcthickness = project2d(md, md.masstransport.spcthickness, md.mesh.numberoflayers)
@@ -812,38 +819,76 @@
         md.thermal.spctemperature = project2d(md, md.thermal.spctemperature, md.mesh.numberoflayers - 1)
 
-        #materials
+        # Hydrologydc variables
+        if md.hydrology.__class__.__name__ == 'hydrologydc':
+            md.hydrology.spcsediment_head = project2d(md, md.hydrology.spcsediment_head, 1)
+            md.hydrology.mask_eplactive_node = project2d(md, md.hydrology.mask_eplactive_node, 1)
+            md.hydrology.sediment_transmitivity = project2d(md, md.hydrology.sediment_transmitivity, 1)
+            md.hydrology.basal_moulin_input = project2d(md, md.hydrology.basal_moulin_input, 1)
+            if md.hydrology.isefficientlayer == 1:
+                md.hydrology.spcepl_head = project2d(md, md.hydrology.spcepl_head, 1)
+            # hydrofields = md.hydrology.__dict__.keys()
+            # for field in hydrofields:
+            #     try:
+            #         isvector = np.logical_or(np.shape(md.hydrology.__dict__[field])[0] == md.mesh.numberofelements,
+            #                                  np.shape(md.hydrology.__dict__[field])[0] == md.mesh.numberofvertices)
+            #     except IndexError:
+            #         isvector = False
+            #     #we collapse only fields that are vertices or element based
+            #     if isvector:
+            #         md.hydrology.__dict__[field] = project2d(md, md.hydrology.__dict__[field], 1)
+
+        # Materials
         md.materials.rheology_B = DepthAverage(md, md.materials.rheology_B)
         md.materials.rheology_n = project2d(md, md.materials.rheology_n, 1)
 
-        #damage:
+        # Damage
         if md.damage.isdamage:
             md.damage.D = DepthAverage(md, md.damage.D)
 
-        #special for thermal modeling:
-        md.basalforcings.groundedice_melting_rate = project2d(md, md.basalforcings.groundedice_melting_rate, 1)
-        md.basalforcings.floatingice_melting_rate = project2d(md, md.basalforcings.floatingice_melting_rate, 1)
+        # Special for thermal modeling
+        if not np.isnan(md.basalforcings.groundedice_melting_rate).all():
+            md.basalforcings.groundedice_melting_rate = project2d(md, md.basalforcings.groundedice_melting_rate, 1)
+        if hasattr(md.basalforcings, 'floatingice_melting_rate') and not np.isnan(md.basalforcings.floatingice_melting_rate).all():
+            md.basalforcings.floatingice_melting_rate = project2d(md, md.basalforcings.floatingice_melting_rate, 1)
         md.basalforcings.geothermalflux = project2d(md, md.basalforcings.geothermalflux, 1)  #bedrock only gets geothermal flux
 
-        #update of connectivity matrix
+        if hasattr(md.calving, 'coeff') and not np.isnan(md.calving.coeff).all():
+            md.calving.coeff = project2d(md, md.calving.coeff, 1)
+        if hasattr(md.frontalforcings, 'meltingrate') and not np.isnan(md.frontalforcings.meltingrate).all():
+            md.frontalforcings.meltingrate = project2d(md, md.frontalforcings.meltingrate, 1)
+
+        # Update of connectivity matrix
         md.mesh.average_vertex_connectivity = 25
 
-        #parameters
+        # Collapse the mesh
+        nodes2d = md.mesh.numberofvertices2d
+        elements2d = md.mesh.numberofelements2d
+
+        # Parameters
         md.geometry.surface = project2d(md, md.geometry.surface, 1)
         md.geometry.thickness = project2d(md, md.geometry.thickness, 1)
         md.geometry.base = project2d(md, md.geometry.base, 1)
-        if isinstance(md.geometry.bed, np.ndarray):
+        if not np.isnan(md.geometry.bed).all():
             md.geometry.bed = project2d(md, md.geometry.bed, 1)
-        md.mask.ocean_levelset = project2d(md, md.mask.ocean_levelset, 1)
-        md.mask.ice_levelset = project2d(md, md.mask.ice_levelset, 1)
-
-        #OutputDefinitions
+        if not np.isnan(md.mask.ocean_levelset).all():
+            md.mask.ocean_levelset = project2d(md, md.mask.ocean_levelset, 1)
+        if not np.isnan(md.mask.ice_levelset).all():
+            md.mask.ice_levelset = project2d(md, md.mask.ice_levelset, 1)
+
+        # Lat/long
+        if np.size(md.mesh.lat) == md.mesh.numberofvertices:
+            md.mesh.lat = project2d(md, md.mesh.lat, 1)
+        if np.size(md.mesh.long) == md.mesh.numberofvertices:
+            md.mesh.long = project2d(md, md.mesh.long, 1)
+
+        # OutputDefinitions
         if md.outputdefinition.definitions:
             for solutionfield, field in list(md.outputdefinition.__dict__.items()):
                 if isinstance(field, list):
-                    #get each definition
+                    # Get each definition
                     for i, fieldi in enumerate(field):
                         if fieldi:
                             fieldr = getattr(md.outputdefinition, solutionfield)[i]
-                            #get subfields
+                            # Get subfields
                             for solutionsubfield, subfield in list(fieldi.__dict__.items()):
                                 if np.size(subfield) == md.mesh.numberofvertices:
@@ -852,5 +897,5 @@
                                     setattr(fieldr, solutionsubfield, project2d(md, subfield, 1))
 
-        #Initialize with the 2d mesh
+        # Initialize he 2d mesh
         mesh = mesh2d()
         mesh.x = md.mesh.x2d
@@ -859,22 +904,23 @@
         mesh.numberofelements = md.mesh.numberofelements2d
         mesh.elements = md.mesh.elements2d
+        # if not np.isnan(md.mesh.vertexonboundary).all():
+        #     mesh.vertexonboundary = project2d(md, md.mesh.vertexonboundary, 1)
+        # if not np.isnan(md.mesh.elementconnectivity).all():
+        #     mesh.elementconnectivity = project2d(md, md.mesh.elementconnectivity, 1)
+        if np.size(md.mesh.lat) == md.mesh.numberofvertices:
+            mesh.lat = project2d(md, md.mesh.lat, 1)
+        if np.size(md.mesh.long) == md.mesh.numberofvertices:
+            mesh.long = project2d(md, md.mesh.long, 1)
+        mesh.epsg = md.mesh.epsg
+        if np.size(md.mesh.scale_factor) == md.mesh.numberofvertices:
+            mesh.scale_factor = project2d(md, md.mesh.scale_factor, 1)
         if not np.isnan(md.mesh.vertexonboundary).all():
-            mesh.vertexonboundary = project2d(md, md.mesh.vertexonboundary, 1)
-        if not np.isnan(md.mesh.elementconnectivity).all():
-            mesh.elementconnectivity = project2d(md, md.mesh.elementconnectivity, 1)
-        if isinstance(md.mesh.lat, np.ndarray):
-            if md.mesh.lat.size == md.mesh.numberofvertices:
-                mesh.lat = project2d(md, md.mesh.lat, 1)
-        if isinstance(md.mesh.long, np.ndarray):
-            if md.mesh.long.size == md.mesh.numberofvertices:
-                md.mesh.long = project2d(md, md.mesh.long, 1)
-        mesh.epsg = md.mesh.epsg
-        if isinstance(md.mesh.scale_factor, np.ndarray):
-            if md.mesh.scale_factor.size == md.mesh.numberofvertices:
-                md.mesh.scale_factor = project2d(md, md.mesh.scale_factor, 1)
+            mesh.vertexonboundary= project2d(md, md.mesh.vertexonboundary, 1)
+        if not np.isnan(md.mesh.elementonboundary).all():
+            mesh.elementonboundary = project2d(md, md.mesh.elementonboundary, 1)
         md.mesh = mesh
         md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
         md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
-        md.mesh.segments = contourenvelope(md)
+        md.mesh.segments = contourenvelope(md.mesh)
 
         return md
Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/sealevelmodel.m	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/sealevelmodel.m	(revision 25499)
@@ -174,5 +174,6 @@
 			
 			%initialize, to avoid issues of having more transitions than meshes. 
-			self.transitions={}; self.eltransitions={};
+			self.transitions={};
+			self.eltransitions={};
 
 			%for elements:
@@ -180,5 +181,5 @@
 			ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
 			ze=self.earth.mesh.z(self.earth.mesh.elements)*[1;1;1]/3;
-
+			
 			for i=1:length(self.icecaps),
 				mdi=self.icecaps{i};
@@ -196,5 +197,4 @@
 				self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force);
 			end
-
 		end % }}}
 		function checkintersections(self) % {{{
@@ -397,5 +397,5 @@
 						eval(['capfieldj= self.icecaps{' num2str(j) '}.' string ';']);
 						if ~isequal(capfieldi(end,:),capfieldj(end,:)),
-							error(['Time stamps for ' string ' field is different between icecaps ' num2str(i) ' and ' num2str(j)]);
+							error(['Time stamps for ' string ' field are different between icecaps ' num2str(i) ' and ' num2str(j)]);
 						end
 					end
Index: /issm/trunk-jpl/src/m/classes/sealevelmodel.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/sealevelmodel.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/sealevelmodel.py	(revision 25499)
@@ -1,7 +1,8 @@
+from copy import deepcopy
 import warnings
 import numpy as np
 from fielddisplay import fielddisplay
-from fieldnames import fieldnames
 from generic import generic
+from getsubattr import getsubattr
 from issmsettings import issmsettings
 from meshintersect3d import meshintersect3d
@@ -10,10 +11,10 @@
 from modelmerge3d import modelmerge3d
 from pairoptions import pairoptions
+from planetradius import planetradius
+from plotmodel import plotmodel
 from private import private
+from setsubattr import setsubattr
 from TwoDToThreeD import TwoDToThreeD
-from plotmodel import plotmodel
-from plot3 import plot3
-from loneedges import loneedges
-from planetradius import planetradius
+
 
 class sealevelmodel(object):
@@ -174,8 +175,7 @@
 
         # For elements
-        onesmatrix = np.array([[1], [1], [1]])
-        xe = self.earth.mesh.x[self.earth.mesh.elements] * onesmatrix / 3
-        ye = self.earth.mesh.y[self.earth.mesh.elements] * onesmatrix / 3
-        ze = self.earth.mesh.z[self.earth.mesh.elements] * onesmatrix / 3
+        xe = np.mean(self.earth.mesh.x[self.earth.mesh.elements - 1], axis=1)
+        ye = np.mean(self.earth.mesh.y[self.earth.mesh.elements - 1], axis=1)
+        ze = np.mean(self.earth.mesh.z[self.earth.mesh.elements - 1], axis=1)
 
         for i in range(len(self.icecaps)):
@@ -184,7 +184,7 @@
 
             # For elements
-            xei = mdi.mesh.x[mdi.mesh.elements] * onesmatrix / 3
-            yei = mdi.mesh.y[mdi.mesh.elements] * onesmatrix / 3
-            zei = mdi.mesh.z[mdi.mesh.elements] * onesmatrix / 3
+            xei = np.mean(mdi.mesh.x[mdi.mesh.elements - 1], axis=1)
+            yei = np.mean(mdi.mesh.y[mdi.mesh.elements - 1], axis=1)
+            zei = np.mean(mdi.mesh.z[mdi.mesh.elements - 1], axis=1)
 
             print('Computing vertex intersections for basin {}'.format(self.basins[i].name))
@@ -300,5 +300,5 @@
 
         # Make 3D model
-        models = self.icecaps
+        models = deepcopy(self.icecaps)
         for i in range(len(models)):
             models[i] = TwoDToThreeD(models[i], self.planet)
@@ -308,5 +308,5 @@
         for i in range(1, len(models)):
             md = modelmerge3d(md, models[i], 'tolerance', tolerance)
-            md.private.bamg.landmask = np.vstack((md.private.bamg.landmask, models[i].private.bamg.landmask))
+            md.private.bamg.landmask = np.hstack((md.private.bamg.landmask, models[i].private.bamg.landmask))
 
         # Look for lone edges if asked for it
@@ -317,4 +317,5 @@
                 ind1 = edges(i, 1)
                 ind2 = edges(i, 2)
+                # TODO: Reconfigure the following in the process of bringing plotting online
                 plot3([md.mesh.x[ind1], md.mesh.x[ind2]], [md.mesh.y[ind1], md.mesh.y[ind2]], [md.mesh.z[ind1], md.mesh.z[ind2]], 'g*-')
 
@@ -323,5 +324,5 @@
 
         # Create mesh radius
-        self.earth.mesh.r = planetradius('earth')
+        self.earth.mesh.r = planetradius('earth') * np.ones((md.mesh.numberofvertices, ))
     #}}}
 
@@ -349,20 +350,21 @@
     def transfer(self, string):  #{{{
         # Recover field size in one icecap
-        n = np.size(getattr(self.icecaps[0], string), 0)
+        n = getsubattr(self.icecaps[0], string).shape[0]
 
         if n == self.icecaps[0].mesh.numberofvertices:
-            setattr(self.earth, string, np.zeros((self.earth.mesh.numberofvertices, )))
+            setsubattr(self.earth, string, np.zeros((self.earth.mesh.numberofvertices, ))) # Assign array of zeros to target attribute
+            earth_attr = getsubattr(self.earth, string) # Retrieve reference to target attribute
             for i in range(len(self.icecaps)):
-                getattr(self.earth, string)[self.transitions[i]] = getattr(self.icecaps[i], string)
-        elif n == (self.self.icecaps[0].mesh.numberofvertices + 1):
+                earth_attr[self.transitions[i]] = getsubattr(self.icecaps[i], string)
+        elif n == (self.icecaps[0].mesh.numberofvertices + 1):
             # Dealing with transient dataset
             # Check that all timetags are similar between all icecaps #{{{
             for i in range(len(self.icecaps)):
-                capfieldi = getattr(self.icecaps[i], string)
-                for j in range(1, len(self.icecaps)):
-                    capfieldj = getattr(self.icecaps[j], string)
+                capfieldi = getsubattr(self.icecaps[i], string)
+                for j in range(i + 1, len(self.icecaps)):
+                    capfieldj = getsubattr(self.icecaps[j], string)
                     if capfieldi[-1, :] != capfieldj[-1, :]:
-                        raise Exception("Time stamps for {} field is different between icecaps {} and {}".format(string, i, j))
-            capfield1 = getattr(self.icecaps[0], string)
+                        raise Exception("Time stamps for {} field are different between icecaps {} and {}".format(string, i, j))
+            capfield1 = getsubattr(self.icecaps[0], string)
             times = capfield1[-1, :]
             nsteps = len(times)
@@ -374,13 +376,14 @@
             # Transfer all the time fields #{{{
             for i in range(len(self.icecaps)):
-                capfieldi = getattr(self.icecaps[i], string)
+                capfieldi = getsubattr(self.icecaps[i], string)
                 for j in range(nsteps):
-                    field[self.transitions[i], j] = capfieldi[0:-2, j]  # Transfer only the values, not the time
-            setattr(self.earth, string, field)  # Do not forget to plug the field variable into its final location
+                    field[self.transitions[i], j] = capfieldi[0:-1, j]  # Transfer only the values, not the time
+            setsubattr(self.earth, string, field)  # Do not forget to plug the field variable into its final location
             #}}}
         elif n == (self.icecaps[0].mesh.numberofelements):
-            setattr(self.earth, string, np.zeros((self.earth.mesh.numberofvertices, )))
+            setsubattr(self.earth, string, np.zeros((self.earth.mesh.numberofelements, ))) # Assign array of zeros to target attribute
+            earth_attr = getsubattr(self.earth, string) # Retrieve reference to target attribute
             for i in range(len(self.icecaps)):
-                getattr(self.earth, string)[self.eltransitions[i]] = getattr(self.icecaps[i], string)
+                earth_attr[self.eltransitions[i]] = getsubattr(self.icecaps[i], string)
         else:
             raise Exception('not supported yet')
Index: /issm/trunk-jpl/src/m/classes/surfaceload.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/surfaceload.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/classes/surfaceload.py	(revision 25499)
@@ -7,10 +7,9 @@
 
 class surfaceload(object):
-    '''
-    SURFACELOAD class definition
+    """SURFACELOAD class definition
 
-        Usage:
-            surfaceload = surfaceload()
-    '''
+    Usage:
+        surfaceload = surfaceload()
+    """
 
     def __init__(self, *args): #{{{
@@ -58,11 +57,11 @@
     def marshall(self, prefix, md, fid): #{{{
         if len(self.icethicknesschange) == 0:
-            self.icethicknesschange = np.zeros(md.mesh.numberofelements + 1)
+            self.icethicknesschange = np.zeros((md.mesh.numberofelements + 1, ))
 
         if len(self.waterheightchange) == 0:
-            self.waterheightchange = np.zeros(md.mesh.numberofelements + 1)
+            self.waterheightchange = np.zeros((md.mesh.numberofelements + 1, ))
 
         if len(self.other) == 0:
-            self.other = np.zeros(md.mesh.numberofelements + 1)
+            self.other = np.zeros((md.mesh.numberofelements + 1, ))
 
         WriteData(fid, prefix, 'object', self, 'fieldname', 'icethicknesschange', 'name', 'md.solidearth.surfaceload.icethicknesschange', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', md.constants.yts)
Index: /issm/trunk-jpl/src/m/consistency/checkfield.py
===================================================================
--- /issm/trunk-jpl/src/m/consistency/checkfield.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/consistency/checkfield.py	(revision 25499)
@@ -8,15 +8,15 @@
 
 def checkfield(md, *args):
-    '''
-    CHECKFIELD - check field consistency
-
-        Used to check model consistency
+    """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:
+    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' 
@@ -34,7 +34,7 @@
         - message: overloaded error message
 
-        Usage:
-            md = checkfield(md, fieldname, options)
-    '''
+    Usage:
+        md = checkfield(md, fieldname, options)
+    """
 
     #get options
@@ -76,12 +76,12 @@
                 #Check that vector size will not be confusing for ModelProcessorx
                 if (md.mesh.numberofvertices == md.mesh.numberofelements):
-                    raise RuntimeError('number of vertices is the same as number of elements')
+                    raise Exception('number of vertices is the same as number of elements')
                 elif (md.mesh.numberofvertices + 1 == md.mesh.numberofelements):
-                    raise RuntimeError('number of vertices + 1 is the same as number of elements')
+                    raise Exception('number of vertices + 1 is the same as number of elements')
                 elif (md.mesh.numberofvertices == md.mesh.numberofelements + 1):
-                    raise RuntimeError('number of vertices is the same as number of elements + 1')
+                    raise Exception('number of vertices is the same as number of elements + 1')
 
                 #Uniform field
-                if (np.size(field, 0) == 1):
+                if (field.shape[0] == 1):
                     if (np.ndim(field) > 1 and np.shape(field)[1] != 1):
                         md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))
@@ -112,5 +112,5 @@
 
             else:
-                raise RuntimeError("fieldsize '{}' not supported yet".format(fieldsize))
+                raise Exception("fieldsize '{}' not supported yet".format(fieldsize))
 
         else:
@@ -155,5 +155,5 @@
                 md = md.checkmessage(options.getfieldvalue('message', "field '{}' value should be '{}'".format(fieldname, fieldvalues[0])))
             elif len(fieldvalues) == 2:
-                md = md.checkmessage(options.getfieldvalue('message', "field '{}' values should be '{}' '{}' or '%s'".format(fieldname, fieldvalues[0], fieldvalues[1])))
+                md = md.checkmessage(options.getfieldvalue('message', "field '{}' values should be '{}' '{}' or '{}'".format(fieldname, fieldvalues[0], fieldvalues[1])))
             else:
                 md = md.checkmessage(options.getfieldvalue('message', "field '{}' should have values in {}".format(fieldname, fieldvalues)))
@@ -166,5 +166,5 @@
         if np.size(lowerbound) > 1:  #checking elementwise
             if any(field < lowerbound):
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values above %d" % (fieldname, lowerbound)))
+                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
         else:
             minval = np.nanmin(field)
@@ -178,5 +178,5 @@
 
             if minval < lowerbound:
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values above %d" % (fieldname, lowerbound)))
+                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
 
     if options.exist('>'):
@@ -186,5 +186,5 @@
         if np.size(lowerbound) > 1:  #checking elementwise
             if any(field <= lowerbound):
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values above %d" % (fieldname, lowerbound)))
+                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
         else:
             minval = np.nanmin(field)
@@ -198,5 +198,5 @@
 
             if minval <= lowerbound:
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values above %d" % (fieldname, lowerbound)))
+                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
 
     #check smaller
@@ -207,5 +207,5 @@
         if np.size(upperbound) > 1:  #checking elementwise
             if any(field > upperbound):
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values below %d" % (fieldname, upperbound)))
+                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
         else:
             maxval = np.nanmax(field)
@@ -220,5 +220,5 @@
                 maxval = field.fov_forward_indices[0]
             if maxval > upperbound:
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values below %d" % (fieldname, upperbound)))
+                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
 
     if options.exist('<'):
@@ -228,5 +228,5 @@
         if np.size(upperbound) > 1:  #checking elementwise
             if any(field >= upperbound):
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values below %d" % (fieldname, upperbound)))
+                md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
 
         else:
@@ -241,41 +241,41 @@
 
                 if maxval >= upperbound:
-                    md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have values below %d" % (fieldname, upperbound)))
+                    md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
 
     #check file
     if options.getfieldvalue('file', 0):
         if not os.path.exists(field):
-            md = md.checkmessage("file provided in '%s': '%s' does not exist" % (fieldname, field))
+            md = md.checkmessage("file provided in {}: {} does not exist".format(fieldname, field))
 
     #Check row of strings
     if options.exist('stringrow'):
         if not isinstance(field, list):
-            md = md.checkmessage(options.getfieldvalue('message', "field '%s' should be a list" % fieldname))
+            md = md.checkmessage(options.getfieldvalue('message', "field {} should be a list".format(fieldname)))
 
     #Check forcings (size and times)
     if options.getfieldvalue('timeseries', 0):
-        if np.size(field, 0) == md.mesh.numberofvertices or np.size(field, 0) == md.mesh.numberofelements:
+        if field.shape[0] == md.mesh.numberofvertices or field.shape[0] == md.mesh.numberofelements:
             if np.ndim(field) > 1 and not np.size(field, 1) == 1:
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname))
-        elif np.size(field, 0) == md.mesh.numberofvertices + 1 or np.size(field, 0) == md.mesh.numberofelements + 1:
+                md = md.checkmessage(options.getfieldvalue('message', "field {} should have only one column as there are md.mesh.numberofvertices lines".format(fieldname)))
+        elif field.shape[0] == md.mesh.numberofvertices + 1 or field.shape[0] == md.mesh.numberofelements + 1:
             if np.ndim(field) > 1 and not all(field[-1, :] == np.sort(field[-1, :])):
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' columns should be sorted chronologically" % fieldname))
+                md = md.checkmessage(options.getfieldvalue('message', "field {} columns should be sorted chronologically".format(fieldname)))
             if np.ndim(field) > 1 and any(field[-1, 0:-1] == field[-1, 1:]):
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' columns must not contain duplicate timesteps" % fieldname))
-        else:
-            md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have md.mesh.numberofvertices or md.mesh.numberofvertices + 1 lines" % fieldname))
+                md = md.checkmessage(options.getfieldvalue('message', "field {} columns must not contain duplicate timesteps".format(fieldname)))
+        else:
+            md = md.checkmessage(options.getfieldvalue('message', "field {} should have md.mesh.numberofvertices or md.mesh.numberofvertices + 1 lines".format(fieldname)))
 
     #Check single value forcings (size and times)
     if options.getfieldvalue('singletimeseries', 0):
-        if np.size(field, 0) == 2:
+        if field.shape[0] == 2:
             if not all(field[-1, :] == np.sort(field[-1, :])):
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' columns should be sorted chronologically" % fieldname))
+                md = md.checkmessage(options.getfieldvalue('message', "field {} columns should be sorted chronologically".format(fieldname)))
             if any(field[-1, 0:-1] == field[-1, 1:]):
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' columns must not contain duplicate timesteps" % fieldname))
-        elif np.size(field, 0) == 1:
+                md = md.checkmessage(options.getfieldvalue('message', "field {} columns must not contain duplicate timesteps".format(fieldname)))
+        elif field.shape[0] == 1:
             if np.ndim(field) > 1 and not np.size(field, 1) == 1:
-                md = md.checkmessage(options.getfieldvalue('message', "field '%s' should be either a scalar or have 2 lines" % fieldname))
-        else:
-            md = md.checkmessage(options.getfieldvalue('message', "field '%s' should have 2 lines or be a scalar" % fieldname))
+                md = md.checkmessage(options.getfieldvalue('message', "field {} should be either a scalar or have 2 lines".format(fieldname)))
+        else:
+            md = md.checkmessage(options.getfieldvalue('message', "field {} should have 2 lines or be a scalar".format(fieldname)))
 
     return md
Index: /issm/trunk-jpl/src/m/geometry/GetAreas3DTria.m
===================================================================
--- /issm/trunk-jpl/src/m/geometry/GetAreas3DTria.m	(revision 25498)
+++ /issm/trunk-jpl/src/m/geometry/GetAreas3DTria.m	(revision 25499)
@@ -2,29 +2,30 @@
 %GETAREAS3DTRIA - compute areas of triangles with 3D coordinates 
 %
-%   compute areas of trianguls with 3D coordinates 
+%   Compute areas of triangles with 3D coordinates 
 %
 %   Usage:
-%      areas  =GetAreas3DTria(index,x,y,z);
+%      areas=GetAreas3DTria(index,x,y,z);
 %
 %   Examples:
-%      areas  =GetAreas3DTria(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z);
+%      areas=GetAreas3DTria(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z);
 
 %get number of elements and number of nodes
-nels=size(index,1); 
-nods=length(x);  
+nels=size(index,1);
+nods=length(x);
 
 %some checks
 if nargout~=1 | (nargin~=3 & nargin~=4),
 	help GetAreas3DTria
-	error('GetAreas error message: bad usage')
+	error('GetAreas3DTria error message: bad usage')
 end
 if ((length(y)~=nods) | (nargin==4 & length(z)~=nods)),
-	error('GetAreas3DTria error message: x,y and z do not have the same length')
+	error('GetAreas3DTria error message: x, y, and z do not have the same length')
 end
+
 if max(index(:))>nods,
 	error(['GetAreas3DTria error message: index should not have values above ' num2str(nods) ])
 end
 if (nargin==4 & size(index,2)~=3),
-	error('GetAreas3DTria error message: index should have 3 columns for 2d meshes.')
+	error('GetAreas3DTria error message: index should have 3 columns for 2d meshes')
 end
 
@@ -37,12 +38,10 @@
 %compute the volume of each element
 if nargin==4,
-   % area of triangles with 3D coordinats
-   for j=1:nels
-		m1=[x1(j) x2(j) x3(j); y1(j) y2(j) y3(j); 1 1 1];
-	   m2=[y1(j) y2(j) y3(j); z1(j) z2(j) z3(j); 1 1 1];
-	   m3=[z1(j) z2(j) z3(j); x1(j) x2(j) x3(j); 1 1 1];
-      areas(j)=sqrt(det(m1)^2 + det(m2)^2 + det(m3)^2)/2;
-   end
-end 
-
-
+	% area of triangles with 3D coordinates
+	for i=1:nels
+		m1=[x1(i) x2(i) x3(i); y1(i) y2(i) y3(i); 1 1 1];
+		m2=[y1(i) y2(i) y3(i); z1(i) z2(i) z3(i); 1 1 1];
+		m3=[z1(i) z2(i) z3(i); x1(i) x2(i) x3(i); 1 1 1];
+		areas(i)=sqrt(det(m1)^2 + det(m2)^2 + det(m3)^2)/2;
+	end
+end
Index: /issm/trunk-jpl/src/m/geometry/GetAreas3DTria.py
===================================================================
--- /issm/trunk-jpl/src/m/geometry/GetAreas3DTria.py	(revision 25499)
+++ /issm/trunk-jpl/src/m/geometry/GetAreas3DTria.py	(revision 25499)
@@ -0,0 +1,62 @@
+import numpy as np
+
+
+def GetAreas3DTria(index, x, y, z, *args):
+    """GETAREAS3DTRIA - compute areas of triangles with 3D coordinates
+
+    Compute areas of triangles with 3D coordinates.
+
+    Usage:
+        areas = GetAreas3DTria(index, x, y, z)
+
+    Examples:
+        areas = GetAreas3DTria(md.mesh.elements, md.mesh.x, md.mesh.y, md.mesh.z)
+
+    TODO:
+    - Determine if *args is needed.
+    """
+
+    # Get number of elements and number of nodes
+    nels = index.shape[0]
+    nods = len(x)
+
+    # Some checks
+    nargs = len(args)
+
+    # TODO: Do we really need this under Python (first 4 arguments are required)?
+    # if nargs != 3 and nargs != 4:
+    #     print(GetAreas3DTria.__doc__)
+    #     raise Exception('GetAreas3DTria error message: bad usage')
+
+    if len(y) != nods or (nargs == 4 and len(z) != nods):
+        print(GetAreas3DTria.__doc__)
+        raise Exception('GetAreas3DTria error message: x, y, and z do not have the same length')
+
+    if np.max(index) > nods:
+        print(GetAreas3DTria.__doc__)
+        raise Exception('GetAreas3DTria error message: index should not have values above {}'.format(nods))
+
+    if nargs == 4 and index.shape[1] != 3:
+        print(GetAreas3DTria.__doc__)
+        raise Exception('GetAreas3DTria error message: index should have 3 columns for 2d meshes')
+
+    # Initialization
+    areas = np.zeros((nels, ))
+    x1 = x[index[:, 0] - 1]
+    x2 = x[index[:, 1] - 1]
+    x3 = x[index[:, 2] - 1]
+    y1 = y[index[:, 0] - 1]
+    y2 = y[index[:, 1] - 1]
+    y3 = y[index[:, 2] - 1]
+    z1 = z[index[:, 0] - 1]
+    z2 = z[index[:, 1] - 1]
+    z3 = z[index[:, 2] - 1]
+
+    # Area of triangles with 3D coordinates
+    for i in range(nels):
+        m1 = np.vstack(([x1[i], x2[i], x3[i]], [y1[i], y2[i], y3[i]], [1, 1, 1]))
+        m2 = np.vstack(([y1[i], y2[i], y3[i]], [z1[i], z2[i], z3[i]], [1, 1, 1]))
+        m3 = np.vstack(([z1[i], z2[i], z3[i]], [x1[i], x2[i], x3[i]], [1, 1, 1]))
+        areas[i] = ((np.linalg.det(m1) ** 2 + np.linalg.det(m2) ** 2 + np.linalg.det(m3) ** 2) ** 0.5) / 2 # NOTE: math.sqrt cannot be applied element-wise to a list/numpy.array
+
+    return areas
Index: /issm/trunk-jpl/src/m/interp/averaging.m
===================================================================
--- /issm/trunk-jpl/src/m/interp/averaging.m	(revision 25498)
+++ /issm/trunk-jpl/src/m/interp/averaging.m	(revision 25499)
@@ -43,5 +43,5 @@
 end
 
-%load some variables (it is much faster if the variabes are loaded from md once for all)
+%load some variables (it is much faster if the variables are loaded from md once for all)
 if layer==0,
 	index=md.mesh.elements;
@@ -87,4 +87,4 @@
 end
 
-%return output as a full matrix (C code do not like sparse matrices)
+%return output as a full matrix (C code does not like sparse matrices)
 average=full(average_node);
Index: /issm/trunk-jpl/src/m/interp/averaging.py
===================================================================
--- /issm/trunk-jpl/src/m/interp/averaging.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/interp/averaging.py	(revision 25499)
@@ -1,4 +1,3 @@
 import numpy as np
-from GetAreas import GetAreas
 try:
     from scipy.sparse import csc_matrix
@@ -6,26 +5,27 @@
     print("could not import scipy, no averaging capabilities enabled")
 
+from GetAreas import GetAreas
+
 
 def averaging(md, data, iterations, layer=0):
-    '''
-    AVERAGING - smooths the input over the mesh
+    """AVERAGING - smooths the input over the mesh
 
-       This routine takes a list over the elements or the nodes in input
-       and return a list over the nodes.
-       For each iterations it computes the average over each element (average
-       of the vertices values) and then computes the average over each node
-       by taking the average of the element around a node weighted by the
-       elements volume
-       For 3d mesh, a last argument can be added to specify the layer to be averaged on.
+    This routine takes a list of the elements or the nodes in input and return 
+    a list over the nodes.
+    For each iterations it computes the average over each element (average of 
+    the vertices values) and then computes the average over each node by taking 
+    the average of the element around a node weighted by the elements volume.
+    For 3d mesh, a last argument can be added to specify the layer to be 
+    averaged on.
 
-       Usage:
-          smoothdata = averaging(md, data, iterations)
-          smoothdata = averaging(md, data, iterations, layer)
+    Usage:
+        smoothdata = averaging(md, data, iterations)
+        smoothdata = averaging(md, data, iterations, layer)
 
-       Examples:
-          velsmoothed = averaging(md, md.initialization.vel, 4)
-          pressure = averaging(md, md.initialization.pressure, 0)
-          temperature = averaging(md, md.initialization.temperature, 1, 1)
-    '''
+    Examples:
+        velsmoothed = averaging(md, md.initialization.vel, 4)
+        pressure = averaging(md, md.initialization.pressure, 0)
+        temperature = averaging(md, md.initialization.temperature, 1, 1)
+    """
 
     if len(data) != md.mesh.numberofelements and len(data) != md.mesh.numberofvertices:
@@ -37,5 +37,5 @@
         layer = 0
 
-    #initialization
+    # Initialization
     if layer == 0:
         weights = np.zeros(md.mesh.numberofvertices, )
@@ -45,5 +45,5 @@
         data = data[(layer - 1) * md.mesh.numberofvertices2d + 1:layer * md.mesh.numberofvertices2d, :]
 
-    #load some variables (it is much faster if the variabes are loaded from md once for all)
+    # Load some variables (it is much faster if the variables are loaded from md once for all)
     if layer == 0:
         index = md.mesh.elements
@@ -55,5 +55,5 @@
         numberofelements = md.mesh.numberofelements2d
 
-    #build some variables
+    # Build some variables
     if md.mesh.dimension() == 3 and layer == 0:
         rep = 6
@@ -66,5 +66,5 @@
         areas = GetAreas(index, md.mesh.x2d, md.mesh.y2d)
 
-    index = index - 1  # since python indexes starting from zero
+    index = index - 1  # Python indexes from zero
     line = index.flatten(1)
     areas = np.vstack(areas).reshape(-1, )
@@ -72,8 +72,8 @@
     linesize = rep * numberofelements
 
-    #update weights that holds the volume of all the element holding the node i
+    # Update weights that holds the volume of all the element holding the node i
     weights = csc_matrix((np.tile(areas, (rep, 1)).reshape(-1, ), (line, np.zeros(linesize, ))), shape=(numberofnodes, 1))
 
-    #initialization
+    # Initialization
     if len(data) == numberofelements:
         average_node = csc_matrix((np.tile(areas * data, (rep, 1)).reshape(-1, ), (line, np.zeros(linesize, ))), shape=(numberofnodes, 1))
@@ -83,5 +83,5 @@
         average_node = csc_matrix(data.reshape(-1, 1))
 
-    #loop over iteration
+    # Loop over iteration
     for i in np.arange(1, iterations + 1):
         average_el = np.asarray(np.dot(average_node.todense()[index].reshape(numberofelements, rep), np.vstack(summation))).reshape(-1, )
@@ -90,5 +90,5 @@
         average_node = csc_matrix(average_node)
 
-    #return output as a full matrix (C code do not like sparse matrices)
+    # Return output as a full matrix (C code does not like sparse matrices)
     average = np.asarray(average_node.todense()).reshape(-1, )
 
Index: /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.m	(revision 25498)
+++ /issm/trunk-jpl/src/m/mesh/TwoDToThreeD.m	(revision 25499)
@@ -1,4 +1,3 @@
 function md=TwoDToThreeD(md,planet)
-
 	%reproject model into lat,long if necessary:
 	if ~strcmpi(md.mesh.proj,epsg2proj(4326)),
Index: /issm/trunk-jpl/src/m/mesh/findsegments.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/findsegments.m	(revision 25498)
+++ /issm/trunk-jpl/src/m/mesh/findsegments.m	(revision 25499)
@@ -66,10 +66,10 @@
 	%'el1' is connected to only one element
 	else
-
-		%find the vertex the 'el1' does not share with 'els2'
+		%find the vertex that 'el1' does not share with 'els2'
 		flag=setdiff(nods1,md.mesh.elements(els2,:));
 
 		for j=1:3,
-			nods=nods1; nods(j)=[];
+			nods=nods1;
+			nods(j)=[];
 			if any(ismember(flag,nods)),
 
Index: /issm/trunk-jpl/src/m/mesh/findsegments.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/findsegments.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/mesh/findsegments.py	(revision 25499)
@@ -35,5 +35,5 @@
     pos = np.nonzero(elementonboundary)[0]
     num_segments = len(pos)
-    segments = np.zeros((num_segments, 3))
+    segments = np.zeros((num_segments, 3)).astype(int)
     count = 0
 
@@ -56,5 +56,5 @@
 
             # Get the vertices on the boundary and build segment
-            nods1 = np.delete(nods1, np.where(nods1 == flag)[0])
+            nods1 = np.delete(nods1, np.where(np.in1d(nods1, flag, assume_unique=True)))
             segments[count, :] = np.append(nods1, el1 + 1)
 
@@ -72,8 +72,26 @@
         # 'el1' is connected to only one element
         else:
+            # NOTE: This block is untested as it does not get touched by test2004 (remove this note once it has been tested)
+
             # Find the vertex that 'el1' does not share with 'els2'
-            flag = els2
-        
-        exit()
+            flag = np.setdiff1d(nods, md.mesh.elements[els2, :])
+
+            for j in range(3):
+                nods = nods1
+                nods = np.delete(nods, j)
+                if np.any(np.in1d(flag, nods)):
+                    segments[count, :] = np.append(nods, el1 + 1)
+
+                    # Swap segment nodes if necessary
+                    ord1 = np.where(nods1[0] == md.mesh.elements[el1, :])[0][0]
+                    ord2 = np.where(nods1[1] == md.mesh.elements[el1, :])[0][0]
+
+                    if ((ord1 == 0 and ord2 == 1) or (ord1 == 1 and ord2 == 2) or (ord1 == 2 and ord2 == 0)):
+                        temp = segments[count, 0]
+                        segments[count, 0] = segments[count, 1]
+                        segments[count, 1] = temp
+
+                    segments[count, 0:2] = np.flip(segments[count, 0:2]) # NOTE: Upper bound of index range is non-inclusive
+                    count = count + 1
 
     return segments
Index: /issm/trunk-jpl/src/m/mesh/meshintersect3d.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/meshintersect3d.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/mesh/meshintersect3d.py	(revision 25499)
@@ -1,4 +1,2 @@
-import math
-
 import matplotlib.pyplot as plt
 from mpl_toolkits.mplot3d import Axes3D
@@ -9,8 +7,7 @@
 
 def meshintersect3d(x, y, z, xs, ys, zs, *args): #{{{
-    '''
-    MESHINTERSECT - returns indices (into x, y, and z) of common values between 
-    (x, y, z) and (xs, ys, zs) (i.e. x(index) = xs; y(index) = ys).
-    '''
+    """MESHINTERSECT - returns indices (into x, y, and z) of common values 
+    between (x, y, z) and (xs, ys, zs) (i.e. x(index) = xs; y(index) = ys).
+    """
 
     # Process options
@@ -19,18 +16,20 @@
     # Retrieve tolerance
     maxtol = options.getfieldvalue('maxtol', 100000) # 100 km
-    tolincrement = options.getfieldvalue('tolincrement', [])
+    tolincrement = options.getfieldvalue('tolincrement', 10)
     force = options.getfieldvalue('force', 0)
 
-    # Go through lats, longs and find within tolerance, the index of the corresponding value in lat, long
-    indices = np.zeros(len(xs))
+    # Go through lats, longs and find, within tolerance, the index of the corresponding value in lat, long
+    indices = np.zeros((len(xs), ))
 
     for i in range(len(xs)):
         tolerance = 0
-        distance = math.sqrt((x - xs[i]) ** 2 + (y - ys[i]) ** 2 + (z - zs[i]) ** 2)
+        distance = ((x - xs[i]) ** 2 + (y - ys[i]) ** 2 + (z - zs[i]) ** 2) ** 0.5 # NOTE: math.sqrt cannot be applied element-wise to a list/numpy.array
+        s = np.where(distance == 0)[0]
 
-        s = np.where(distance == 0)[0]
-        if s.size():
-            if s.size() > 1:
+        if len(s):
+            if len(s) > 1:
                 # We have two vertices that are coincident! Not good
+                #
+                # TODO: Reconfigure the following in the process of bringing plotting online
                 for j in range(len(s)):
                     plot(x[s[j]], y[s[j]], z[s[j]], c='cyan', s=40)
@@ -46,24 +45,25 @@
             # We could not find a 0 distance, find the lowest tolerance that generates a find
             count = 1
-            while not s.size():
+            while not len(s):
                 if count > 1000:
-                    print('could not find a vertex matching vertex %i of input mesh!' % i)
-                    print('Might think anbout changing tolerance increment')
-                    raise RuntimeError('')
+                    raise Exception('could not find a vertex matching vertex {} of input mesh!\nMight think anbout changing tolerance increment'.format(i))
                 tolerance = tolerance + tolincrement
-                s = np.where(distance < tolerance)
+                s = np.where(distance < tolerance)[0]
                 count = count + 1
             if tolerance > maxtol:
-                print('found matching vertices %i in output mesh for input mesh vertex %i' % (s, i))
-                print('however, these vertices are farther than the max tolerance allowed!')
-                raise RuntimeError('')
+                raise Exception('found matching vertices {} in output mesh for input mesh vertex {}\nhowever, these vertices are farther than the max tolerance allowed!'.format(s, i))
 
             # Recover minimum distance
             sf = distance[s]
-            pos = np.where(sf == sf.min())
+            pos = np.where(sf == np.min(sf))[0]
             s = s[pos]
             indices[i] = s
 
-    if np.where(indices == 0).size():
+    if len(np.where(indices == 0)[0]) > 1: # NOTE: This check is different than the corresponding one under MATLAB as one index may indeed be '0'
         raise RuntimeError('issue with transition vector having one empty slot')
+
+    # Convert results to type 'int' to avoid modifying structures to which they are assigned
+    indices = indices.astype(int)
+
+    return indices
 #}}}
Index: /issm/trunk-jpl/src/m/mesh/modelmerge3d.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/modelmerge3d.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/mesh/modelmerge3d.py	(revision 25499)
@@ -112,7 +112,6 @@
 
     # Vertex on boundary
-    md.mesh.vertexconnectivity = np.zeros((md.mesh.numberofvertices, ))
-    md.mesh.vertexonboundary[md.mesh.segments[:, 0:1]] = 1
-    print(md.mesh.vertexonboundary)
+    md.mesh.vertexonboundary = np.zeros((md.mesh.numberofvertices, ))
+    md.mesh.vertexonboundary[md.mesh.segments[:, 0:2] - 1] = 1
 
     # Some checks
@@ -120,5 +119,5 @@
         raise Exception('issue in modelmerge, one of the element ids is > number of vertices!')
 
-    exit()
+    return md
 
 
Index: /issm/trunk-jpl/src/m/miscellaneous/getsubattr.py
===================================================================
--- /issm/trunk-jpl/src/m/miscellaneous/getsubattr.py	(revision 25499)
+++ /issm/trunk-jpl/src/m/miscellaneous/getsubattr.py	(revision 25499)
@@ -0,0 +1,27 @@
+def getsubattr(obj, string):
+    """GETSUBATTR - Returns a reference to the desired attribute of 'obj' based 
+    on 'string'.
+
+    If 'string' represents a single attribute, this function works just like 
+    'getattr'.
+
+    Usage:
+        attr = getsubattr(obj, string)
+
+        where 'object' is a structure that can be addressed with dot notation 
+        and 'string' is a string with one or more attributes separated by '.'.
+
+    Example:
+        attr = getsubattr(md, 'mask.land_levelset')
+        attr = getsubattr(md, 'mask') # Works just like getattr(md, 'mask')
+    """
+
+    attrs = string.split('.')
+
+    # Recurse until we get desired attribute
+    if len(attrs) > 1:
+        attr = getsubattr(getattr(obj, attrs[0]), '.'.join(attrs[1:]))
+    else:
+        attr = getattr(obj, string)
+    
+    return attr
Index: /issm/trunk-jpl/src/m/miscellaneous/setsubattr.py
===================================================================
--- /issm/trunk-jpl/src/m/miscellaneous/setsubattr.py	(revision 25499)
+++ /issm/trunk-jpl/src/m/miscellaneous/setsubattr.py	(revision 25499)
@@ -0,0 +1,26 @@
+def setsubattr(obj, string, value):
+    """SETSUBATTR - Sets the desired attribute of 'obj' based on 'string' to 
+    the value supplied to 'value'.
+
+    If 'string' represents a single attribute, this function works just like 
+    'setattr'.
+
+    Usage:
+        setsubattr(obj, string, value)
+
+        where 'object' is a structure that can be addressed with dot notation, 
+        'string' is a string with one or more attributes separated by '.', and 
+        'value' is any value.
+
+    Example:
+        attr = getsubattr(md, 'mask.land_levelset')
+        attr = getsubattr(md, 'mask') # Works just like getattr(md, 'mask')
+    """
+
+    attrs = string.split('.')
+
+    # Recurse until we get desired attribute
+    if len(attrs) > 1:
+        setsubattr(getattr(obj, attrs[0]), '.'.join(attrs[1:]), value)
+    else:
+        setattr(obj, string, value)
Index: /issm/trunk-jpl/src/m/parameterization/contourenvelope.m
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/contourenvelope.m	(revision 25498)
+++ /issm/trunk-jpl/src/m/parameterization/contourenvelope.m	(revision 25499)
@@ -27,5 +27,5 @@
 		isfile=0;
 	else
-		error('contourenvelope error message:  second argument should be a file or an elements flag');
+		error('contourenvelope error message: second argument should be a file or an elements flag');
 	end
 end
Index: /issm/trunk-jpl/src/m/parameterization/contourenvelope.py
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/contourenvelope.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/parameterization/contourenvelope.py	(revision 25499)
@@ -1,27 +1,31 @@
+import copy
 import os.path
 import numpy as np
-import copy
+
+from ContourToMesh import ContourToMesh
+from ElementConnectivity import ElementConnectivity
+import MatlabFuncs as m
 from NodeConnectivity import NodeConnectivity
-from ElementConnectivity import ElementConnectivity
-from ContourToMesh import ContourToMesh
-import MatlabFuncs as m
 
 
-def contourenvelope(md, *args):
+def contourenvelope(mh, *args):
     """CONTOURENVELOPE - build a set of segments enveloping a contour .exp
 
     Usage:
-        segments = contourenvelope(md, varargin)
+        segments = contourenvelope(mh, *args)
 
     Example:
-        segments = contourenvelope(md, 'Stream.exp')
-        segments = contourenvelope(md)
+        segments = contourenvelope(mh, 'Stream.exp')
+        segments = contourenvelope(mh)
     """
 
-    #some checks
-    if len(args) > 1:
-        raise RuntimeError("contourenvelope error message: bad usage")
+    # Some checks
+    nargs = len(args)
 
-    if len(args) == 1:
+    if nargs > 1:
+        print(contourenvelope.__doc__)
+        raise Exception("contourenvelope error message: bad usage")
+
+    if nargs == 1:
         flags = args[0]
 
@@ -35,40 +39,39 @@
             isfile = 0
         else:
-            raise TypeError("contourenvelope error message:  second argument should be a file or an elements flag")
+            raise TypeError("contourenvelope error message: second argument should be a file or an elements flag")
 
-    #Now, build the connectivity tables for this mesh.
-    #Computing connectivity
-    if np.size(md.mesh.vertexconnectivity, axis=0) != md.mesh.numberofvertices and np.size(md.mesh.vertexconnectivity, axis=0) != md.mesh.numberofvertices2d:
-        md.mesh.vertexconnectivity = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
-    if np.size(md.mesh.elementconnectivity, axis=0) != md.mesh.numberofelements and np.size(md.mesh.elementconnectivity, axis=0) != md.mesh.numberofelements2d:
-        md.mesh.elementconnectivity = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
+    # Now, build the connectivity tables for this mesh
+    # Computing connectivity
+    if np.size(mh.vertexconnectivity, axis=0) != mh.numberofvertices and np.size(mh.vertexconnectivity, axis=0) != mh.numberofvertices2d:
+        mh.vertexconnectivity = NodeConnectivity(mh.elements, mh.numberofvertices)
+    if np.size(mh.elementconnectivity, axis=0) != mh.numberofelements and np.size(mh.elementconnectivity, axis=0) != mh.numberofelements2d:
+        mh.elementconnectivity = ElementConnectivity(mh.elements, mh.vertexconnectivity)
 
     #get nodes inside profile
-    elementconnectivity = copy.deepcopy(md.mesh.elementconnectivity)
-    if md.mesh.dimension() == 2:
-        elements = copy.deepcopy(md.mesh.elements)
-        x = copy.deepcopy(md.mesh.x)
-        y = copy.deepcopy(md.mesh.y)
-        numberofvertices = copy.deepcopy(md.mesh.numberofvertices)
-        numberofelements = copy.deepcopy(md.mesh.numberofelements)
+    elementconnectivity = copy.deepcopy(mh.elementconnectivity)
+    if mh.dimension() == 2:
+        elements = copy.deepcopy(mh.elements)
+        x = copy.deepcopy(mh.x)
+        y = copy.deepcopy(mh.y)
+        numberofvertices = copy.deepcopy(mh.numberofvertices)
+        numberofelements = copy.deepcopy(mh.numberofelements)
     else:
-        elements = copy.deepcopy(md.mesh.elements2d)
-        x = copy.deepcopy(md.mesh.x2d)
-        y = copy.deepcopy(md.mesh.y2d)
-        numberofvertices = copy.deepcopy(md.mesh.numberofvertices2d)
-        numberofelements = copy.deepcopy(md.mesh.numberofelements2d)
+        elements = copy.deepcopy(mh.elements2d)
+        x = copy.deepcopy(mh.x2d)
+        y = copy.deepcopy(mh.y2d)
+        numberofvertices = copy.deepcopy(mh.numberofvertices2d)
+        numberofelements = copy.deepcopy(mh.numberofelements2d)
 
     if len(args) == 1:
-
         if isfile:
-            #get flag list of elements and nodes inside the contour
+            # Get flag list of elements and nodes inside the contour
             nodein = ContourToMesh(elements, x, y, file, 'node', 1)
             elemin = (np.sum(nodein(elements), axis=1) == np.size(elements, axis=1))
-    #modify element connectivity
+            # Modify element connectivity
             elemout = np.nonzero(np.logical_not(elemin))[0]
             elementconnectivity[elemout, :] = 0
             elementconnectivity[np.nonzero(m.ismember(elementconnectivity, elemout + 1))] = 0
         else:
-            #get flag list of elements and nodes inside the contour
+            # Get flag list of elements and nodes inside the contour
             nodein = np.zeros(numberofvertices)
             elemin = np.zeros(numberofelements)
@@ -78,11 +81,11 @@
             nodein[elements[pos, :] - 1] = 1
 
-    #modify element connectivity
+            # Modify element connectivity
             elemout = np.nonzero(np.logical_not(elemin))[0]
             elementconnectivity[elemout, :] = 0
             elementconnectivity[np.nonzero(m.ismember(elementconnectivity, elemout + 1))] = 0
 
-    #Find element on boundary
-    #First: find elements on the boundary of the domain
+    # Find element on boundary
+    # First: find elements on the boundary of the domain
     flag = copy.deepcopy(elementconnectivity)
     if len(args) == 1:
@@ -90,5 +93,5 @@
     elementonboundary = np.logical_and(np.prod(flag, axis=1) == 0, np.sum(flag, axis=1) > 0)
 
-    #Find segments on boundary
+    # Find segments on boundary
     pos = np.nonzero(elementonboundary)[0]
     num_segments = np.size(pos)
Index: /issm/trunk-jpl/src/m/solve/marshall.m
===================================================================
--- /issm/trunk-jpl/src/m/solve/marshall.m	(revision 25498)
+++ /issm/trunk-jpl/src/m/solve/marshall.m	(revision 25499)
@@ -45,5 +45,5 @@
 % % Uncomment the following to make a copy of the binary input file for debugging 
 % % purposes (can be fed into scripts/ReadBin.py).
-% copyfile([md.miscellaneous.name '.bin'], [md.miscellaneous.name '.m.bin'])
+copyfile([md.miscellaneous.name '.bin'], [md.miscellaneous.name '.m.bin'])
 if st==-1,
 	error(['marshall error message: could not close file ' [md.miscellaneous.name '.bin']]);
Index: /issm/trunk-jpl/src/m/solve/marshall.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/marshall.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/solve/marshall.py	(revision 25499)
@@ -44,6 +44,6 @@
         # # Uncomment the following to make a copy of the binary input file for 
         # # debugging purposes (can be fed into scripts/ReadBin.py).
-        # copy_cmd = "cp {}.bin {}.py.bin".format(md.miscellaneous.name, md.miscellaneous.name)
-        # subprocess.call(copy_cmd, shell=True)
+        copy_cmd = "cp {}.bin {}.py.bin".format(md.miscellaneous.name, md.miscellaneous.name)
+        subprocess.call(copy_cmd, shell=True)
     except IOError as e:
         print("marshall error message: could not close file '{}.bin' due to:".format(md.miscellaneous.name), e)
Index: /issm/trunk-jpl/src/m/solve/solve.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/solve.py	(revision 25498)
+++ /issm/trunk-jpl/src/m/solve/solve.py	(revision 25499)
@@ -12,10 +12,9 @@
 
 def solve(md, solutionstring, *args):
-    '''
-    SOLVE - apply solution sequence for this model
+    """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:
@@ -43,8 +42,8 @@
                                   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')
+    """
 
     #recover and process solve options
Index: /issm/trunk-jpl/test/NightlyRun/test2004.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2004.m	(revision 25498)
+++ /issm/trunk-jpl/test/NightlyRun/test2004.m	(revision 25499)
@@ -93,5 +93,6 @@
 
 	%miscellaneous: 
-	md.mesh.proj=bas.proj; md.miscellaneous.name=bas.name;
+	md.mesh.proj=bas.proj;
+	md.miscellaneous.name=bas.name;
 
 	%recover mask where we have land: 
@@ -324,5 +325,5 @@
 sl.caticecaps('tolerance',tolerance,'loneedgesdetect',loneedgesdetect);
 
-%figure out how each icecap's mesh connects to the larger earth mesh: 
+%figure out how each icecap's mesh connects to the larger Earth mesh: 
 sl.intersections('force',1);