Index: /issm/trunk-jpl/test/NightlyRun/test2003.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2003.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2003.py	(revision 22134)
@@ -0,0 +1,113 @@
+#Test Name: EarthSlr_rotationalFeedback
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from gmshplanet import *
+from love_numbers import *
+from paterson import *
+from maskpsl import *
+from gmtmask import *
+
+#mesh earth:
+md = model()
+md.mesh = gmshplanet('radius',6.371012*1e3,'resolution',1000) #500 km resolution mesh
+
+#parameterize slr solution:
+#slr loading:  {{{
+md.slr.deltathickness = np.zeros((md.mesh.numberofelements,))
+md.slr.sealevel = np.zeros((md.mesh.numberofvertices,))
+md.slr.steric_rate = np.zeros((md.mesh.numberofvertices,))
+
+#antarctica
+#Access every element in lat using the indices in elements
+#-1 to convert to base 0 indexing, 1 (not 2, in matlab) to sum over rows
+late = sum(md.mesh.lat[md.mesh.elements - 1],1) / 3
+longe = sum(md.mesh.long[md.mesh.elements - 1],1) / 3
+pos = np.intersect1d(np.where(late < -75), np.where(longe < 0))
+md.slr.deltathickness[pos] = -1
+
+#elastic loading from love numbers: 
+nlov = 1000
+md.slr.love_h = np.array(love_numbers('h'))
+md.slr.love_h = np.resize(md.slr.love_h, nlov + 1)
+md.slr.love_k = np.array(love_numbers('k'))
+md.slr.love_k = np.resize(md.slr.love_k, nlov + 1)
+md.slr.love_l = np.array(love_numbers('l'))
+md.slr.love_l = np.resize(md.slr.love_l, nlov + 1)
+#}}}
+
+#mask:  {{{
+md.mask = maskpsl() # use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset 
+mask = gmtmask(md.mesh.lat,md.mesh.long) 
+
+icemask = np.ones((md.mesh.numberofvertices,))
+pos = np.where(mask == 0)
+#pos[0] because np.where(mask==0) returns a 2d array, the latter parts of which are all array/s of 0s
+icemask[pos[0]] = -1
+pos = np.where(sum(mask[md.mesh.elements - 1],1) < 3)
+icemask[md.mesh.elements[pos,:] - 1] = -1
+md.mask.ice_levelset = icemask
+md.mask.ocean_levelset = np.zeros((md.mesh.numberofvertices,))
+pos = np.where(md.mask.ice_levelset == 1)
+md.mask.ocean_levelset[pos] = 1
+
+#make sure that the ice level set is all inclusive:
+md.mask.land_levelset = np.zeros((md.mesh.numberofvertices,))
+md.mask.groundedice_levelset = -np.ones((md.mesh.numberofvertices,))
+
+#make sure wherever there is an ice load, that the mask is set to ice: 
+pos = np.where(md.slr.deltathickness)
+md.mask.ice_levelset[md.mesh.elements[pos,:]] = -1
+# }}}
+
+# use model representation of ocea area (not the ture area) 
+md.slr.ocean_area_scaling = 0
+
+#geometry
+di = md.materials.rho_ice / md.materials.rho_water
+md.geometry.thickness = np.ones((md.mesh.numberofvertices,))
+md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices,))
+md.geometry.base = md.geometry.surface - md.geometry.thickness
+md.geometry.bed = md.geometry.base
+
+#materials
+md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_B = paterson(md.initialization.temperature)
+md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements,))
+
+#Miscellaneous
+md.miscellaneous.name = 'test2003'
+
+#Solution parameters
+md.slr.reltol = float('NaN')
+md.slr.abstol = 1e-3
+
+#eustatic + rigid + elastic run: 
+md.slr.rigid = 1
+md.slr.elastic = 1
+md.slr.rotation = 0
+md.cluster = generic('name',gethostname(),'np',3)
+#md.verbose = verbose('111111111')
+#print md.calving
+#print md.gia
+#print md.love
+#print md.esa
+#print md.autodiff
+md = solve(md,'Sealevelrise')
+SnoRotation = md.results.SealevelriseSolution.Sealevel
+
+#eustatic + rigid + elastic + rotation run: 
+md.slr.rigid = 1
+md.slr.elastic = 1
+md.slr.rotation = 1
+md.cluster = generic('name',gethostname(),'np',3)
+#md.verbose = verbose('111111111')
+md = solve(md,'Sealevelrise')
+SRotation = md.results.SealevelriseSolution.Sealevel
+
+#Fields and tolerances to track changes
+field_names = ['noRotation', 'Rotation']
+field_tolerances = [1e-13, 1e-13]
+field_values = [SnoRotation, SRotation]
+
Index: /issm/trunk-jpl/test/NightlyRun/test2010.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2010.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2010.py	(revision 22134)
@@ -0,0 +1,110 @@
+#Test Name: MomentOfIntertia
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from gmshplanet import *
+from love_numbers import *
+from paterson import *
+from maskpsl import *
+from gmtmask import *
+
+#mesh earth: 
+md = model()
+rad_e = 6.371012*1e3 # mean radius of Earth, km
+md.mesh = gmshplanet('radius',rad_e,'resolution',1000.0)  # km resolution 
+
+#parameterize slr solution:
+#slr loading:  {{{
+late = sum(md.mesh.lat[md.mesh.elements - 1],1) / 3
+longe = sum(md.mesh.long[md.mesh.elements - 1],1) / 3
+
+md.slr.deltathickness = np.zeros((md.mesh.numberofelements,))
+pos = np.intersect1d(np.where(late < -75), np.where(longe > 0))
+#python does not include last element in array slices, (6:7) -> [5:7]
+md.slr.deltathickness[pos[5:7]] = -1
+
+md.slr.sealevel = np.zeros((md.mesh.numberofvertices,))
+md.slr.steric_rate = np.zeros((md.mesh.numberofvertices,))
+md.slr.ocean_area_scaling = 1
+
+#elastic loading from love numbers: 
+nlov = 1000
+md.slr.love_h = np.array(love_numbers('h'))
+md.slr.love_h = np.resize(md.slr.love_h, nlov + 1)
+md.slr.love_k = np.array(love_numbers('k'))
+md.slr.love_k = np.resize(md.slr.love_k, nlov + 1)
+md.slr.love_l = np.array(love_numbers('l'))
+md.slr.love_l = np.resize(md.slr.love_l, nlov + 1)
+
+#}}}
+#mask:  {{{
+md.mask = maskpsl() # use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset 
+mask = gmtmask(md.mesh.lat,md.mesh.long)
+
+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)
+icemask[md.mesh.elements[pos,:] - 1] = -1
+md.mask.ice_levelset = icemask
+md.mask.ocean_levelset = np.zeros((md.mesh.numberofvertices,))
+pos = np.where(md.mask.ice_levelset == 1)
+md.mask.ocean_levelset[pos] = 1
+
+#make sure that the ice level set is all inclusive:
+md.mask.land_levelset = np.zeros((md.mesh.numberofvertices,))
+md.mask.groundedice_levelset = -np.ones((md.mesh.numberofvertices,))
+
+#make sure wherever there is an ice load, that the mask is set to ice: 
+#md.mask.ice_levelset = np.ones((md.mesh.numberofvertices,))
+pos = np.where(md.slr.deltathickness)
+md.mask.ice_levelset[md.mesh.elements[pos,:]] = -1
+# }}}
+#geometry {{{
+di = md.materials.rho_ice / md.materials.rho_water
+md.geometry.thickness = np.ones((md.mesh.numberofvertices,))
+md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices,))
+md.geometry.base = md.geometry.surface - md.geometry.thickness
+md.geometry.bed = md.geometry.base
+# }}}
+#materials {{{
+md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_B = paterson(md.initialization.temperature)
+md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements,))
+# }}}
+#Miscellaneous {{{
+md.miscellaneous.name = 'test2010'
+# }}}
+#Solution parameters {{{
+md.slr.reltol = float('NaN')
+md.slr.abstol = 1e-3
+# }}}
+
+#eustatic + rigid + elastic run: 
+md.slr.rigid = 1
+md.slr.elastic = 1
+md.slr.rotation = 1
+md.cluster = generic('name',gethostname(),'np',3)
+
+# uncomment following 2 lines for 
+md = solve(md,'Sealevelrise')
+eus = md.results.SealevelriseSolution.SealevelEustatic
+slr = md.results.SealevelriseSolution.Sealevel
+moixz = md.results.SealevelriseSolution.SealevelInertiaTensorXZ
+moiyz = md.results.SealevelriseSolution.SealevelInertiaTensorYZ
+moizz = md.results.SealevelriseSolution.SealevelInertiaTensorZZ 
+
+# analytical moi => just checking FOR ICE only!!! {{{
+# ...have to mute ** slr induced MOI in Tria.cpp ** prior to the comparison 
+#rad_e = rad_e*1e3 # now in meters 
+#areas=GetAreasSphericalTria(md.mesh.elements,md.mesh.lat,md.mesh.long,rad_e)
+#lat=late*pi/180 lon=longe*pi/180
+#moi_xz = sum(-md.materials.rho_freshwater.*md.slr.deltathickness.*areas.*rad_e^2.*sin(lat).*cos(lat).*cos(lon))
+#moi_yz = sum(-md.materials.rho_freshwater.*md.slr.deltathickness.*areas.*rad_e^2.*sin(lat).*cos(lat).*sin(lon))
+# }}}
+
+#Fields and tolerances to track changes
+field_names     = ['eus','slr','moixz','moiyz','moizz']
+field_tolerances= [1e-13,1e-13,1e-13,1e-13,1e-13]
+field_values = [eus,slr,moixz,moiyz,moizz]
Index: /issm/trunk-jpl/test/NightlyRun/test2051.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2051.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2051.py	(revision 22134)
@@ -0,0 +1,37 @@
+#Test Name: GiaIvinsBenchmarksAB2dA1
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from triangle import *
+from setmask import *
+from parameterize import *
+
+# Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.) 
+md = triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/GiaIvinsBenchmarksAB.py')
+
+# indicate what you want to compute 
+md.gia.cross_section_shape = 1 # for square-edged x-section 
+
+# define loading history 
+md.timestepping.start_time = 2002100 # after 2 kyr of deglaciation 
+md.timestepping.final_time = 2500000 # 2,500 kyr
+md.geometry.thickness = np.array([
+	np.append(md.geometry.thickness * 0.0 ,0.0),
+	np.append(md.geometry.thickness, 1000),
+	np.append(md.geometry.thickness, 2000000),
+	np.append(md.geometry.thickness * 0.0, 2000100),
+	np.append(md.geometry.thickness * 0.0, md.timestepping.start_time)
+	]).T
+
+# solve for GIA deflection 
+md.cluster=generic('name',gethostname(),'np',3)
+md.verbose = verbose('1111111')
+md = solve(md,'Gia')
+
+#Fields and tolerances to track changes
+field_names      = ['GiaW','GiadWdt']
+field_tolerances = [1e-13,1e-13]
+field_values = [md.results.GiaSolution.GiaW, md.results.GiaSolution.GiadWdt]
Index: /issm/trunk-jpl/test/NightlyRun/test2052.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2052.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2052.py	(revision 22134)
@@ -0,0 +1,37 @@
+#Test Name: GiaIvinsBenchmarksAB2dA2
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from triangle import *
+from setmask import *
+from parameterize import *
+
+#Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.) 
+md = triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/GiaIvinsBenchmarksAB.py')
+
+#indicate what you want to compute 
+md.gia.cross_section_shape = 1    # for square-edged x-section 
+
+#define loading history 
+md.timestepping.start_time = 2005100 # after 5 kyr of deglaciation 
+md.timestepping.final_time = 2500000 # 2,500 kyr
+md.geometry.thickness = np.array([
+	np.append(md.geometry.thickness * 0.0 ,0.0),
+	np.append(md.geometry.thickness, 1000),
+	np.append(md.geometry.thickness, 2000000),
+	np.append(md.geometry.thickness * 0.0, 2000100),
+	np.append(md.geometry.thickness * 0.0, md.timestepping.start_time)
+	]).T
+
+#solve for GIA deflection 
+md.cluster = generic('name',gethostname(),'np',3)
+md.verbose = verbose('1111111')
+md = solve(md,'Gia')
+
+#Fields and tolerances to track changes
+field_names      = ['GiaW','GiadWdt']
+field_tolerances = [1e-13,1e-13]
+field_values = [md.results.GiaSolution.GiaW,md.results.GiaSolution.GiadWdt]
Index: /issm/trunk-jpl/test/NightlyRun/test2053.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2053.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2053.py	(revision 22134)
@@ -0,0 +1,37 @@
+#Test Name: GiaIvinsBenchmarksAB2dA3
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from triangle import *
+from setmask import *
+from parameterize import *
+
+#Benchmark experiments (Figure A2a Ivins and James, 1999, Geophys. J. Int.) 
+md = triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/GiaIvinsBenchmarksAB.py')
+
+#indicate what you want to compute 
+md.gia.cross_section_shape = 1    # for square-edged x-section 
+
+#define loading history 
+md.timestepping.start_time = 2010100 # after 10 kyr of deglaciation 
+md.timestepping.final_time = 2500000 # 2,500 kyr
+md.geometry.thickness = np.array([
+	np.append(md.geometry.thickness * 0.0 ,0.0),
+	np.append(md.geometry.thickness, 1000),
+	np.append(md.geometry.thickness, 2000000),
+	np.append(md.geometry.thickness * 0.0, 2000100),
+	np.append(md.geometry.thickness * 0.0, md.timestepping.start_time)
+	]).T
+
+#solve for GIA deflection 
+md.cluster = generic('name',gethostname(),'np',3)
+md.verbose = verbose('1111111')
+md = solve(md,'Gia')
+
+#Fields and tolerances to track changes
+field_names      = ['GiaW','GiadWdt']
+field_tolerances = [1e-13,1e-13]
+field_values = [md.results.GiaSolution.GiaW,md.results.GiaSolution.GiadWdt]
Index: /issm/trunk-jpl/test/NightlyRun/test2071.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2071.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2071.py	(revision 22134)
@@ -0,0 +1,36 @@
+#Test Name: GiaIvinsBenchmarksAB2dC1
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from triangle import *
+from setmask import *
+from parameterize import *
+
+#Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.) 
+md = triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/GiaIvinsBenchmarksCD.py')
+
+#indicate what you want to compute 
+md.gia.cross_section_shape = 1 # for square-edged x-section 
+
+#define loading history 
+md.timestepping.start_time = 0.3 # for t \approx 0 kyr : to get eleastic response!   
+md.timestepping.final_time = 2500000 # 2,500 kyr
+md.geometry.thickness = np.array([
+	np.append(md.geometry.thickness * 0.0, 0.0),
+	np.append(md.geometry.thickness / 2.0, 0.1),
+	np.append(md.geometry.thickness, 0.2),
+	np.append(md.geometry.thickness, md.timestepping.start_time)
+	]).T
+
+#solve for GIA deflection 
+md.cluster = generic('name',gethostname(),'np',3)
+md.verbose = verbose('1111111')
+md = solve(md,'Gia')
+
+#Fields and tolerances to track changes
+field_names      = ['GiaW','GiadWdt']
+field_tolerances = [1e-13,1e-13]
+field_values = [md.results.GiaSolution.GiaW, md.results.GiaSolution.GiadWdt]
Index: /issm/trunk-jpl/test/NightlyRun/test2072.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2072.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2072.py	(revision 22134)
@@ -0,0 +1,36 @@
+#Test Name: GiaIvinsBenchmarksAB2dC2
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from triangle import *
+from setmask import *
+from parameterize import *
+
+#Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.) 
+md = triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/GiaIvinsBenchmarksCD.py')
+
+#indicate what you want to compute 
+md.gia.cross_section_shape = 1 # for square-edged x-section 
+
+#define loading history 
+md.timestepping.start_time = 1000.3 # for t \approx 1 kyr   
+md.timestepping.final_time = 2500000 # 2,500 kyr
+md.geometry.thickness = np.array([
+	np.append(md.geometry.thickness * 0.0, 0.0),
+	np.append(md.geometry.thickness / 2.0, 0.1),
+	np.append(md.geometry.thickness, 0.2),
+	np.append(md.geometry.thickness, md.timestepping.start_time)
+	]).T
+
+#solve for GIA deflection 
+md.cluster = generic('name',gethostname(),'np',3)
+md.verbose = verbose('1111111')
+md = solve(md,'Gia')
+
+#Fields and tolerances to track changes
+field_names      = ['GiaW','GiadWdt']
+field_tolerances = [1e-13,1e-13]
+field_values = [md.results.GiaSolution.GiaW, md.results.GiaSolution.GiadWdt]
Index: /issm/trunk-jpl/test/NightlyRun/test2073.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2073.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2073.py	(revision 22134)
@@ -0,0 +1,36 @@
+#Test Name: GiaIvinsBenchmarksAB2dC3
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from triangle import *
+from setmask import *
+from parameterize import *
+
+#Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.) 
+md = triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/GiaIvinsBenchmarksCD.py')
+
+#indicate what you want to compute 
+md.gia.cross_section_shape = 1 # for square-edged x-section 
+
+#define loading history 
+md.timestepping.start_time = 2400000 # for t \approx \infty  
+md.timestepping.final_time = 2500000 # 2,500 kyr
+md.geometry.thickness = np.array([
+	np.append(md.geometry.thickness * 0.0, 0.0),
+	np.append(md.geometry.thickness / 2.0, 0.1),
+	np.append(md.geometry.thickness, 0.2),
+	np.append(md.geometry.thickness, md.timestepping.start_time)
+	]).T
+
+#solve for GIA deflection 
+md.cluster = generic('name',gethostname(),'np',3)
+md.verbose = verbose('1111111')
+md = solve(md,'Gia')
+
+#Fields and tolerances to track changes
+field_names      = ['GiaW','GiadWdt']
+field_tolerances = [1e-13,1e-13]
+field_values = [md.results.GiaSolution.GiaW, md.results.GiaSolution.GiadWdt]
Index: /issm/trunk-jpl/test/NightlyRun/test2081.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2081.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2081.py	(revision 22134)
@@ -0,0 +1,36 @@
+#Test Name: GiaIvinsBenchmarksAB2dD1
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from triangle import *
+from setmask import *
+from parameterize import *
+
+#Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.) 
+md = triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/GiaIvinsBenchmarksCD.py')
+
+#indicate what you want to compute 
+md.gia.cross_section_shape = 2 # for square-edged x-section 
+
+#define loading history 
+md.timestepping.start_time = 0.3 # for t \approx 0 kyr : to get eleastic response!  
+md.timestepping.final_time = 2500000 # 2,500 kyr
+md.geometry.thickness = np.array([
+	np.append(md.geometry.thickness * 0.0, 0.0),
+	np.append(md.geometry.thickness / 2.0, 0.1),
+	np.append(md.geometry.thickness, 0.2),
+	np.append(md.geometry.thickness, md.timestepping.start_time)
+	]).T
+
+#solve for GIA deflection 
+md.cluster = generic('name',gethostname(),'np',3)
+md.verbose = verbose('1111111')
+md = solve(md,'Gia')
+
+#Fields and tolerances to track changes
+field_names      = ['GiaW','GiadWdt']
+field_tolerances = [1e-13,1e-13]
+field_values = [md.results.GiaSolution.GiaW, md.results.GiaSolution.GiadWdt]
Index: /issm/trunk-jpl/test/NightlyRun/test2082.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2082.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2082.py	(revision 22134)
@@ -0,0 +1,36 @@
+#Test Name: GiaIvinsBenchmarksAB2dD2
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from triangle import *
+from setmask import *
+from parameterize import *
+
+#Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.) 
+md = triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/GiaIvinsBenchmarksCD.py')
+
+#indicate what you want to compute 
+md.gia.cross_section_shape = 2 # for square-edged x-section 
+
+#define loading history 
+md.timestepping.start_time = 1000.3 # for t \approx 1 kyr  
+md.timestepping.final_time = 2500000 # 2,500 kyr
+md.geometry.thickness = np.array([
+	np.append(md.geometry.thickness * 0.0, 0.0),
+	np.append(md.geometry.thickness / 2.0, 0.1),
+	np.append(md.geometry.thickness, 0.2),
+	np.append(md.geometry.thickness, md.timestepping.start_time)
+	]).T
+
+#solve for GIA deflection 
+md.cluster = generic('name',gethostname(),'np',3)
+md.verbose = verbose('1111111')
+md = solve(md,'Gia')
+
+#Fields and tolerances to track changes
+field_names      = ['GiaW','GiadWdt']
+field_tolerances = [1e-13,1e-13]
+field_values = [md.results.GiaSolution.GiaW, md.results.GiaSolution.GiadWdt]
Index: /issm/trunk-jpl/test/NightlyRun/test2083.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2083.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2083.py	(revision 22134)
@@ -0,0 +1,36 @@
+#Test Name: GiaIvinsBenchmarksAB2dD3
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from triangle import *
+from setmask import *
+from parameterize import *
+
+#Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.) 
+md = triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000)
+md = setmask(md,'','')
+md = parameterize(md,'../Par/GiaIvinsBenchmarksCD.py')
+
+#indicate what you want to compute 
+md.gia.cross_section_shape = 2 # for square-edged x-section 
+
+#define loading history 
+md.timestepping.start_time = 2400000 # for t \approx \infty  
+md.timestepping.final_time = 2500000 # 2,500 kyr
+md.geometry.thickness = np.array([
+	np.append(md.geometry.thickness * 0.0, 0.0),
+	np.append(md.geometry.thickness / 2.0, 0.1),
+	np.append(md.geometry.thickness, 0.2),
+	np.append(md.geometry.thickness, md.timestepping.start_time)
+	]).T
+
+#solve for GIA deflection 
+md.cluster = generic('name',gethostname(),'np',3)
+md.verbose = verbose('1111111')
+md = solve(md,'Gia')
+
+#Fields and tolerances to track changes
+field_names      = ['GiaW','GiadWdt']
+field_tolerances = [1e-13,1e-13]
+field_values = [md.results.GiaSolution.GiaW, md.results.GiaSolution.GiadWdt]
Index: /issm/trunk-jpl/test/NightlyRun/test2101.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2101.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2101.py	(revision 22134)
@@ -0,0 +1,82 @@
+#Test Name: EarthEsa
+#Elastostatic adjustment for an elemental ice unloading
+
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from gmshplanet import *
+from love_numbers import *
+from paterson import *
+from maskpsl import *
+from gmtmask import *
+
+#mesh earth: 
+md = model()
+md.mesh = gmshplanet('radius',6.371012*1e3,'resolution',1000)
+
+#define load 
+md.esa.deltathickness = np.zeros((md.mesh.numberofelements,))
+pos = 450
+md.esa.deltathickness[pos] = -100   # this is the only "icy" element
+
+#love numbers:
+nlov = 10000
+md.esa.love_h = np.array(love_numbers('h'))	#Originally had CM arg
+md.esa.love_h = np.resize(md.esa.love_h, nlov + 1)
+md.esa.love_l = np.array(love_numbers('l'))	#Originally had CM arg
+md.esa.love_l = np.resize(md.esa.love_l, nlov + 1)
+
+#mask:  {{{
+md.mask = maskpsl() # use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset 
+md.mask.ocean_levelset = gmtmask(md.mesh.lat,md.mesh.long) 
+
+#make sure wherever there is an ice load, that the mask is set to ice: 
+md.mask.ice_levelset = np.ones((md.mesh.numberofvertices,))
+pos = np.where(md.esa.deltathickness)
+md.mask.ice_levelset[md.mesh.elements[pos,:]] = -1
+
+#is ice grounded? 
+md.mask.groundedice_levelset = -np.ones((md.mesh.numberofvertices,))
+pos = np.where(md.mask.ice_levelset <= 0)
+md.mask.groundedice_levelset[pos] = 1
+
+#make sure ice domain is on the continent: 
+#pos = np.where(md.mask.ice_levelset <= 0)
+#md.mask.ocean_levelset[pos] = 0
+
+#land mask 
+md.mask.land_levelset = 1 - md.mask.ocean_levelset
+
+# }}}
+#geometry:  {{{
+di = md.materials.rho_ice / md.materials.rho_water
+md.geometry.thickness = np.ones((md.mesh.numberofvertices,))
+md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices,))
+md.geometry.base = md.geometry.surface - md.geometry.thickness
+md.geometry.bed = md.geometry.base
+# }}}
+#materials:  {{{
+md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_B = paterson(md.initialization.temperature)
+md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements,))
+# }}}
+#Miscellaneous: {{{
+md.miscellaneous.name = 'test2101'
+# }}}
+
+#solve esa 
+md.esa.requested_outputs = ['EsaUmotion','EsaNmotion','EsaEmotion']
+md.cluster = generic('name',gethostname(),'np',3)
+md.verbose = verbose('111111111')
+md = solve(md,'Esa')
+
+#Fields and tolerances to track changes
+field_names     = ['EsaUmotion','EsaNmotion','EsaEmotion']
+field_tolerances= [1e-13,1e-13,1e-13]
+field_values = [
+	md.results.EsaSolution.EsaUmotion,
+	md.results.EsaSolution.EsaNmotion,
+	md.results.EsaSolution.EsaEmotion,
+	]
+
Index: /issm/trunk-jpl/test/NightlyRun/test2110.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2110.py	(revision 22134)
+++ /issm/trunk-jpl/test/NightlyRun/test2110.py	(revision 22134)
@@ -0,0 +1,85 @@
+#Test Name: Esa2Dsurface
+#Elastostatic adjustment for an elemental ice unloading
+
+import numpy as np
+from model import *
+from socket import gethostname
+from solve import *
+from roundmesh import *
+#from gmshplanet import *
+from love_numbers import *
+from paterson import *
+#from maskpsl import *
+#from gmtmask import *
+
+#mesh earth: 
+md = model()
+md = roundmesh(md,50000,2000)     # radius and element size (meters)
+
+#define load 
+md.esa.deltathickness = np.zeros((md.mesh.numberofelements,))
+disc_radius = 20 # km
+index = md.mesh.elements
+x_element = np.mean(md.mesh.x[index - 1], 1)
+y_element = np.mean(md.mesh.y[index - 1], 1)
+rad_dist = np.sqrt(x_element**2 + y_element**2) / 1000  # radial distance in km
+md.esa.deltathickness[np.where(rad_dist <= disc_radius)] = -1   # 1 m water withdrawl
+
+#love numbers:
+nlov = 10000	# horizontal displacements do not work for low degree truncation, e.g., 101
+md.esa.love_h = np.array(love_numbers('h'))
+md.esa.love_h = np.resize(md.esa.love_h, nlov + 1)
+md.esa.love_l = np.array(love_numbers('l'))
+md.esa.love_l = np.resize(md.esa.love_l, nlov + 1)
+
+#mask:  {{{
+#make sure wherever there is an ice load, that the mask is set to ice: 
+md.mask.ice_levelset = np.ones((md.mesh.numberofvertices,))
+pos = np.where(md.esa.deltathickness)
+md.mask.ice_levelset[md.mesh.elements[pos,:]] = -1
+
+#is ice grounded? 
+md.mask.groundedice_levelset = -np.ones((md.mesh.numberofvertices,))
+pos = np.where(md.mask.ice_levelset <= 0)
+md.mask.groundedice_levelset[pos] = 1
+
+# }}}
+#geometry:  {{{
+di = md.materials.rho_ice / md.materials.rho_water
+md.geometry.thickness = np.ones((md.mesh.numberofvertices,))
+md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices,))
+md.geometry.base = md.geometry.surface - md.geometry.thickness
+md.geometry.bed = md.geometry.base
+# }}}
+#materials:  {{{
+md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_B = paterson(md.initialization.temperature)
+md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements,))
+# }}}
+#Miscellaneous: {{{
+md.miscellaneous.name = 'test2110'
+# }}}
+
+md.esa.degacc = 0.01
+
+#solve esa 
+md.esa.requested_outputs = ['EsaUmotion','EsaNmotion','EsaEmotion',
+	'EsaStrainratexx','EsaStrainratexy','EsaStrainrateyy','EsaRotationrate']
+md.cluster = generic('name',gethostname(),'np',3)
+md.verbose = verbose('111111111')
+md = solve(md,'Esa')
+
+#Fields and tolerances to track changes
+field_names     = ['EsaUmotion','EsaNmotion','EsaEmotion',
+	'EsaStrainratexx','EsaStrainratexy','EsaStrainrateyy','EsaRotationrate']
+field_tolerances = [1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]
+field_values = [
+	md.results.EsaSolution.EsaUmotion,
+	md.results.EsaSolution.EsaNmotion,
+	md.results.EsaSolution.EsaEmotion,
+	md.results.EsaSolution.EsaStrainratexx,
+	md.results.EsaSolution.EsaStrainratexy,
+	md.results.EsaSolution.EsaStrainrateyy,
+	md.results.EsaSolution.EsaRotationrate,
+	]
+
Index: /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.py
===================================================================
--- /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.py	(revision 22134)
+++ /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.py	(revision 22134)
@@ -0,0 +1,84 @@
+#Geometry specific to Experiments A and B
+
+import os.path
+import inspect
+import numpy as np
+from arch import *
+from InterpFromMeshToMesh2d import *
+from paterson import *
+from verbose import *
+from SetIceSheetBC import *
+
+rad = 600000.
+nv = md.mesh.numberofvertices
+if (np.isnan(md.geometry.thickness)):
+	md.geometry.thickness = np.zeros((md.mesh.numberofvertices,))
+for i in range(nv):
+	dist = np.sqrt(md.mesh.x[i]**2 + md.mesh.y[i]**2)
+	if (dist <= rad):
+		md.geometry.thickness[i] = 2000.0
+	else:
+		md.geometry.thickness[i] = 1.0 # non-zero thickness
+
+md.geometry.thickness = md.geometry.thickness.reshape(-1,1)
+md.geometry.base = np.zeros((md.mesh.numberofvertices,))
+md.geometry.surface = md.geometry.thickness + md.geometry.base.reshape(-1,1) #would otherwise create a 91x91 matrix
+
+#Ice density used for benchmarking, not 917 kg/m^3
+md.materials.rho_ice = 1000 #kg m^3
+
+#GIA parameters specific to Experiments A  and B
+md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices,))		#in Pa.s
+md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices,))	#in km
+md.materials.lithosphere_shear_modulus = 6.7*1e10				#in Pa
+md.materials.lithosphere_density = 3.36						#in g/cm^3
+md.materials.mantle_shear_modulus = 1.45*1e11					#in Pa
+md.materials.mantle_density = 3.38						#in g/cm^3
+
+#Initial velocity 
+x     = archread('../Data/SquareSheetConstrained.arch','x')
+y     = archread('../Data/SquareSheetConstrained.arch','y')
+vx    = archread('../Data/SquareSheetConstrained.arch','vx')
+vy    = archread('../Data/SquareSheetConstrained.arch','vy')
+index = archread('../Data/SquareSheetConstrained.arch','index').astype(int)
+
+md.initialization.vx = np.array(InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)).reshape(-1,1)
+md.initialization.vy = np.array(InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y)).reshape(-1,1)
+vx    = None
+vy    = None
+x     = None
+y     = None
+index = None
+md.initialization.vz = np.zeros((md.mesh.numberofvertices,))
+md.initialization.pressure = np.zeros((md.mesh.numberofvertices,))
+
+#Materials
+md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_B = paterson(md.initialization.temperature)
+md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements,))
+
+#Friction
+md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices,))
+md.friction.coefficient[np.where(md.mask.groundedice_levelset < 0.)] = 0.
+md.friction.p = np.ones((md.mesh.numberofelements,))
+md.friction.q = np.ones((md.mesh.numberofelements,))
+
+#Numerical parameters
+md.stressbalance.viscosity_overshoot = 0.0
+md.masstransport.stabilization = 1
+md.thermal.stabilization = 1.
+md.verbose = verbose(0)
+md.settings.waitonlock = 30.
+md.stressbalance.restol = 0.05
+md.steadystate.reltol = 0.05
+md.stressbalance.reltol = 0.05
+md.stressbalance.abstol = float('NaN')
+md.timestepping.time_step = 1.
+md.timestepping.final_time = 3.
+
+#Boundary conditions:
+md = SetIceSheetBC(md)
+
+#Change name so that no test have the same name
+if len(inspect.stack()) > 2:
+	md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]
Index: /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksCD.py
===================================================================
--- /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksCD.py	(revision 22134)
+++ /issm/trunk-jpl/test/Par/GiaIvinsBenchmarksCD.py	(revision 22134)
@@ -0,0 +1,84 @@
+#Geometry specific to Experiments C and D
+
+import os.path
+import inspect
+import numpy as np
+from arch import *
+from InterpFromMeshToMesh2d import *
+from paterson import *
+from verbose import *
+from SetIceSheetBC import *
+
+rad = 800000
+nv = md.mesh.numberofvertices
+if (np.isnan(md.geometry.thickness)):
+	md.geometry.thickness = np.zeros((md.mesh.numberofvertices,))
+for i in range(nv):
+	dist = np.sqrt(md.mesh.x[i]**2 + md.mesh.y[i]**2)
+	if (dist <= rad):
+		md.geometry.thickness[i] = 3000.0
+	else:
+		md.geometry.thickness[i] = 1.0 # non-zero thickness
+
+md.geometry.thickness = md.geometry.thickness.reshape(-1,1)
+md.geometry.base = np.zeros((md.mesh.numberofvertices,))
+md.geometry.surface = md.geometry.thickness + md.geometry.base.reshape(-1,1) #would otherwise create a 91x91 matrix
+
+#Ice density used for benchmarking, not 917 kg/m^3
+md.materials.rho_ice = 1000 #kg m^3
+
+#GIA parameters specific to Experiments A  and B
+md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices,))		#in Pa.s
+md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices,))	#in km
+md.materials.lithosphere_shear_modulus = 6.7*1e10				#in Pa
+md.materials.lithosphere_density = 3.32						#in g/cm^3
+md.materials.mantle_shear_modulus = 1.45*1e11					#in Pa
+md.materials.mantle_density = 3.34						#in g/cm^3
+
+#Initial velocity 
+x     = archread('../Data/SquareSheetConstrained.arch','x')
+y     = archread('../Data/SquareSheetConstrained.arch','y')
+vx    = archread('../Data/SquareSheetConstrained.arch','vx')
+vy    = archread('../Data/SquareSheetConstrained.arch','vy')
+index = archread('../Data/SquareSheetConstrained.arch','index').astype(int)
+
+md.initialization.vx = np.array(InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)).reshape(-1,1)
+md.initialization.vy = np.array(InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y)).reshape(-1,1)
+vx    = None
+vy    = None
+x     = None
+y     = None
+index = None
+md.initialization.vz = np.zeros((md.mesh.numberofvertices,))
+md.initialization.pressure = np.zeros((md.mesh.numberofvertices,))
+
+#Materials
+md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_B = paterson(md.initialization.temperature)
+md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements,))
+
+#Friction
+md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices,))
+md.friction.coefficient[np.where(md.mask.groundedice_levelset < 0.)] = 0.
+md.friction.p = np.ones((md.mesh.numberofelements,))
+md.friction.q = np.ones((md.mesh.numberofelements,))
+
+#Numerical parameters
+md.stressbalance.viscosity_overshoot = 0.0
+md.masstransport.stabilization = 1.
+md.thermal.stabilization = 1.
+md.verbose = verbose(0)
+md.settings.waitonlock = 30.
+md.stressbalance.restol = 0.05
+md.steadystate.reltol = 0.05
+md.stressbalance.reltol = 0.05
+md.stressbalance.abstol = float('NaN')
+md.timestepping.time_step = 1.
+md.timestepping.final_time = 3.
+
+#Boundary conditions:
+md = SetIceSheetBC(md)
+
+#Change name so that no test have the same name
+if len(inspect.stack()) > 2:
+	md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]