Index: /issm/trunk-jpl/src/m/classes/dsl.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/dsl.m	(revision 26065)
+++ /issm/trunk-jpl/src/m/classes/dsl.m	(revision 26066)
@@ -65,5 +65,5 @@
 		end % }}}
 		function savemodeljs(self,fid,modelname) % {{{
-		
+
 			writejs1Darray(fid,[modelname '.dsl.global_average_thermosteric_sea_level'],self.global_average_thermosteric_sea_level);
 			writejs1Darray(fid,[modelname '.dsl.sea_surface_height_above_geoid'],self.sea_surface_height_above_geoid);
@@ -85,5 +85,4 @@
 				disp('      no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero');
 			end
-
 		end % }}}
 	
Index: /issm/trunk-jpl/src/m/classes/dsl.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/dsl.py	(revision 26065)
+++ /issm/trunk-jpl/src/m/classes/dsl.py	(revision 26066)
@@ -28,12 +28,20 @@
     #}}}
 
-    def extrude(self, md): #{{{
-        self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node')
-        self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node')
-        return self
+    def defaultoutputs(self, md): #{{{
+        return []
     #}}}
 
-    def defaultoutputs(self, md): #{{{
-        return []
+    def initialize(self, md): #{{{
+        if np.isnan(self.global_average_thermosteric_sea_level):
+            self.global_average_thermosteric_sea_level = np.array([0, 0]).reshape(-1, 1)
+            print('      no dsl.global_average_thermosteric_sea_level specified: transient values set at zero')
+
+        if np.isnan(self.sea_surface_height_above_geoid):
+            self.sea_surface_height_above_geoid = np.array(np.zeros(md.mesh.numberofvertices)).reshape(-1, 1)
+            print('      no dsl.sea_surface_height_above_geoid specified: transient values set at zero')
+
+        if np.isnan(self.sea_water_pressure_at_sea_floor):
+            self.sea_water_pressure_at_sea_floor = np.array(np.zeros(md.mesh.numberofvertices)).reshape(-1, 1)
+            print('      no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero')
     #}}}
 
@@ -53,4 +61,10 @@
     # }}}
 
+    def extrude(self, md): #{{{
+        self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node')
+        self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node')
+        return self
+    #}}}
+
     def marshall(self, prefix, md, fid):   #{{{
         yts = md.constants.yts
Index: /issm/trunk-jpl/src/m/parameterization/parameterize.py
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/parameterize.py	(revision 26065)
+++ /issm/trunk-jpl/src/m/parameterization/parameterize.py	(revision 26066)
@@ -4,17 +4,16 @@
 
 def parameterize(md, parametername):
-    '''
-    PARAMETERIZE - parameterize a model
+    """PARAMETERIZE - parameterize a model
 
-        From a parameter Python file, start filling in all the model fields 
-        that were not filled in by the mesh.py and mask.py model methods.
-        Warning: the parameter file must be able to be run in Python
+    From a parameter Python file, start filling in all the model fields that 
+    were not filled in by the mesh.py and mask.py model methods. Warning: the 
+    parameter file must be able to be run in Python
 
-        Usage:
-            md = parameterize(md, parametername)
+    Usage:
+        md = parameterize(md, parametername)
 
-        Example:
-            md = parameterize(md, 'Square.par')
-    '''
+    Example:
+        md = parameterize(md, 'Square.py')
+    """
 
     #some checks
Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.par
===================================================================
--- /issm/trunk-jpl/test/Par/SquareSheetConstrained.par	(revision 26065)
+++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.par	(revision 26066)
@@ -60,5 +60,4 @@
 md.groundingline.migration='None';
 
-
 %Boundary conditions:
 md=SetIceSheetBC(md);
Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareSheetConstrained.py	(revision 26065)
+++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.py	(revision 26066)
@@ -19,4 +19,5 @@
 md.geometry.thickness = hmax + (hmin - hmax) * (md.mesh.y - ymin) / (ymax - ymin) + 0.1 * (hmin - hmax) * (md.mesh.x - xmin) / (xmax - xmin)
 md.geometry.base = -md.materials.rho_ice / md.materials.rho_water * md.geometry.thickness + 20.
+md.geometry.bed = md.geometry.base
 md.geometry.surface = md.geometry.base + md.geometry.thickness
 
@@ -61,12 +62,4 @@
 md.groundingline.migration = 'None'
 
-#GIA:
-md.gia.lithosphere_thickness = 100. * np.ones((md.mesh.numberofvertices))  # in km
-md.gia.mantle_viscosity = 1.e21 * np.ones((md.mesh.numberofvertices))  # in Pa.s
-md.materials.lithosphere_shear_modulus = 6.7e10  # in Pa
-md.materials.lithosphere_density = 3.32  # in g / cm^ - 3
-md.materials.mantle_shear_modulus = 1.45e11  # in Pa
-md.materials.mantle_density = 3.34  # in g / cm^ - 3
-
 #Boundary conditions:
 md = SetIceSheetBC(md)