Index: /issm/trunk-jpl/src/m/classes/geometry.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/geometry.py	(revision 23830)
+++ /issm/trunk-jpl/src/m/classes/geometry.py	(revision 23831)
@@ -57,4 +57,10 @@
 			if solution=='TransientSolution' and md.transient.isgroundingline:
 				md = checkfield(md,'fieldname','geometry.bed','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
+				if np.any(self.bed > self.base):
+				    md.checkmessage('base<bed on one or more vertex')
+				pos = np.where(md.mask.groundedice_levelset > 0)
+				if np.any(np.abs(self.bed[pos]-self.base[pos]>10**-9):
+				    md.checkmessage('equality base=bed on grounded ice violated')
+				md = checkfield(md,'fieldname','geometry.bed','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
 
 		return md
Index: /issm/trunk-jpl/test/NightlyRun/runme.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/runme.py	(revision 23830)
+++ /issm/trunk-jpl/test/NightlyRun/runme.py	(revision 23831)
@@ -252,5 +252,5 @@
 		args = parser.parse_args()
 
-		md = runme([args.id, args.include_name], [args.exclude, args.exclude_name], args.benchmark, args.procedure, args.output, args.rank, args.numprocs
+		md = runme([args.id, args.include_name], [args.exclude, args.exclude_name], args.benchmark, args.procedure, args.output, args.rank, args.numprocs)
 
 		exit(md)
Index: /issm/trunk-jpl/test/Par/Pig.py
===================================================================
--- /issm/trunk-jpl/test/Par/Pig.py	(revision 23830)
+++ /issm/trunk-jpl/test/Par/Pig.py	(revision 23831)
@@ -18,4 +18,5 @@
 surface   = numpy.array(archread('../Data/Pig.arch','surface'))
 thickness = numpy.array(archread('../Data/Pig.arch','thickness'))
+bed       = numpy.array(archread('../Data/Pig.arch','bed'))
 
 md.inversion.vx_obs   =InterpFromMeshToMesh2d(index,x,y,vx_obs,md.mesh.x,md.mesh.y)[0][:,0]
@@ -24,4 +25,7 @@
 md.geometry.thickness=InterpFromMeshToMesh2d(index,x,y,thickness,md.mesh.x,md.mesh.y)[0][:,0]
 md.geometry.base=md.geometry.surface-md.geometry.thickness
+md.geometry.bed =md.geometry.base
+pos = np.where(md.mask.groundedice_levelset<0)
+md.geometry.bed[pos] =InterpFromMeshToMesh2d(index,x,y,bed,md.mesh.x[pos],md.mesh.y[pos])
 md.initialization.vx=md.inversion.vx_obs
 md.initialization.vy=md.inversion.vy_obs