Index: /issm/trunk-jpl/src/m/plot/applyoptions.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/applyoptions.py	(revision 27086)
+++ /issm/trunk-jpl/src/m/plot/applyoptions.py	(revision 27087)
@@ -11,4 +11,5 @@
 from plot_contour import plot_contour
 from plot_streamlines import plot_streamlines
+from plot_edgeoverlay import plot_edgeoverlay
 
 
@@ -156,12 +157,9 @@
     # {{{ ShowBasins TODO
     # }}}
-    # {{{ clim
-    if options.exist('clim'):
-        lims = options.getfieldvalue('clim')
-        assert len(lims) == 2, 'error, clim should be passed as a list of length 2'
-    elif options.exist('caxis'):
+    # {{{ caxis
+    if options.exist('caxis'):
         lims = options.getfieldvalue('caxis')
         assert len(lims) == 2, 'error, caxis should be passed as a list of length 2'
-        options.addfielddefault('clim', lims)
+        options.addfielddefault('caxis', lims)
     else:
         if len(data) > 0:
@@ -192,4 +190,9 @@
     if options.exist('contourlevels'):
         plot_contour(md, data, options, ax)
+    # }}}
+    # {{{ edgeoverlay
+    if options.exist('edgeoverlay'):
+        edgedata = options.getfieldvalue('edgeoverlay')
+        plot_edgeoverlay(md, edgedata, options, ax)
     # }}}
     # {{{ wrapping TODO
Index: /issm/trunk-jpl/src/m/plot/plot_unit.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 27086)
+++ /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 27087)
@@ -52,8 +52,8 @@
         cmap = getcolormap(options)
     if options.exist('cmap_set_over'):
-        over = options.getfieldvalue('cmap_set_over', '0.5')
+        over = options.getfieldvalue('cmap_set_over', 'k')
         cmap.set_over(over)
     if options.exist('cmap_set_under'):
-        under = options.getfieldvalue('cmap_set_under', '0.5')
+        under = options.getfieldvalue('cmap_set_under', 'k')
         cmap.set_under(under)
     options.addfield('colormap', cmap)
@@ -75,7 +75,5 @@
         limextent = 0.
 
-    if options.exist('clim'):
-        lims = options.getfieldvalue('clim', [np.nanmin(data), np.nanmax(data)])
-    elif options.exist('caxis'):
+    if options.exist('caxis'):
         lims = options.getfieldvalue('caxis', [np.nanmin(data), np.nanmax(data)])
     else:
@@ -93,8 +91,4 @@
     else:
         norm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1])
-    if options.exist('log'):
-        norm = mpl.colors.LogNorm(vmin=lims[0], vmax=lims[1])
-    else:
-        norm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1])
     options.addfield('colornorm', norm)
     # }}}
@@ -102,5 +96,4 @@
     # Plot depending on the datatype
     # {{{ data are on elements
-
     if datatype == 1:
         if is2d:
Index: /issm/trunk-jpl/src/m/plot/plotmodel.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plotmodel.py	(revision 27086)
+++ /issm/trunk-jpl/src/m/plot/plotmodel.py	(revision 27087)
@@ -86,9 +86,9 @@
         #       taken from https://matplotlib.org/api/_as_gen/mpl_toolkits.axes_grid1.axes_grid.ImageGrid.html
         #
-        direction = options.list[0].getfieldvalue('direction', 'row') # {"row", "column"}, default: "row"
-        axes_pad = options.list[0].getfieldvalue('axes_pad', 0.25) # float or (float, float), default : 0.02; Padding or (horizonal padding, vertical padding) between axes, in inches
-        add_all = options.list[0].getfieldvalue('add_all', True) # bool, default: True
-        share_all = options.list[0].getfieldvalue('share_all', True) # bool, default: False
-        label_mode = options.list[0].getfieldvalue('label_mode', 'L') # {"L", "1", "all"}, default: "L"; Determines which axes will get tick labels: "L": All axes on the left column get vertical tick labels; all axes on the bottom row get horizontal tick labels;. "1": Only the bottom left axes is labelled. "all": all axes are labelled.
+        direction = options.list[0].getfieldvalue('direction', 'row')  # {"row", "column"}, default: "row"
+        axes_pad = options.list[0].getfieldvalue('axes_pad', 0.25)  # float or (float, float), default : 0.02; Padding or (horizonal padding, vertical padding) between axes, in inches
+        add_all = options.list[0].getfieldvalue('add_all', True)  # bool, default: True
+        share_all = options.list[0].getfieldvalue('share_all', True)  # bool, default: False
+        label_mode = options.list[0].getfieldvalue('label_mode', 'L')  # {"L", "1", "all"}, default: "L"; Determines which axes will get tick labels: "L": All axes on the left column get vertical tick labels; all axes on the bottom row get horizontal tick labels;. "1": Only the bottom left axes is labelled. "all": all axes are labelled.
 
         # Translate MATLAB colorbar mode to matplotlib
Index: /issm/trunk-jpl/src/m/plot/processdata.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/processdata.py	(revision 27086)
+++ /issm/trunk-jpl/src/m/plot/processdata.py	(revision 27087)
@@ -60,8 +60,22 @@
     #  quiver plot {{{
     if datasize[1] > 1 and datasize[0] != numberofvertices + 1:
-        if datasize[0] == numberofvertices and datasize[1] == 2:
+        if datasize[0] == numberofvertices and datasize[1] in [2, 3]:
             datatype = 3
+            if md.mesh.dimension() == 3:
+                if datasize[1] == 2:
+                    data = np.hstack(data, np.zeros((datasize[0])))
+                elif datasize[1] > 3:
+                    raise ValueError('plotmodel error message: data should have two or three columns of length md.mesh.numberofvertices for a quiver plot')
         else:
-            raise ValueError('plotmodel error message: data should have two columns of length md.mesh.numberofvertices for a quiver plot')
+            #we should have a patch
+            print("Assuming that data provided is a patch")
+            datatype = 4
+            index = md.mesh.elements
+            if np.shape(data)[1] < np.shape(index)[1]:
+                    raise ValueError('plotmodel error message: data should have more columns than vertices per elements to plot a patch')
+            procdata = np.zeros((numberofvertices))
+            procdata[md.mesh.elements -1] = data[:, 0:np.shape(index)[1]]
+            datasize = [numberofvertices, 1]
+
     # }}}
 
@@ -142,16 +156,24 @@
     nanfill = options.getfieldvalue('nan', -9999)
     if np.any(np.isnan(procdata)):
-        lb = np.nanmin(procdata)
-        ub = np.nanmax(procdata)
-        if lb == ub:
-            lb = lb - 0.5
-            ub = ub + 0.5
-            nanfill = lb - 1
+        if options.exist('caxis'):
+            [lb, ub] = options.getfieldvalue('caxis')
+        else:
+            lb = np.nanmin(procdata)
+            ub = np.nanmax(procdata)
+            if lb == ub:
+                lb = lb - 0.5
+                ub = ub + 0.5
+                nanfill = lb - 1
+            options.addfielddefault('caxis', [lb, ub])
+
         procdata[np.isnan(procdata)] = nanfill
         procdata = np.ma.array(procdata, mask=np.isnan(procdata))
-        #clim looks to be deprecated and replaced by caxis
-        #options.addfielddefault('clim', [lb, ub])
-        options.addfielddefault('cmap_set_under', '1')
-        print(("WARNING: nan's treated as", nanfill, "by default.  Change using pairoption 'nan', nan_fill_value in plotmodel call"))
+        print('from nan processing {} and {}'.format(lb, ub))
+        if nanfill < lb:
+            options.addfielddefault('cmap_set_under', 'k')
+        elif nanfill > ub:
+            options.addfielddefault('cmap_set_over', 'k')
+        if nanfill < ub and nanfill > lb:
+            print(("WARNING: nan's treated as", nanfill, "by default. Which is in your data interval, change it with ['nan', value] in plotmodel options"))
     # }}}