Index: /issm/trunk-jpl/src/m/classes/flowequation.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/flowequation.py	(revision 26445)
+++ /issm/trunk-jpl/src/m/classes/flowequation.py	(revision 26446)
@@ -27,9 +27,9 @@
         self.fe_HO = ''
         self.fe_FS = ''
-        self.augmented_lagrangian_r = 1.
-        self.augmented_lagrangian_rhop = 1.
-        self.augmented_lagrangian_rlambda = 1.
-        self.augmented_lagrangian_rholambda = 1.
-        self.XTH_theta = 0.
+        self.augmented_lagrangian_r = 1
+        self.augmented_lagrangian_rhop = 1
+        self.augmented_lagrangian_rlambda = 1
+        self.augmented_lagrangian_rholambda = 1
+        self.XTH_theta = 0
         self.vertex_equation = np.nan
         self.element_equation = np.nan
@@ -107,6 +107,9 @@
         md = checkfield(md, 'fieldname', 'flowequation.XTH_theta', 'numel', [1], '>=', 0., '<', .5)
         if m.strcmp(md.mesh.domaintype(), '2Dhorizontal'):
-            md = checkfield(md, 'fieldname', 'flowequation.vertex_equation', 'size', [md.mesh.numberofvertices], 'values', [1, 2])
-            md = checkfield(md, 'fieldname', 'flowequation.element_equation', 'size', [md.mesh.numberofelements], 'values', [1, 2])
+            md = checkfield(md, 'fieldname', 'flowequation.vertex_equation', 'size', [md.mesh.numberofvertices], 'values', [1, 2, 4])
+            md = checkfield(md, 'fieldname', 'flowequation.element_equation', 'size', [md.mesh.numberofelements], 'values', [1, 2, 4])
+        elif m.strcmp(md.mesh.domaintype(), '3Dsurface'):
+            md = checkfield(md, 'fieldname', 'flowequation.vertex_equation', 'size', [md.mesh.numberofvertices], 'values', np.arange(1, 2 + 1))
+            md = checkfield(md, 'fieldname', 'flowequation.element_equation', 'size', [md.mesh.numberofelements], 'values', np.arange(1, 2 + 1))
         elif m.strcmp(md.mesh.domaintype(), '2Dvertical'):
             md = checkfield(md, 'fieldname', 'flowequation.vertex_equation', 'size', [md.mesh.numberofvertices], 'values', [2, 5, 6])
@@ -117,4 +120,5 @@
         else:
             raise RuntimeError('Case not supported yet')
+
         if not (self.isSIA or self.isSSA or self.isL1L2 or self.isMLHO or self.isHO or self.isFS):
             md.checkmessage("no element types set for this model")
