Index: /issm/trunk-jpl/src/m/parameterization/setflowequation.py
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/setflowequation.py	(revision 24289)
+++ /issm/trunk-jpl/src/m/parameterization/setflowequation.py	(revision 24290)
@@ -46,4 +46,5 @@
     FSflag = FlagElements(md, options.getfieldvalue('FS', ''))
     filltype = options.getfieldvalue('fill', 'none')
+
     #Flag the elements that have not been flagged as filltype
     if 'SIA' in filltype:
@@ -91,5 +92,5 @@
     #First modify FSflag to get rid of elements contrained everywhere (spc + border with HO or SSA)
     if any(FSflag):
-        fullspcnodes = np.logical_or(~np.isnan(md.stressbalance.spcvx) + ~np.isnan(md.stressbalance.spcvy) + ~np.isnan(md.stressbalance.spcvz), np.logical_and(nodeonHO, nodeonFS))  #find all the nodes on the boundary of the domain without icefront
+        fullspcnodes = np.logical_or(~np.isnan(md.stressbalance.spcvx) + ~np.isnan(md.stressbalance.spcvy) + ~np.isnan(md.stressbalance.spcvz) == 3, np.logical_and(nodeonHO, nodeonFS))  #find all the nodes on the boundary of the domain without icefront
         fullspcelems = np.sum(fullspcnodes[md.mesh.elements - 1], axis=1) == 6  #find all the nodes on the boundary of the domain without icefront
         FSflag[np.where(fullspcelems.reshape(-1))] = False
Index: /issm/trunk-jpl/src/m/plot/plot_BC.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_BC.py	(revision 24289)
+++ /issm/trunk-jpl/src/m/plot/plot_BC.py	(revision 24290)
@@ -46,4 +46,8 @@
 
         for key in spc_dict:
+            mark = spc_dict[str(key)][1]
+            color = spc_dict[str(key)][2]
+            size = spc_dict[str(key)][3]
+            name = spc_dict[str(key)][4]
             #first reduce vectors if layer is used
             if options.getfieldvalue('layer', 0) >= 1:
@@ -52,11 +56,10 @@
                 fulldata = md.__dict__[str(spc_dict[str(key)][0])].__dict__[str(key)]
                 data = fulldata[(plotlayer - 1) * slicesize:plotlayer * slicesize]
-                mark = spc_dict[str(key)][1]
-                color = spc_dict[str(key)][2]
-                size = spc_dict[str(key)][3]
-                name = spc_dict[str(key)][4]
-                ax.scatter(x[np.where(~np.isnan(data))],
-                           y[np.where(~np.isnan(data))],
-                           marker=mark, c=color, s=size, label=name, linewidth=0)
+            else:
+                data = md.__dict__[str(spc_dict[str(key)][0])].__dict__[str(key)]
+            ax.scatter(x[np.where(~np.isnan(data))],
+                       y[np.where(~np.isnan(data))],
+                       marker=mark, c=color, s=size, label=name, linewidth=0)
+
     ax.set_xlim(XLims)
     ax.set_ylim(YLims)
Index: /issm/trunk-jpl/src/m/plot/plot_icefront.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_icefront.py	(revision 24289)
+++ /issm/trunk-jpl/src/m/plot/plot_icefront.py	(revision 24290)
@@ -16,4 +16,7 @@
     x, y, z, elements, is2d, isplanet = processmesh(md, [], options)
 
+    icenodes = md.mask.ice_levelset < 0
+    iceelement = np.sum(icenodes[elements], axis=1)
+
     if options.exist('layer'):
         nodes_per_elt = np.shape(md.mesh.elements2d)[1]
@@ -21,5 +24,5 @@
         nodes_per_elt = np.shape(md.mesh.elements)[1]
     #icefront check
-    icefront = np.where(np.abs(np.sum(md.mask.ice_levelset[elements], 1)) != nodes_per_elt)
+    icefront = np.where(np.logical_and(iceelement != nodes_per_elt, iceelement != 0))
 
     #plot mesh