Index: /issm/trunk-jpl/src/m/plot/plot_BC.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_BC.py	(revision 21587)
+++ /issm/trunk-jpl/src/m/plot/plot_BC.py	(revision 21588)
@@ -8,6 +8,7 @@
 from applyoptions import applyoptions
 from plot_icefront import plot_icefront
+from mpl_toolkits.mplot3d import Axes3D
 
-def plot_BC(md,options,fig,ax):
+def plot_BC(md,options,fig,axgrid,gridindex):
 	'''
 	PLOT_BC - plot model boundary conditions
@@ -18,8 +19,15 @@
 		See also: PLOTMODEL
 	'''
+	x,y,z,elements,is2d,isplanet=processmesh(md,[],options)
 	
+	ax=axgrid[gridindex]
+	fig.delaxes(axgrid.cbar_axes[gridindex])
+
+	if not is2d:
+		ax=inset_locator.inset_axes(axgrid[gridindex],width='100%',height='100%',loc=3,borderpad=0,axes_class=Axes3D)
+
 	#plot neuman
 	plot_icefront(md,options,fig,ax)
-	x,y,z,elements,is2d,isplanet=processmesh(md,[],options)
+
 	XLims=[np.min(x),np.max(x)]
 	YLims=[np.min(y),np.max(y)]
@@ -36,11 +44,18 @@
 							 y[np.where(~np.isnan(md.stressbalance.spcvz))],
 							 marker='o',c='y',s=80,label='vz Dirichlet',linewidth=0)
-
+		ax.scatter(x[np.where(~np.isnan(md.hydrology.spcepl_head))],
+							 y[np.where(~np.isnan(md.hydrology.spcepl_head))],
+							 marker='v',c='r',s=240,label='EPL Head',linewidth=0)
+		ax.scatter(x[np.where(~np.isnan(md.hydrology.spcsediment_head))],
+							 y[np.where(~np.isnan(md.hydrology.spcsediment_head))],
+							 marker='^',c='b',s=240,label='IDS head',linewidth=0)
+		
 		ax.set_xlim(XLims)
 		ax.set_ylim(YLims)
-	ax.legend()
+	ax.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
+						ncol=3, mode="expand", borderaxespad=0.)
 	#apply options
 	options.addfielddefault('title','Boundary conditions')
 	options.addfielddefault('colorbar','off')
-	applyoptions(md,[],options,fig,ax)
+	applyoptions(md,[],options,fig,axgrid,gridindex)
 	
Index: /issm/trunk-jpl/src/m/plot/plot_icefront.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_icefront.py	(revision 21587)
+++ /issm/trunk-jpl/src/m/plot/plot_icefront.py	(revision 21588)
@@ -16,18 +16,24 @@
 #process mesh and data
 	x,y,z,elements,is2d,isplanet=processmesh(md,[],options)
-	icefront=np.where(np.logical_and(np.sum(md.mask.ice_levelset[elements],1)<3,np.sum(md.mask.ice_levelset[elements],1)>-3)) 
+
+	#icefront check
+	icefront=np.where(np.abs(np.sum(md.mask.ice_levelset[elements],1))!=3) 
 	onlyice=np.where(np.sum(md.mask.ice_levelset[elements],1)==-3)
 	noice=np.where(np.sum(md.mask.ice_levelset[elements],1)==3)
 
+	#hydro neumann
+	hydro_neumann=np.where(md.hydrology.neumannflux!=0)
 	#plot mesh
 	ax.triplot(x,y,elements)
 
 	#highlight elements on neumann
-	colors=np.asarray([0.5 for element in elements[icefront]])
-	
-	ax.tripcolor(x,y,elements[icefront],facecolors=colors,alpha=0.5,label='elements on ice front')#,facecolor='b'
+	if len(icefront[0])>0:
+		colors=np.asarray([0.5 for element in elements[icefront]])
+		ax.tripcolor(x,y,elements[icefront],facecolors=colors,alpha=0.5,label='elements on ice front')
+	if len(hydro_neumann[0])>0:
+		colors=np.asarray([0.5 for element in elements[hydro_neumann]])
+		ax.tripcolor(x,y,elements[hydro_neumann],facecolors=colors,alpha=0.5,label='non zero neumann flux for the hydrology')
 
 	#apply options
 	options.addfielddefault('title','Neumann boundary conditions')
 	options.addfielddefault('colorbar','off')
-	applyoptions(md,[],options,fig,ax)
Index: /issm/trunk-jpl/src/m/plot/plot_mesh.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_mesh.py	(revision 21587)
+++ /issm/trunk-jpl/src/m/plot/plot_mesh.py	(revision 21588)
@@ -3,9 +3,8 @@
 except ImportError:
 	print "could not import pylab, matplotlib has not been installed, no plotting capabilities enabled"
+
 import numpy as np
-import matplotlib as mpl
 from processmesh import processmesh
 from applyoptions import applyoptions
-from matplotlib.patches import Polygon
 from mpl_toolkits.mplot3d.art3d import Line3DCollection
 from mpl_toolkits.axes_grid1 import inset_locator
@@ -29,5 +28,4 @@
 	else:
 		ax=inset_locator.inset_axes(axgrid[gridindex],width='100%',height='100%',loc=3,borderpad=0,axes_class=Axes3D)
-		
 		
 		AB=elements[:,0:2]