Index: /issm/trunk-jpl/src/m/plot/applyoptions.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/applyoptions.py	(revision 17755)
+++ /issm/trunk-jpl/src/m/plot/applyoptions.py	(revision 17756)
@@ -1,4 +1,5 @@
 import numpy as npy
 from cmaptools import truncate_colormap
+from plot_contour import plot_contour
 
 try:
@@ -23,4 +24,8 @@
 	#if not isnan(md.mesh.hemisphere):
 	#	options.addfielddefault('hemisphere',md.mesh.hemisphere)
+
+	# get handle to current figure and axes instance
+	fig = p.gcf()
+	ax=p.gca()
 
 	#font {{{
@@ -113,10 +118,9 @@
 
 	#ticklabel notation {{{
-	p.gca().ticklabel_format(style='sci',scilimits=(0,0))
+	#ax.ticklabel_format(style='sci',scilimits=(0,0))
 	#}}}
 
 	#ticklabelfontsize {{{
 	if options.exist('ticklabelfontsize'):
-		ax=p.gca()
 		for label in ax.get_xticklabels() + ax.get_yticklabels():
 			label.set_fontsize(options.getfieldvalue('ticklabelfontsize'))
@@ -131,7 +135,7 @@
 	if options.exist('axis'):
 		if options.getfieldvalue('axis',True)=='off':
-			p.gca().ticklabel_format(style='plain')
-			p.setp(p.gca().get_xticklabels(), visible=False)
-			p.setp(p.gca().get_yticklabels(), visible=False)
+			ax.ticklabel_format(style='plain')
+			p.setp(ax.get_xticklabels(), visible=False)
+			p.setp(ax.get_yticklabels(), visible=False)
 	# }}}
 
@@ -177,4 +181,16 @@
 	defaultmap=truncate_colormap(mpl.cm.gnuplot2,0.1,0.9,128)
 	cmap=options.getfieldvalue('colormap',defaultmap)
+	norm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1])
+	options.addfield('colornorm',norm)
+	cbar_extend=0
+	if options.exist('cmap_set_over'):
+		over=options.getfieldvalue('cmap_set_over','0.5')
+		cmap.set_over(over)
+		cbar_extend+=1
+	if options.exist('cmap_set_under'):
+		under=options.getfieldvalue('cmap_set_under','0.5')
+		cmap.set_under(under)
+		cbar_extend+=2
+	options.addfield('colormap',cmap)
 	#}}}
 
@@ -183,19 +199,7 @@
 	#colorbar {{{
 	if options.getfieldvalue('colorbar',1)==1:
-		fig = p.gcf()
-		ax = p.gca()
 		divider = make_axes_locatable(ax)
 		cax = divider.new_horizontal("5%", pad=0.05, axes_class=mpl.axes.Axes)
 		fig.add_axes(cax) 
-		norm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1])
-		cbar_extend=0
-		if options.exist('cmap_set_over'):
-			over=options.getfieldvalue('cmap_set_over','0.5')
-			cmap.set_over(over)
-			cbar_extend+=1
-		if options.exist('cmap_set_under'):
-			under=options.getfieldvalue('cmap_set_under','0.5')
-			cmap.set_under(under)
-			cbar_extend+=2
 		if cbar_extend==0:
 			extend='neither'
@@ -226,5 +230,8 @@
 	#streamlines
 
-	#contours
+	#contours {{{
+	if options.exist('contourlevels'):
+		plot_contour(md,data,options)
+	#}}}
 
 	#axis positions
Index: /issm/trunk-jpl/src/m/plot/plot_contour.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_contour.py	(revision 17756)
+++ /issm/trunk-jpl/src/m/plot/plot_contour.py	(revision 17756)
@@ -0,0 +1,39 @@
+from averaging import averaging
+import matplotlib.pyplot as plt
+from processmesh import processmesh
+from processdata import processdata
+
+def plot_contour(md,datain,options):
+	'''
+	plot contours of a given field (called within plotmodel)
+
+	Usage:
+		plot_contour(md,data,options)
+
+	See also: plotmodel
+	'''
+
+	x,y,z,elements,is2d,isplanet=processmesh(md,datain,options)
+	data,datatype=processdata(md,datain,options)
+	ax=plt.gca()
+
+	# process data: must be on nodes
+	if datatype==1: # element data
+		data=averaging(md,data,0)
+	elif datatype==2:
+		pass
+	elif datatype==3: # quiver (vector) data
+		data=npy.sqrt(datain**2)
+	else:
+		raise ValueError('datatype not supported in call to plot_contour')
+
+	# contouronly TODO (cla will also clear an overlay image)
+
+	# retrieve necessary options
+	levels=options.getfieldvalue('contourlevels')
+	colors=options.getfieldvalue('contourcolors')
+	norm=options.getfieldvalue('colornorm')
+	linestyles=options.getfieldvalue('contourlinestyles')
+	linewidths=options.getfieldvalue('contourlinewidths')
+
+	ax.tricontour(x,y,elements,data,levels,colors=colors,norm=norm,linestyles=linestyles,linewidths=linewidths)
Index: /issm/trunk-jpl/src/m/plot/plot_manager.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_manager.py	(revision 17755)
+++ /issm/trunk-jpl/src/m/plot/plot_manager.py	(revision 17756)
@@ -30,4 +30,7 @@
 	#initialize plot handle variable
 	handle=None
+
+	# initialize subplot
+	p.subplot(nlines,ncols,i,aspect='equal')
 
 	##basemap plot
@@ -70,6 +73,6 @@
 
 	#standard plot
-	if not handle:
-		p.subplot(nlines,ncols,i,aspect='equal')
+	#if not handle:
+	#	p.subplot(nlines,ncols,i,aspect='equal')
 
 	#plot unit
Index: /issm/trunk-jpl/src/m/plot/plot_overlay.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_overlay.py	(revision 17755)
+++ /issm/trunk-jpl/src/m/plot/plot_overlay.py	(revision 17756)
@@ -52,5 +52,5 @@
 	gtif=gdal.Open(inputname)
 	arr=gtif.ReadAsArray()
-	os.system('rm -rf ./temp.tif')
+	#os.system('rm -rf ./temp.tif')
 	
 	if gtif.RasterCount>=3:  # RGB array
@@ -62,4 +62,18 @@
 	# normalize array
 	arr=arr/npy.float(npy.max(arr.ravel()))
+
+	if options.getfieldvalue('overlayhist',0)==1:
+		ax=plt.gca()
+		num=2
+		while True:
+			if not plt.fignum_exists(num):
+				break
+			else:
+				num+=1
+		plt.figure(num)
+		plt.hist(arr.flatten(),bins=256,range=(0.,1.))
+		plt.title('histogram of overlay image, use for setting overlaylims')
+		plt.sca(ax) # return to original axes/figure
+		
 
 	# get parameters from cropped geotiff
@@ -89,5 +103,6 @@
 
 	norm=mpl.colors.Normalize(vmin=overlaylims[0],vmax=overlaylims[1])
+
 	handle.pcolormesh(xgmap, ygmap, npy.flipud(arr), cmap=mpl.cm.Greys, norm=norm)
-	plt.axes().set_aspect('equal','box')
+	handle.set_aspect('equal','box')
 	return handle
Index: /issm/trunk-jpl/src/m/plot/plotmodel.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plotmodel.py	(revision 17755)
+++ /issm/trunk-jpl/src/m/plot/plotmodel.py	(revision 17756)
@@ -56,6 +56,7 @@
 
 		if not hold: # TODO need to also check whether figurenumber is a new plot so that old plots are not mistakenly cleared
-			p.clf()
+			p.cla()
 
+		#TODO fig, axarray = plt.subplots(nrows,ncols), then pass fix and axarr to plot_manager
 		#if figsize specified
 		if options.list[0].exist('figsize'):