Index: /issm/trunk-jpl/src/m/plot/plot_BC.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_BC.py	(revision 27087)
+++ /issm/trunk-jpl/src/m/plot/plot_BC.py	(revision 27088)
@@ -4,4 +4,5 @@
 from plot_icefront import plot_icefront
 from hydrologydc import hydrologydc
+from hydrologyglads import hydrologyglads
 from mpl_toolkits.mplot3d import Axes3D
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
@@ -36,5 +37,6 @@
         #define what to plot with plot style
         spc_dict = {'spcvx': ['stressbalance', 'o', 'r', 240, 'vx Dirichlet'],
-                    'spcvy': ['stressbalance', 'o', 'b', 160, 'vy Dirichlet']}
+                    'spcvy': ['stressbalance', 'o', 'b', 160, 'vy Dirichlet'],
+                    'spcthickness': ['masstransport', 'o', 'k', 40, 'Thickness']}
         if not is2d:
             spc_dict['spcvz'] = ['stressbalance', 'o', 'y', 80, 'vy Dirichlet']
@@ -44,4 +46,7 @@
             if md.hydrology.isefficientlayer:
                 spc_dict['spcsediment_head'] = ['hydrology', '^', 'b', 240, 'IDS Head']
+
+        if isinstance(md.hydrology, hydrologyglads):
+            spc_dict['spcphi'] = ['hydrology', 'v', 'r', 240, 'phi']
 
         for key in spc_dict:
Index: /issm/trunk-jpl/src/m/plot/plot_edgeoverlay.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_edgeoverlay.py	(revision 27088)
+++ /issm/trunk-jpl/src/m/plot/plot_edgeoverlay.py	(revision 27088)
@@ -0,0 +1,44 @@
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+from processmesh import processmesh
+
+
+def plot_edgeoverlay(md, datain, options, ax):
+    '''
+    plot_channels - plot channels area for GLADS
+    Usage:
+        plot_channels(md, options, fig, axes)
+
+    See also: PLOTMODEL'''
+
+    x, y, z, elements, is2d, isplanet = processmesh(md, [], options)
+
+    flags = datain > 6.7e-5  # this is appropriate for channel Area (perhaps)
+    hide = np.invert(flags)
+
+    NodeMask = np.zeros(np.shape(md.mesh.x), dtype=bool)
+    HideElt = md.mesh.edges[np.where(hide), 0] - 1
+    NodeMask[HideElt] = True
+    MaskX = np.ma.array(x, mask=NodeMask)
+    MaskY = np.ma.array(y, mask=NodeMask)
+
+    EdgeEnd = md.mesh.edges[:, 1] - 1
+    EdgeStart = md.mesh.edges[:, 0] - 1
+    quiverU = MaskX[EdgeEnd] - MaskX[EdgeStart]
+    quiverV = MaskY[EdgeEnd] - MaskY[EdgeStart]
+
+    Masked = np.ma.masked_array(datain, mask=hide)
+    if options.exist('cedgelim'):
+        lims = options.getfieldvalue('cedgelim', [Masked.min(), Masked.max()])
+        edgenorm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1])
+    else:
+        edgenorm = mpl.colors.Normalize(vmin=Masked.min(), vmax=Masked.max())
+    edgemap = plt.cm.get_cmap('inferno')
+
+    ax.quiver(MaskX[EdgeStart], MaskY[EdgeStart], quiverU, quiverV, datain,
+              units="xy", angles="xy", scale_units="xy", scale=1,
+              headwidth=0, headlength=0, width=100, headaxislength=0,
+              norm=edgenorm, cmap=edgemap)
+
+    plt.colorbar(plt.cm.ScalarMappable(norm=edgenorm, cmap=edgemap), ax=ax)
Index: /issm/trunk-jpl/src/m/plot/plot_icefront.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_icefront.py	(revision 27087)
+++ /issm/trunk-jpl/src/m/plot/plot_icefront.py	(revision 27088)
@@ -1,4 +1,5 @@
 import numpy as np
 from processmesh import processmesh
+import matplotlib as mpl
 from mpl_toolkits.mplot3d.art3d import Line3DCollection
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
@@ -26,4 +27,10 @@
     icefront = np.where(np.logical_and(iceelement != nodes_per_elt, iceelement != 0))
 
+    oceannodes = md.mask.ocean_levelset < 0
+    oceanelement = np.sum(oceannodes[elements], axis=1)
+
+    #icefront check
+    groundingline = np.where(np.logical_and(oceanelement != nodes_per_elt, oceanelement != 0))
+
     #plot mesh
     if is2d:
@@ -32,6 +39,11 @@
         #highlight elements on neumann
         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')
+            colors = np.ones(np.shape(elements[icefront])[0])
+            cmap = mpl.colors.ListedColormap("navy")
+            ax.tripcolor(x, y, elements[icefront], facecolors=colors, edgecolor='k', label='elements on ice front', cmap=cmap)
+        if len(groundingline[0]) > 0:
+            colors = np.ones(np.shape(elements[groundingline])[0])
+            cmap = mpl.colors.ListedColormap("limegreen")
+            ax.tripcolor(x, y, elements[groundingline], facecolors=colors, edgecolor='k', label='elements on grounding line', cmap=cmap)
     else:
         ax = inset_axes(ax, width='100%', height='100%', loc=3, borderpad=0, axes_class=Axes3D)