Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 17765)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 17766)
@@ -32,9 +32,8 @@
 	else:
 		#Guess where the ice front is
-		vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices,))
-		pos=numpy.nonzero(numpy.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0,axis=1) > 0)[0]
-		#vertexonfloatingice[md.mesh.elements[pos].astype(int)-1]=1.
-		vertexonfloatingice[md.mesh.elements[pos]-1]=1.
-		vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice>0.)
+		vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices,1))
+		pos=numpy.nonzero(numpy.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0.,axis=1) >0.)[0]
+		vertexonfloatingice[md.mesh.elements[pos].astype(int)-1]=1.
+		vertexonicefront=numpy.logical_and(numpy.reshape(md.mesh.vertexonboundary,(-1,1)),vertexonfloatingice>0.)
 
 #	pos=find(md.mesh.vertexonboundary & ~vertexonicefront);
@@ -43,7 +42,7 @@
 		print "SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually."
 
-	md.stressbalance.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices,))
-	md.stressbalance.spcvy=float('nan')*numpy.ones((md.mesh.numberofvertices,))
-	md.stressbalance.spcvz=float('nan')*numpy.ones((md.mesh.numberofvertices,))
+	md.stressbalance.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	md.stressbalance.spcvy=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	md.stressbalance.spcvz=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
 	md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
 	md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
@@ -69,7 +68,7 @@
 	else:
 		pos=numpy.nonzero(md.mesh.vertexonboundary)[0]
-	md.stressbalance.spcvx[pos]=0
-	md.stressbalance.spcvy[pos]=0
-	md.stressbalance.spcvz[pos]=0
+	md.stressbalance.spcvx[pos[:]]=0
+	md.stressbalance.spcvy[pos[:]]=0
+	md.stressbalance.spcvz[pos[:]]=0
 
 	#Dirichlet Values
@@ -91,21 +90,22 @@
 	#Deal with other boundary conditions
 	if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
-		md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,))
+		md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
 		print "      no balancethickness.thickening_rate specified: values set as zero"
 
-	md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,))
-	md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,))
-	md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,))
+	md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
 
 	if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
-		md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,))
+		md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
 		if hasattr(md.mesh,'vertexonsurface'):
 			pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
 			md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
 		if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:
-			md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,))
-			md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)[0]]=50.*10.**-3    #50mW/m2
+			md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1))
+			md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)]=50.*10.**-3    #50mW/m2
 	else:
 		print "      no thermal boundary conditions created: no observed temperature found"
 
 	return md
+