Index: /issm/trunk-jpl/src/m/plot/plot_manager.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_manager.py	(revision 21439)
+++ /issm/trunk-jpl/src/m/plot/plot_manager.py	(revision 21440)
@@ -24,5 +24,5 @@
 	from plot_overlay import plot_overlay
 
-def plot_manager(md,options,fig,ax):
+def plot_manager(md,options,fig,axgrid,gridindex):
 	'''
 	PLOT_MANAGER - distribute the plots called by plotmodel
@@ -46,4 +46,5 @@
 	options.addfielddefault('ticklabels','on')
 
+	ax=axgrid[gridindex]
 	# {{{ basemap plot TOFIX
 	#if options.exist('basemap'):
@@ -60,5 +61,6 @@
 	if isinstance(data,(str,unicode)):
 		if data=='mesh': 
-			plot_mesh(md,options,fig,ax)
+			plot_mesh(md,options,fig,axgrid,gridindex)
+
 			#fig.delaxes(fig.axes[1]) # hack to remove colorbar after the fact
 			return
Index: /issm/trunk-jpl/src/m/plot/plot_mesh.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_mesh.py	(revision 21439)
+++ /issm/trunk-jpl/src/m/plot/plot_mesh.py	(revision 21440)
@@ -3,9 +3,11 @@
 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
-
-def plot_mesh(md,options,fig,ax):
+from matplotlib.patches import Polygon
+from mpl_toolkits.mplot3d.art3d import Line3DCollection
+def plot_mesh(md,options,fig,axgrid,gridindex):
 	'''
 	PLOT_MESH - plot model mesh
@@ -16,12 +18,40 @@
 		See also: PLOTMODEL
 	'''
+	x,y,z,elements,is2d,isplanet=processmesh(md,'mesh',options)
 
-	x,y,z,elements,is2d,isplanet=processmesh(md,[],options)
-
+	ax=axgrid[gridindex]
+	fig.delaxes(axgrid.cbar_axes[gridindex])
+	
 	if is2d:
 		ax.triplot(x,y,elements)
 	else:
-		print 'WARNING: only 2D mesh plot is currently implemented'
-	
+		fig.delaxes(ax)
+		geom=int(str(axgrid.get_geometry()[0])+str(axgrid.get_geometry()[1])+str(gridindex+1))
+		ax = fig.add_subplot(geom,projection='3d')
+
+		AB=elements[:,0:2]
+		BC=elements[:,1:3]
+		CA=np.vstack((elements[:,2],elements[:,0])).T
+		DE=elements[:,3:5]
+		EF=elements[:,4:]
+		FD=np.vstack((elements[:,5],elements[:,3])).T
+		AD=np.vstack((elements[:,0],elements[:,3])).T
+		BE=np.vstack((elements[:,1],elements[:,4])).T
+		CF=np.vstack((elements[:,2],elements[:,5])).T
+
+		tmpa=np.vstack((AB,BC,CA,DE,EF,FD,AD,BE,CF))
+		#deleting segments that appear multiple times
+		tmpb = np.ascontiguousarray(tmpa).view(np.dtype((np.void, tmpa.dtype.itemsize * tmpa.shape[1])))
+		_, idx = np.unique(tmpb, return_index=True)
+		triel= tmpa[idx]
+		
+		for triangle in triel:
+			tri=zip(x[triangle],y[triangle],z[triangle])
+			pl3=Line3DCollection([tri],edgecolor='r')
+			ax.add_collection3d(pl3)
+
+		ax.set_xlim([min(x),max(x)])
+		ax.set_ylim([min(y),max(y)])
+		ax.set_zlim([min(z),max(z)])
 	#apply options
 	options.addfielddefault('title','Mesh')
Index: /issm/trunk-jpl/src/m/plot/plotmodel.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plotmodel.py	(revision 21439)
+++ /issm/trunk-jpl/src/m/plot/plotmodel.py	(revision 21440)
@@ -2,14 +2,14 @@
 from plotoptions import plotoptions
 from plotdoc import plotdoc
+from plot_manager import plot_manager
+from math import ceil, sqrt
 
 try:
 	import pylab as p
 	import matplotlib.pyplot as plt
-        from mpl_toolkits.axes_grid1 import ImageGrid, AxesGrid
+	from mpl_toolkits.axes_grid1 import ImageGrid, AxesGrid
+	from mpl_toolkits.mplot3d import Axes3D
 except ImportError:
 	print "could not import pylab, matplotlib has not been installed, no plotting capabilities enabled"
-
-from plot_manager import plot_manager
-from math import ceil, sqrt
 
 def plotmodel(md,*args):
@@ -83,5 +83,5 @@
 		cbar_pad=options.list[0].getfieldvalue('colorbarpad','2.5%') # None or %
 
-		axgrid=ImageGrid(fig, 111,
+		axgrid=ImageGrid(fig,111,
 				nrows_ncols=(nrows,ncols),
 				ngrids=plotnum,
@@ -94,6 +94,5 @@
 				cbar_location=cbar_location,
 				cbar_size=cbar_size,
-				cbar_pad=cbar_pad
-				)
+				cbar_pad=cbar_pad)
 
 		if cbar_mode=='None':
@@ -101,7 +100,6 @@
 				fig._axstack.remove(ax)
 
-		for i in xrange(numberofplots):
-			plot_manager(options.list[i].getfieldvalue('model',md),options.list[i],fig,axgrid[i])
-
+		for i,ax in enumerate(axgrid.axes_all):
+			plot_manager(options.list[i].getfieldvalue('model',md),options.list[i],fig,axgrid,i)
 		fig.show()
 	else:
Index: /issm/trunk-jpl/src/m/plot/processmesh.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/processmesh.py	(revision 21439)
+++ /issm/trunk-jpl/src/m/plot/processmesh.py	(revision 21440)
@@ -15,51 +15,33 @@
 	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')
+	if md.mesh.numberofvertices==md.mesh.numberofelements:
+		raise ValueError('processmesh error: the number of elements is the same as the number of nodes')
 	# }}}
-	# {{{ treating non data plots mesh
-	if len(data)==0 or not isinstance(data,dict):
-		if 'latlon' not in options.getfieldvalue('coord','xy').lower(): #convert to lower case for comparison
-			try:
-				x=md.mesh.x2d
-			except AttributeError:
-				x=md.mesh.x
-			try:
-				y=md.mesh.y2d
-			except AttributeError:				
-				y=md.mesh.y
-		else:
-			x=md.mesh.long
-			y=md.mesh.lat
-		try:
-			z=md.mesh.z
-		except AttributeError:
-			z=np.zeros_like(md.mesh.x)
+  # {{{ treating coordinates
+
+	try:
+		z=md.mesh.z
+	except AttributeError:
+		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
+
+	#is it a 2D plot?
+	if md.mesh.dimension()==2 or options.getfieldvalue('layer',0)>=1:
+		is2d=1
+	else:
+		is2d=0
 		
-		try:
-			elements=md.mesh.elements2d-1
-		except AttributeError:
-			elements=md.mesh.elements-1
-
-		#is it a 2D plot?
-		if md.mesh.dimension()==2 or options.getfieldvalue('layer',0)>=1:
-			is2d=1
-		else:
-			is2d=0
-
-		#layer projection?
-		if options.getfieldvalue('layer',0)>=1:
-			 if 'latlon' in options.getfieldvalue('coord','xy').lower():
-				 raise ValueError('processmesh error: cannot work with 3D mesh in lat-lon coords')
-			 #we modify the mesh temporarily to a 2D mesh from which the 3D mesh was extruded
-			 z=np.zeros(np.size(x))
-			 elements=elements
-	else:
-		#Process mesh for plotting 
-		if md.mesh.dimension()==2:
-			is2d=1
-		else:
-			# process polycollection here for 3D plot
-			is2d=0
 	#units
 	if options.exist('unit'):