Index: /issm/trunk-jpl/src/m/classes/slr.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/slr.py	(revision 25162)
+++ /issm/trunk-jpl/src/m/classes/slr.py	(revision 25163)
@@ -10,14 +10,14 @@
 
 class slr(object):
-    """
+    '''
     SLR class definition
 
         Usage:
           slr = slr()
-    """
+    '''
     def __init__(self):  # {{{
-        self.deltathickness = float('NaN')
-        self.sealevel = float('NaN')
-        self.spcthickness = float('NaN')
+        self.deltathickness = np.nan
+        self.sealevel = np.nan
+        self.spcthickness = np.nan
         self.maxiter = 0
         self.reltol = 0
@@ -42,6 +42,4 @@
         self.degacc = 0
         self.horiz = 0
-        self.Ngia = float('NaN')
-        self.Ugia = float('NaN')
         self.planetradius = planetradius('earth')
         self.requested_outputs = []
@@ -88,6 +86,6 @@
     def setdefaultparameters(self):  # {{{
         #Convergence criterion: absolute, relative and residual
-        self.reltol = 0.01  #default
-        self.abstol = float('NaN')  #1 mm of sea level rise
+        self.reltol = 0.01 #default
+        self.abstol = np.nan #1 mm of sea level rise
         #maximum of non - linear iterations.
         self.maxiter = 5
@@ -110,6 +108,5 @@
         #numerical discretization accuracy
         self.degacc = 0.01
-        #steric:
-        self.steric_rate = 0
+        #hydro
         self.hydro_rate = 0
         #how many time steps we skip before we run SLR solver during transient
@@ -148,11 +145,8 @@
         md = checkfield(md, 'fieldname', 'slr.maxiter', 'size', [1, 1], '>=', 1)
         md = checkfield(md, 'fieldname', 'slr.geodetic_run_frequency', 'size', [1, 1], '>=', 1)
-        md = checkfield(md, 'fieldname', 'slr.steric_rate', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
         md = checkfield(md, 'fieldname', 'slr.hydro_rate', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
         md = checkfield(md, 'fieldname', 'slr.degacc', 'size', [1, 1], '>=', 1e-10)
         md = checkfield(md, 'fieldname', 'slr.requested_outputs', 'stringrow', 1)
         md = checkfield(md, 'fieldname', 'slr.horiz', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
-        md = checkfield(md, 'fieldname', 'slr.Ngia', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
-        md = checkfield(md, 'fieldname', 'slr.Ugia', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
 
         #check that love numbers are provided at the same level of accuracy:
@@ -199,8 +193,5 @@
         WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'format', 'Boolean')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'geodetic_run_frequency', 'format', 'Integer')
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'steric_rate', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1e-3 / md.constants.yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'hydro_rate', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1e-3 / md.constants.yts)
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'Ngia', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1e-3 / md.constants.yts)
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'Ugia', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1e-3 / md.constants.yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'format', 'Double')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'transitions', 'format', 'MatArray')
Index: /issm/trunk-jpl/src/m/coordsystems/epsg2proj.py
===================================================================
--- /issm/trunk-jpl/src/m/coordsystems/epsg2proj.py	(revision 25162)
+++ /issm/trunk-jpl/src/m/coordsystems/epsg2proj.py	(revision 25163)
@@ -21,5 +21,5 @@
 
     #First, get GDAL version
-    subproc_args = "gdalsrsinfo --version | awk '{print $2}' | cut -d '.' -f1"
+    subproc_args = shlex.split("gdalsrsinfo --version | awk '{print $2}' | cut -d '.' -f1")
     subproc = subprocess.Popen(subproc_args, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
     outs, errs = subproc.communicate()
@@ -29,5 +29,5 @@
     version_major=int(outs)
 
-    subproc_args = "gdalsrsinfo epsg:{} | grep PROJ.4 | tr -d '\n' | sed 's/PROJ.4 : //'".format(epsg)
+    subproc_args = shlex.split("gdalsrsinfo epsg:{} | grep PROJ.4 | tr -d '\n' | sed 's/PROJ.4 : //'".format(epsg))
     subproc = subprocess.Popen(subproc_args, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
     outs, errs = subproc.communicate()
Index: /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py
===================================================================
--- /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py	(revision 25162)
+++ /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py	(revision 25163)
@@ -34,5 +34,5 @@
 
     TODO:
-    - Remove shlex and follow subprocess pattern implemented in src/m/coordsystems/epsg2proj.py
+    - Follow subprocess pattern implemented in src/m/coordsystems/epsg2proj.py
     '''
 
@@ -46,6 +46,6 @@
     file_in.close()
 
-    args = shlex.split('gdaltransform -s_srs {} -t_srs {} < {} > {}'.format(proj_in, proj_out, filename_in, filename_out))
-    subprocess.check_call(args) # Will raise CalledProcessError if return code is not 0
+    subproc_args = shlex.split('gdaltransform -s_srs {} -t_srs {} < {} > {}'.format(proj_in, proj_out, filename_in, filename_out))
+    subprocess.check_call(subproc_args) # Will raise CalledProcessError if return code is not 0
 
     A = loadvars(filename_out)
Index: /issm/trunk-jpl/src/m/modules/ExpToLevelSet.py
===================================================================
--- /issm/trunk-jpl/src/m/modules/ExpToLevelSet.py	(revision 25163)
+++ /issm/trunk-jpl/src/m/modules/ExpToLevelSet.py	(revision 25163)
@@ -0,0 +1,36 @@
+from helpers import fileparts
+from shpread import shpread
+
+def ExpToLevelSet(x, y, contourname): #{{{
+    '''
+    EXPTOLEVELSET - Determine levelset distance between a contour and a cloud 
+    of points
+
+        Usage:
+            distance = ExpToLevelSet(x, y, contourname)
+
+        x, y:           cloud point
+        contourname:    name of .exp file containing the contours
+        distance:       distance vector representing a levelset where the 0 
+                        level is one of the contour segments
+
+        Example:
+            distance = ExpToLevelSet(md.mesh.x, md.mesh.y, 'Contour.exp')
+
+        TODO:
+        - Need to compile Python version of ExpToLevelSet_matlab for this 
+        to work as intended (see src/m/modules/ExpToLevelSet.m)
+    '''
+
+    if isinstance(contourname, basestring):
+        path, name, ext = fileparts(contourname)
+        if ext == '.shp':
+            #read contour from shapefile
+            contourname = shpread(contourname)
+
+    # NOTE: This library does not currently exist! See TODO list in function 
+    #       header.
+    distance = ExpToLevelSet_python(x, y, contourname)
+
+    return distance
+#}}}
Index: /issm/trunk-jpl/src/m/partition/AreaAverageOntoPartition.py
===================================================================
--- /issm/trunk-jpl/src/m/partition/AreaAverageOntoPartition.py	(revision 25162)
+++ /issm/trunk-jpl/src/m/partition/AreaAverageOntoPartition.py	(revision 25163)
@@ -1,6 +1,7 @@
+import copy
+
 import numpy as np
 
 from adjacency import adjacency
-import copy
 from project2d import project2d
 from qmupart2npart import qmupart2npart
Index: /issm/trunk-jpl/src/m/plot/processmesh.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/processmesh.m	(revision 25162)
+++ /issm/trunk-jpl/src/m/plot/processmesh.m	(revision 25163)
@@ -15,5 +15,5 @@
 end
 
-%special case for mesg 2dvertical
+%special case for mesh 2dvertical
 if strcmp(domaintype(md.mesh),'2Dvertical'),
 	[x y z elements is2d isplanet] = processmesh(md.mesh,options);
@@ -64,12 +64,8 @@
 
 %is it a 2d plot?
-if md.mesh.dimension()==2,
+if md.mesh.dimension()==2 | getfieldvalue(options,'layer',0)>=1,
 	is2d=1;
 else
-	if getfieldvalue(options,'layer',0)>=1,
-		is2d=1;
-	else
-		is2d=0;
-	end
+	is2d=0;
 end
 
@@ -79,5 +75,5 @@
 		error('processmesh error message: cannot work with 3D meshes for now');
 	end
-	%we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded. 
+	%we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded
 	x=x2d;
 	y=y2d;
Index: /issm/trunk-jpl/src/m/plot/processmesh.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/processmesh.py	(revision 25162)
+++ /issm/trunk-jpl/src/m/plot/processmesh.py	(revision 25163)
@@ -1,40 +1,68 @@
+#import matplotlib.delaunay as delaunay
 import numpy as np
 
 
 def processmesh(md, data, options):
-    """
-    PROCESSMESH - process the mesh for plotting
+    '''
+    PROCESSMESH - process mesh to be plotted
 
     Usage:
-    x, y, z, elements, is2d = processmech(md, data, options)
+        x, y, z, elements, is2d = processmech(md, data, options)
 
     See also: PLOTMODEL, PROCESSDATA
-    """
 
-    # {{{ check mesh size parameters
+    TODO:
+    - Test application of matlplotlib.delaunay
+    - Test that output of delaunay matches output of delaunay in MATLAB (see 
+    src/m/plot/processmesh.m)
+    '''
+
+    #some checks
     if md.mesh.numberofvertices == 0:
         raise ValueError('processmesh error: mesh is empty')
     if md.mesh.numberofvertices == md.mesh.numberofelements:
         raise ValueError('processmesh error: the number of elements is the same as the number of nodes')
-    # }}}
-    # {{{ treating coordinates
 
-    try:
+    #special case for mesh 2dvertical
+    if md.mesh.domaintype() == '2Dvertical':
+        x, y, z, elements, is2d, isplanet = processmesh(md.mesh, options)
+        return x, y, z, elements, is2d, isplanet
+
+    # # special case for mesh 3dsurface
+    # if md.mesh.domaintype() == '3Dsurface':
+    #     x, y, z, elements, is2d, isplanet = processmesh(md.mesh, options)
+    #     if options.getfieldvalue('coord', 'xy') == 'latlon':
+    #         x = md.mesh.long
+    #         y = md.mesh.lat
+    #         elements = delaunay(x, y)
+    #         z = md.mesh.lat
+    #         z[:] = 0
+    #     return x, y, z, elements, is2d, isplanet
+
+    if hasattr(md.mesh, 'elements2d'):
+        elements2d = md.mesh.elements2d
+
+    if options.exist('amr'):
+        step = options.getfieldvalue('amr')
+        x = md.results.TransientSolution[step].MeshX
+        y = md.results.TransientSolution[step].MeshY
+        elements = md.results.TransientSolution[step].MeshElements
+    else:
+        elements = md.mesh.elements
+        if options.getfieldvalue('coord', 'xy') != 'latlon':
+            x = md.mesh.x
+            if hasattr(md.mesh, 'x2d'):
+                x2d = md.mesh.x2d
+            y = md.mesh.y
+            if hasattr(md.mesh, 'y2d'):
+                y2d = md.mesh.y2d
+
+    if hasattr(md.mesh, 'z'):
         z = md.mesh.z
-    except AttributeError:
+    else:
         z = np.zeros(np.shape(md.mesh.x))
-    elements = md.mesh.elements - 1
-
-    if options.getfieldvalue('layer', 0) >= 1:
-        x = md.mesh.x2d
-        y = md.mesh.y2d
-        z = np.zeros(np.shape(x))
-        elements = md.mesh.elements2d - 1
-    elif 'latlon' in options.getfieldvalue('coord', 'xy'):
-        x = md.mesh.long
-        y = md.mesh.lat
-    else:
-        x = md.mesh.x
-        y = md.mesh.y
+    z = options.getfieldvalue('z', z)
+    if isinstance(z, basestring):
+        z = getattr(md, z)
 
     #is it a 2D plot?
@@ -43,4 +71,16 @@
     else:
         is2d = 0
+
+    elements = md.mesh.elements - 1
+
+    #layer projection?
+    if options.getfieldvalue('layer', 0) >= 1:
+        if options.getfieldvalue('coord', 'xy') == 'latlon':
+            raise Exception('processmesh error message: cannot work with 3D meshes for now')
+        #we modify the mesh temporarily to a 2d mesh from which the 3d mesh was extruded
+        x = x2d
+        y = y2d
+        z = np.zeros(np.shape(x2d))
+        elements = md.mesh.elements2d - 1
 
     #units
@@ -51,5 +91,8 @@
         z = z * unit
 
-    #is model a member of planet class? (workaround until planet class defined)
+    #is model a member of planet class? 
+    #
+    # TODO: Change this when planet class defined (see src/m/plot/processmesh.m)
+    #
     if md.__class__.__name__ != 'model':
         isplanet = 1
Index: /issm/trunk-jpl/src/m/qmu/dakota_in_data.py
===================================================================
--- /issm/trunk-jpl/src/m/qmu/dakota_in_data.py	(revision 25162)
+++ /issm/trunk-jpl/src/m/qmu/dakota_in_data.py	(revision 25163)
@@ -9,37 +9,34 @@
 def dakota_in_data(dmeth, variables, responses, dparams, filei, *args):
     '''
-  define the data to write the dakota .in and .m files.
+    DAKOTA_IN_DATA - define the data to write the dakota .in and .m files.
+    
+    Usage:
+        [] = dakota_in_data(dmeth, variables, responses, dparams, filei, *args)
 
-  [] = dakota_in_data(dmeth, variables, responses, dparams, filei, *args)
+    where the required input is:
+        dmeth         (dakota_method, method class object)
+        variables     (structure array, variable class objects)
+        responses     (structure array, response class objects)
+        dparams       (structure array, method - independent parameters)
+        filei         (character, name of .in and .m files)
 
-  where the required input is:
-    dmeth         (dakota_method, method class object)
-    variables     (structure array, variable class objects)
-    responses     (structure array, response class objects)
-    dparams       (structure array, method - independent parameters)
-    filei         (character, name of .in and .m files)
+    Params may be empty, in which case defaults will be used.
 
-  params may be empty, in which case defaults will be used.
+    The optional args are passed directly through to the QmuUpdateFunctions 
+    brancher to be used by the analysis package. For example, this could be 
+    model information.
 
-  the optional args are passed directly through to the
-  QmuUpdateFunctions brancher to be used by the analysis
-  package.  for example, this could be model information.
+    This function defines the data to write the dakota .in and .m files. It is 
+    necessary for multiple reasons. First, it collects the parameters and 
+    applies some defaults that are unique to the environment. Second, some 
+    analysis package ariables and/or responses may be treated differently by
+    Dakota. For example, an analysis package variable may be defined as an 
+    array, so the QmuSetupDesign brancher will create Dakota variables for each 
+    element of the array. Finally, it calls the functions to write the .in and 
+    .m files. this function is independent of the particular analysis package.
 
-  this function defines the data to write the dakota .in and
-  .m files.  it is necessary for multiple reasons.  first,
-  it collects the parameters and applies some defaults that
-  are unique to the environment.  second, some analysis package
-  variables and / or responses may be treated differently by
-  dakota.  for example, an analysis package variable may be
-  defined as an array, so the QmuSetupDesign brancher will
-  create dakota variables for each element of the array.
-  finally it calls the functions to write the .in and .m files.
-  this function is independent of the particular analysis
-  package.
-
-  this data would typically be generated by a matlab script
-  for a specific model, using the method, variable, and
-  response class objects.
-'''
+    This data would typically be generated by a MATLAB script for a specific 
+    model, using the method, variable, and response class objects.
+    '''
 
     #  parameters
Index: /issm/trunk-jpl/src/m/qmu/postqmu.m
===================================================================
--- /issm/trunk-jpl/src/m/qmu/postqmu.m	(revision 25162)
+++ /issm/trunk-jpl/src/m/qmu/postqmu.m	(revision 25163)
@@ -1,5 +1,7 @@
 function md=postqmu(md)
-%INPUT function md=postqmu(md)
-%Deal with dakota output results in files.
+%POSTQMU - Deal with Dakota output results in files
+%
+%   Usage:
+%      md = postqmu(md)
 
 %  check to see if dakota returned errors in the err file
@@ -14,5 +16,5 @@
            fline=fgetl(fide);
        end
-       warning(['Dakota returned error in ''' qmuerrfile ' file.  qmu directory retained.'])
+       warning(['Dakota returned error in ''' qmuerrfile ' file. qmu directory retained.'])
    end
    status=fclose(fide);
@@ -46,5 +48,4 @@
 end
 
-
 %put dakotaresults in their right location.
 md.results.dakota=dakotaresults;
Index: /issm/trunk-jpl/src/m/qmu/postqmu.py
===================================================================
--- /issm/trunk-jpl/src/m/qmu/postqmu.py	(revision 25162)
+++ /issm/trunk-jpl/src/m/qmu/postqmu.py	(revision 25163)
@@ -1,16 +1,22 @@
 from os import getpid, stat
 from os.path import isfile
-from subprocess import Popen
+import shlex
+from subprocess import call
+
 from dakota_out_parse import *
 from helpers import *
+from loadresultsfromdisk import *
 
 
 def postqmu(md):
     '''
-    Deal with dakota output results in files.
+    POSTQMU - Deal with Dakota output results in files
+    
+        Usage:
+            md = postqmu(md)
 
-    INPUT function
-    md = postqmu(md)
-
+        TODO:
+        - Run Dakota test to check that updates from 6/26 are working
+        - Add checks to Popen
     '''
 
@@ -41,4 +47,14 @@
         dakotaresults.dresp_dat = dresp_dat
 
+    if md.qmu.output:
+        if md.qmu.method.method == 'nond_sampling':
+            dakotaresults.modelresults = []
+            md2 = copy.deepcopy(md)
+            md2.qmu.isdakota = 0
+            for i in range(md2.qmu.method.params.samples):
+                print('reading qmu file {}.outbin.{}'.format(md2.miscellaneous.name, i))
+                md2 = loadresultsfromdisk(md2, '{}.outbin.{}'.format(md2.miscellaneous.name, i))
+                dakotaresults.modelresults.append(md2.results)
+
     # put dakotaresults in their right location.
     md.results.dakota = dakotaresults
@@ -46,7 +62,10 @@
     # move all the individual function evalutations into zip files
     if not md.qmu.isdakota:
-        Popen('zip - mq params.in.zip params.in.[1 - 9] * ', shell=True)
-        Popen('zip - mq results.out.zip results.out.[1 - 9] * ', shell=True)
-        Popen('zip - mq matlab.out.zip matlab * .out.[1 - 9] * ', shell=True)
+        subproc_args = shlex.split('zip -mq params.in.zip params.in.[1-9]*')
+        call(subproc_args, shell=True)
+        subproc_args = shlex.split('zip -mq results.out.zip results.out.[1-9]*')
+        call(subproc_args, shell=True)
+        subproc_args = shlex.split('zip -mq matlab.out.zip matlab*.out.[1-9]*')
+        call(subproc_args, shell=True)
 
     return md
Index: /issm/trunk-jpl/src/m/shp/shpwrite.py
===================================================================
--- /issm/trunk-jpl/src/m/shp/shpwrite.py	(revision 25162)
+++ /issm/trunk-jpl/src/m/shp/shpwrite.py	(revision 25163)
@@ -3,4 +3,5 @@
 except ImportError:
     print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled")
+
 
 def shpwrite(shp, filename): #{{{
Index: /issm/trunk-jpl/test/NightlyRun/test2010.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2010.m	(revision 25162)
+++ /issm/trunk-jpl/test/NightlyRun/test2010.m	(revision 25163)
@@ -29,6 +29,8 @@
 mask=gmtmask(md.mesh.lat,md.mesh.long);
 icemask=ones(md.mesh.numberofvertices,1);
-pos=find(mask==0);  icemask(pos)=-1;
-pos=find(sum(mask(md.mesh.elements),2)<3);   icemask(md.mesh.elements(pos,:))=-1;
+pos=find(mask==0);
+icemask(pos)=-1;
+pos=find(sum(mask(md.mesh.elements),2)<3);
+icemask(md.mesh.elements(pos,:))=-1;
 md.mask.ice_levelset=icemask;
 md.mask.ocean_levelset=-icemask;
@@ -65,7 +67,7 @@
 % }}}
 
-
 %eustatic + rigid + elastic run:
-md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1;
+md.solidearth.settings.rigid=1;
+md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1;
 md.cluster=generic('name',oshostname(),'np',3);
 
Index: /issm/trunk-jpl/test/NightlyRun/test2010.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2010.py	(revision 25162)
+++ /issm/trunk-jpl/test/NightlyRun/test2010.py	(revision 25163)
@@ -22,13 +22,13 @@
 longe = sum(md.mesh.long[md.mesh.elements - 1], 1) / 3
 
-md.solidearth.deltathickness = np.zeros((md.mesh.numberofelements, ))
+md.solidearth.surfaceload.icethicknesschange = np.zeros(md.mesh.numberofelements)
 pos = np.intersect1d(np.array(np.where(late < -75)), np.array(np.where(longe > 0)))
 #python does not include last element in array slices, (6:7) -> [5:7]
-md.solidearth.deltathickness[pos[5:7]] = -1
+md.solidearth.surfaceload.icethicknesschange[pos[5:7]] = -1
 
-md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, ))
-md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, ))
-md.dsl.sea_surface_height_change_above_geoid=np.zeros((md.mesh.numberofvertices+1, ))
-md.dsl.sea_water_pressure_change_at_sea_floor=np.zeros((md.mesh.numberofvertices+1, ))
+md.solidearth.sealevel = np.zeros(md.mesh.numberofvertices)
+md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1))
+md.dsl.sea_surface_height_change_above_geoid = np.zeros(md.mesh.numberofvertices + 1)
+md.dsl.sea_water_pressure_change_at_sea_floor = np.zeros(md.mesh.numberofvertices + 1)
 
 md.solidearth.ocean_area_scaling = 1
@@ -40,35 +40,37 @@
 #mask:  {{{
 mask = gmtmask(md.mesh.lat, md.mesh.long)
-icemask = np.ones((md.mesh.numberofvertices, ))
+icemask = np.ones(md.mesh.numberofvertices)
 pos = np.where(mask == 0)
 icemask[pos[0]] = -1
-pos = np.where(sum(mask[md.mesh.elements - 1], 1) < 3)
+pos = np.where(np.sum(mask[md.mesh.elements - 1], 1) < 3)
 icemask[md.mesh.elements[pos, :] - 1] = -1
 md.mask.ice_levelset = icemask
 md.mask.ocean_levelset = -icemask
 
-#make sure that the ice level set is all inclusive:
-md.mask.land_levelset = np.zeros((md.mesh.numberofvertices, ))
-md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, ))
-
-#make sure that the elements that have loads are fully grounded:
-pos = np.nonzero(md.solidearth.deltathickness)[0]
-md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1
+#make sure that the elements that have loads are fully grounded
+pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0]
+md.mask.ocean_levelset = -np.ones(md.mesh.numberofvertices)
 
 #make sure wherever there is an ice load, that the mask is set to ice:
+#pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # Do we need to do this twice?
+md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = 1
+
+#make sure wherever there is an ice load, that the mask is set to ice
+md.mask.ice_levelset = np.ones(md.mesh.numberofvertices)
 icemask[md.mesh.elements[pos, :] - 1] = -1
 md.mask.ice_levelset = icemask
 # }}}
+
 #geometry {{{
 di = md.materials.rho_ice / md.materials.rho_water
-md.geometry.thickness = np.ones((md.mesh.numberofvertices, ))
-md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, ))
+md.geometry.thickness = np.ones(md.mesh.numberofvertices)
+md.geometry.surface = (1 - di) * np.zeros(md.mesh.numberofvertices)
 md.geometry.base = md.geometry.surface - md.geometry.thickness
 md.geometry.bed = md.geometry.base
 # }}}
 #materials {{{
-md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, ))
+md.initialization.temperature = 273.25 * np.ones(md.mesh.numberofvertices)
 md.materials.rheology_B = paterson(md.initialization.temperature)
-md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, ))
+md.materials.rheology_n = 3 * np.ones(md.mesh.numberofelements)
 # }}}
 #Miscellaneous {{{
@@ -76,19 +78,13 @@
 # }}}
 #Solution parameters {{{
-md.solidearth.reltol = float('NaN')
-md.solidearth.abstol = 1e-3
-md.solidearth.geodetic = 1
+md.solidearth.settings.reltol = np.nan
+md.solidearth.settings.abstol = 1e-3
+md.solidearth.settings.geodetic = 1
 # }}}
 
-#New stuff
-md.solidearth.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, ))
-md.solidearth.Ngia = np.zeros((md.mesh.numberofvertices, ))
-md.solidearth.Ugia = np.zeros((md.mesh.numberofvertices, ))
-md.solidearth.hydro_rate = np.zeros((md.mesh.numberofvertices, ))
-
 #eustatic + rigid + elastic run:
-md.solidearth.rigid = 1
-md.solidearth.elastic = 1
-md.solidearth.rotation = 1
+md.solidearth.settings.rigid = 1
+md.solidearth.settings.elastic = 1
+md.solidearth.settings.rotation = 1
 md.cluster = generic('name', gethostname(), 'np', 3)
 
@@ -96,5 +92,5 @@
 md = solve(md, 'Sealevelrise')
 eus = md.results.SealevelriseSolution.SealevelRSLEustatic
-solidearth = md.results.SealevelriseSolution.Sealevel
+slr = md.results.SealevelriseSolution.Sealevel
 moixz = md.results.SealevelriseSolution.SealevelInertiaTensorXZ
 moiyz = md.results.SealevelriseSolution.SealevelInertiaTensorYZ
@@ -111,5 +107,5 @@
 
 #Fields and tolerances to track changes
-field_names = ['eus', 'solidearth', 'moixz', 'moiyz', 'moizz']
+field_names = ['eus', 'slr', 'moixz', 'moiyz', 'moizz']
 field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, 1e-13]
 field_values = [eus, solidearth, moixz, moiyz, moizz]