Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 25608)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 25609)
@@ -146,4 +146,5 @@
 syn keyword cConstant FlowequationIsHOEnum
 syn keyword cConstant FlowequationIsL1L2Enum
+syn keyword cConstant FlowequationIsMLHOEnum
 syn keyword cConstant FlowequationIsSIAEnum
 syn keyword cConstant FlowequationIsSSAEnum
@@ -1133,4 +1134,5 @@
 syn keyword cConstant JEnum
 syn keyword cConstant L1L2ApproximationEnum
+syn keyword cConstant MLHOApproximationEnum
 syn keyword cConstant L2ProjectionBaseAnalysisEnum
 syn keyword cConstant L2ProjectionEPLAnalysisEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 25608)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 25609)
@@ -140,4 +140,5 @@
 	FlowequationIsHOEnum,
 	FlowequationIsL1L2Enum,
+	FlowequationIsMLHOEnum,
 	FlowequationIsSIAEnum,
 	FlowequationIsSSAEnum,
@@ -1132,4 +1133,5 @@
 	JEnum,
 	L1L2ApproximationEnum,
+	MLHOApproximationEnum,
 	L2ProjectionBaseAnalysisEnum,
 	L2ProjectionEPLAnalysisEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 25608)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 25609)
@@ -148,4 +148,5 @@
 		case FlowequationIsHOEnum : return "FlowequationIsHO";
 		case FlowequationIsL1L2Enum : return "FlowequationIsL1L2";
+		case FlowequationIsMLHOEnum : return "FlowequationIsMLHO";
 		case FlowequationIsSIAEnum : return "FlowequationIsSIA";
 		case FlowequationIsSSAEnum : return "FlowequationIsSSA";
@@ -1135,4 +1136,5 @@
 		case JEnum : return "J";
 		case L1L2ApproximationEnum : return "L1L2Approximation";
+		case MLHOApproximationEnum : return "MLHOApproximation";
 		case L2ProjectionBaseAnalysisEnum : return "L2ProjectionBaseAnalysis";
 		case L2ProjectionEPLAnalysisEnum : return "L2ProjectionEPLAnalysis";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 25608)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 25609)
@@ -151,4 +151,5 @@
 	      else if (strcmp(name,"FlowequationIsHO")==0) return FlowequationIsHOEnum;
 	      else if (strcmp(name,"FlowequationIsL1L2")==0) return FlowequationIsL1L2Enum;
+	      else if (strcmp(name,"FlowequationIsMLHO")==0) return FlowequationIsMLHOEnum;
 	      else if (strcmp(name,"FlowequationIsSIA")==0) return FlowequationIsSIAEnum;
 	      else if (strcmp(name,"FlowequationIsSSA")==0) return FlowequationIsSSAEnum;
@@ -259,9 +260,9 @@
 	      else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum;
 	      else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
-	      else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
+	      if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
+	      else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
 	      else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
 	      else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
@@ -382,9 +383,9 @@
 	      else if (strcmp(name,"SmbIsd18opd")==0) return SmbIsd18opdEnum;
 	      else if (strcmp(name,"SmbIsdelta18o")==0) return SmbIsdelta18oEnum;
-	      else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum;
+	      if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum;
+	      else if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum;
 	      else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum;
 	      else if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum;
@@ -505,9 +506,9 @@
 	      else if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum;
 	      else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum;
-	      else if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
+	      if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum;
+	      else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
 	      else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
 	      else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
@@ -628,9 +629,9 @@
 	      else if (strcmp(name,"HydrologydcEplThicknessTransient")==0) return HydrologydcEplThicknessTransientEnum;
 	      else if (strcmp(name,"HydrologydcMaskEplactiveElt")==0) return HydrologydcMaskEplactiveEltEnum;
-	      else if (strcmp(name,"HydrologydcMaskEplactiveNode")==0) return HydrologydcMaskEplactiveNodeEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum;
+	      if (strcmp(name,"HydrologydcMaskEplactiveNode")==0) return HydrologydcMaskEplactiveNodeEnum;
+	      else if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum;
 	      else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum;
 	      else if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum;
@@ -751,9 +752,9 @@
 	      else if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum;
 	      else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
-	      else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
+	      if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum;
+	      else if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
 	      else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
 	      else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
@@ -874,9 +875,9 @@
 	      else if (strcmp(name,"Vx")==0) return VxEnum;
 	      else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
-	      else if (strcmp(name,"VxObs")==0) return VxObsEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
+	      if (strcmp(name,"VxObs")==0) return VxObsEnum;
+	      else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
 	      else if (strcmp(name,"Vy")==0) return VyEnum;
 	      else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
@@ -997,9 +998,9 @@
 	      else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
 	      else if (strcmp(name,"AdaptiveTimestepping")==0) return AdaptiveTimesteppingEnum;
-	      else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
+	      if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
+	      else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
 	      else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
 	      else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
@@ -1120,9 +1121,9 @@
 	      else if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum;
 	      else if (strcmp(name,"GroundinglineMassFlux")==0) return GroundinglineMassFluxEnum;
-	      else if (strcmp(name,"Gset")==0) return GsetEnum;
          else stage=10;
    }
    if(stage==10){
-	      if (strcmp(name,"Gsl")==0) return GslEnum;
+	      if (strcmp(name,"Gset")==0) return GsetEnum;
+	      else if (strcmp(name,"Gsl")==0) return GslEnum;
 	      else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
 	      else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum;
@@ -1162,4 +1163,5 @@
 	      else if (strcmp(name,"J")==0) return JEnum;
 	      else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum;
+	      else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum;
 	      else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
 	      else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
@@ -1242,10 +1244,10 @@
 	      else if (strcmp(name,"P0DG")==0) return P0DGEnum;
 	      else if (strcmp(name,"P1DG")==0) return P1DGEnum;
-	      else if (strcmp(name,"P1P1")==0) return P1P1Enum;
-	      else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
          else stage=11;
    }
    if(stage==11){
-	      if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
+	      if (strcmp(name,"P1P1")==0) return P1P1Enum;
+	      else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
+	      else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
 	      else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
 	      else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
@@ -1365,10 +1367,10 @@
 	      else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
 	      else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;
-	      else if (strcmp(name,"EtaAbsGradient")==0) return EtaAbsGradientEnum;
-	      else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
          else stage=12;
    }
    if(stage==12){
-	      if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
+	      if (strcmp(name,"EtaAbsGradient")==0) return EtaAbsGradientEnum;
+	      else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
+	      else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
 	      else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
 	      else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;
Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp	(revision 25608)
+++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp	(revision 25609)
@@ -311,9 +311,10 @@
 		case 2: return SSAApproximationEnum;
 		case 3: return L1L2ApproximationEnum;
-		case 4: return HOApproximationEnum;
-		case 5: return FSApproximationEnum;
-		case 6: return SSAHOApproximationEnum;
-		case 7: return HOFSApproximationEnum;
-		case 8: return SSAFSApproximationEnum;
+		case 4: return MLHOApproximationEnum;
+		case 5: return HOApproximationEnum;
+		case 6: return FSApproximationEnum;
+		case 7: return SSAHOApproximationEnum;
+		case 8: return HOFSApproximationEnum;
+		case 9: return SSAFSApproximationEnum;
 		default: _error_("Marshalled vertex equation code \""<<enum_in<<"\" not supported yet.");
 	}
@@ -325,9 +326,10 @@
 		case 2: return SSAApproximationEnum;
 		case 3: return L1L2ApproximationEnum;
-		case 4: return HOApproximationEnum;
-		case 5: return FSApproximationEnum;
-		case 6: return SSAHOApproximationEnum;
-		case 7: return SSAFSApproximationEnum;
-		case 8: return HOFSApproximationEnum;
+		case 4: return MLHOApproximationEnum;
+		case 5: return HOApproximationEnum;
+		case 6: return FSApproximationEnum;
+		case 7: return SSAHOApproximationEnum;
+		case 8: return SSAFSApproximationEnum;
+		case 9: return HOFSApproximationEnum;
 		default: _error_("Marshalled element equation code \""<<enum_in<<"\" not supported yet.");
 	}
Index: /issm/trunk-jpl/src/m/classes/flowequation.js
===================================================================
--- /issm/trunk-jpl/src/m/classes/flowequation.js	(revision 25608)
+++ /issm/trunk-jpl/src/m/classes/flowequation.js	(revision 25609)
@@ -22,4 +22,5 @@
 		fielddisplay(this,'isSSA','is the Shelfy-Stream Approximation (SSA) used ?');
 		fielddisplay(this,'isL1L2','is the L1L2 approximation used ?');
+		fielddisplay(this,'isMLHO','is the MLHO approximation used ?');
 		fielddisplay(this,'isHO','is the Higher-Order (HO) approximation used ?');
 		fielddisplay(this,'isFS','are the Full-FS (FS) equations used ?');
@@ -56,4 +57,5 @@
 			checkfield(md,'fieldname','flowequation.isSSA','numel',[1],'values',[0, 1]);
 			checkfield(md,'fieldname','flowequation.isL1L2','numel',[1],'values',[0, 1]);
+			checkfield(md,'fieldname','flowequation.isMLHO','numel',[1],'values',[0, 1]);
 			checkfield(md,'fieldname','flowequation.isHO','numel',[1],'values',[0, 1]);
 			checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0, 1]);
@@ -82,10 +84,10 @@
 			}
 			else if (md.mesh.domaintype() =='3D'){
-				checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices, 1],'values',[0,1,2,3,4,5,6,7,8]);
-				checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements, 1],'values',[0,1,2,3,4,5,6,7,8]);
+				checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices, 1],'values',[0,1,2,3,4,5,6,7,8,9]);
+				checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements, 1],'values',[0,1,2,3,4,5,6,7,8,9]);
 			}
 			else throw Error('Case not supported yet');
 			
-			if (!(this.isSIA | this.isSSA | this.isL1L2 | this.isHO | this.isFS)){
+			if (!(this.isSIA | this.isSSA | this.isL1L2 | this.isMLHO | this.isHO | this.isFS)){
 				md = checkmessage(md,['no element types set for this model']);
 			}
@@ -102,4 +104,5 @@
 			WriteData(fid,prefix,'object',this,'fieldname','isSSA','format','Boolean');
 			WriteData(fid,prefix,'object',this,'fieldname','isL1L2','format','Boolean');
+			WriteData(fid,prefix,'object',this,'fieldname','isMLHO','format','Boolean');
 			WriteData(fid,prefix,'object',this,'fieldname','isHO','format','Boolean');
 			WriteData(fid,prefix,'object',this,'fieldname','isFS','format','Boolean');
@@ -129,4 +132,5 @@
 	this.isSSA                          = 0;
 	this.isL1L2                         = 0;
+	this.isMLHO                         = 0;
 	this.isHO                           = 0;
 	this.isFS                           = 0;
Index: /issm/trunk-jpl/src/m/classes/flowequation.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/flowequation.m	(revision 25608)
+++ /issm/trunk-jpl/src/m/classes/flowequation.m	(revision 25609)
@@ -9,4 +9,5 @@
 		isSSA                          = 0;
 		isL1L2                         = 0;
+		isMLHO                         = 0;
 		isHO                           = 0;
 		isFS                           = 0;
@@ -98,4 +99,5 @@
 			md = checkfield(md,'fieldname','flowequation.isSSA','numel',[1],'values',[0 1]);
 			md = checkfield(md,'fieldname','flowequation.isL1L2','numel',[1],'values',[0 1]);
+			md = checkfield(md,'fieldname','flowequation.isMLHO','numel',[1],'values',[0 1]);
 			md = checkfield(md,'fieldname','flowequation.isHO','numel',[1],'values',[0 1]);
 			md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0 1]);
@@ -123,10 +125,10 @@
 				md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[2,4,5]);
 			elseif strcmp(domaintype(md.mesh),'3D'),
-				md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[0:8]);
-				md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[0:8]);
+				md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices 1],'values',[0:9]);
+				md = checkfield(md,'fieldname','flowequation.element_equation','size',[md.mesh.numberofelements 1],'values',[0:9]);
 			else
 				error('Case not supported yet');
 			end
-			if ~(self.isSIA || self.isSSA || self.isL1L2 || self.isHO || self.isFS),
+			if ~(self.isSIA || self.isSSA || self.isL1L2 || self.isMLHO || self.isHO || self.isFS),
 				md = checkmessage(md,['no element types set for this model']);
 			end
@@ -146,4 +148,5 @@
 			fielddisplay(self,'isSSA','is the Shelfy-Stream Approximation (SSA) used?');
 			fielddisplay(self,'isL1L2','is the L1L2 approximation used?');
+			fielddisplay(self,'isMLHO','is the Mono-Layer Higher-Order approximation used?');
 			fielddisplay(self,'isHO','is the Higher-Order (HO) approximation used?');
 			fielddisplay(self,'isFS','are the Full-FS (FS) equations used?');
@@ -164,4 +167,5 @@
 			WriteData(fid,prefix,'object',self,'fieldname','isSSA','format','Boolean');
 			WriteData(fid,prefix,'object',self,'fieldname','isL1L2','format','Boolean');
+			WriteData(fid,prefix,'object',self,'fieldname','isMLHO','format','Boolean');
 			WriteData(fid,prefix,'object',self,'fieldname','isHO','format','Boolean');
 			WriteData(fid,prefix,'object',self,'fieldname','isFS','format','Boolean');
@@ -188,4 +192,5 @@
 			writejsdouble(fid,[modelname '.flowequation.isSSA'],self.isSSA);
 			writejsdouble(fid,[modelname '.flowequation.isL1L2'],self.isL1L2);
+			writejsdouble(fid,[modelname '.flowequation.isMLHO'],self.isMLHO);
 			writejsdouble(fid,[modelname '.flowequation.isHO'],self.isHO);
 			writejsdouble(fid,[modelname '.flowequation.isFS'],self.isFS);
Index: /issm/trunk-jpl/src/m/classes/flowequation.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/flowequation.py	(revision 25608)
+++ /issm/trunk-jpl/src/m/classes/flowequation.py	(revision 25609)
@@ -19,4 +19,5 @@
         self.isSSA = 0
         self.isL1L2 = 0
+        self.isMLHO = 0
         self.isHO = 0
         self.isFS = 0
@@ -47,4 +48,5 @@
         s += "{}\n".format(fielddisplay(self, 'isSSA', "is the Shelfy-Stream Approximation (SSA) used?"))
         s += "{}\n".format(fielddisplay(self, 'isL1L2', "are L1L2 equations used?"))
+        s += "{}\n".format(fielddisplay(self, 'isMLHO', "are Mono-layer Higher-Order equations used?"))
         s += "{}\n".format(fielddisplay(self, 'isHO', "is the Higher-Order (HO) approximation used?"))
         s += "{}\n".format(fielddisplay(self, 'isFS', "are the Full-FS (FS) equations used?"))
@@ -94,4 +96,5 @@
         md = checkfield(md, 'fieldname', 'flowequation.isSSA', 'numel', [1], 'values', [0, 1])
         md = checkfield(md, 'fieldname', 'flowequation.isL1L2', 'numel', [1], 'values', [0, 1])
+        md = checkfield(md, 'fieldname', 'flowequation.isMLHO', 'numel', [1], 'values', [0, 1])
         md = checkfield(md, 'fieldname', 'flowequation.isHO', 'numel', [1], 'values', [0, 1])
         md = checkfield(md, 'fieldname', 'flowequation.isFS', 'numel', [1], 'values', [0, 1])
@@ -116,9 +119,9 @@
             md = checkfield(md, 'fieldname', 'flowequation.element_equation', 'size', [md.mesh.numberofelements], 'values', [2, 4, 5])
         elif m.strcmp(md.mesh.domaintype(), '3D'):
-            md = checkfield(md, 'fieldname', 'flowequation.vertex_equation', 'size', [md.mesh.numberofvertices], 'values', np.arange(0, 8 + 1))
-            md = checkfield(md, 'fieldname', 'flowequation.element_equation', 'size', [md.mesh.numberofelements], 'values', np.arange(0, 8 + 1))
+            md = checkfield(md, 'fieldname', 'flowequation.vertex_equation', 'size', [md.mesh.numberofvertices], 'values', np.arange(0, 9 + 1))
+            md = checkfield(md, 'fieldname', 'flowequation.element_equation', 'size', [md.mesh.numberofelements], 'values', np.arange(0, 9 + 1))
         else:
             raise RuntimeError('mesh type not supported yet')
-        if not (self.isSIA or self.isSSA or self.isL1L2 or self.isHO or self.isFS):
+        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")
 
@@ -135,4 +138,5 @@
         WriteData(fid, prefix, 'object', self, 'fieldname', 'isSSA', 'format', 'Boolean')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'isL1L2', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'isMLHO', 'format', 'Boolean')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'isHO', 'format', 'Boolean')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'isFS', 'format', 'Boolean')
Index: /issm/trunk-jpl/src/m/parameterization/setflowequation.js
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/setflowequation.js	(revision 25608)
+++ /issm/trunk-jpl/src/m/parameterization/setflowequation.js	(revision 25609)
@@ -3,5 +3,5 @@
 //
 //   This routine works like plotmodel: it works with an even number of inputs
-//   'SIA','SSA','L1L2','HO','FS' and 'fill' are the possible options
+//   'SIA','SSA','L1L2','MLHO','HO','FS' and 'fill' are the possible options
 //   that must be followed by the corresponding exp file or flags list
 //   It can either be a domain file (argus type, .exp extension), or an array of element flags. 
@@ -11,5 +11,5 @@
 //   a string 'all' will be considered as the entire domain
 //   You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling'
-//   NB: L1L2 cannot currently be coupled to any other ice flow model
+//   NB: L1L2 and MLHO cannot currently be coupled to any other ice flow model
 //
 //   Usage:
@@ -38,4 +38,5 @@
 	HOflag   = FlagElements(md,options.getfieldvalue('HO',''));
 	L1L2flag = FlagElements(md,options.getfieldvalue('L1L2',''));
+	MLHOflag = FlagElements(md,options.getfieldvalue('MLHO',''));
 	FSflag   = FlagElements(md,options.getfieldvalue('FS',''));
 	filltype = options.getfieldvalue('fill','none');
@@ -54,9 +55,9 @@
 
 	//check that each element has at least one flag
-	for(var i=0;i<md.mesh.numberofelements;i++)if((SIAflag[i] + SSAflag[i] + HOflag[i] + L1L2flag[i] + FSflag[i])==0)
+	for(var i=0;i<md.mesh.numberofelements;i++)if((SIAflag[i] + SSAflag[i] + HOflag[i] + L1L2flag[i] + MLHOflag[i] + FSflag[i])==0)
 	throw Error("elements type not assigned, supported models are 'SIA','SSA','HO' and 'FS'");
 
 	//check that each element has only one flag
-	if (ArrayAnyAboveStrict(ArrayXPY(SIAflag,SSAflag,HOflag,L1L2flag),1)){
+	if (ArrayAnyAboveStrict(ArrayXPY(SIAflag,SSAflag,HOflag,L1L2flag,MLHOflag),1)){
 		console.log('setflowequation warning message: some elements have several types, higher order type is used for them')
 
@@ -70,4 +71,5 @@
 	//check that L1L2 is not coupled to any other model for now
 	if (ArrayAnyEqual(L1L2flag,1) & ArrayAnyEqual(ArrayOr(SIAflag,SSAflag,HOflag,FSflag),1)) throw Error('L1L2 cannot be coupled to any other model');
+	if (ArrayAnyEqual(MLHOflag,1) & ArrayAnyEqual(ArrayOr(SIAflag,SSAflag,HOflag,FSflag),1)) throw Error('MLHO cannot be coupled to any other model');
 
 	//Check that no HO or FS for 2d mesh
@@ -97,4 +99,8 @@
 	pos=ArrayFind(L1L2flag,1);
 	for(var i=0;i<pos.length;i++) for(var j=0;j<md.mesh.elements[0].length;j++) nodeonL1L2[md.mesh.elements[pos[i]][j]-1]=1;
+
+	nodeonMLHO=NewArrayFill(md.mesh.numberofvertices,0);
+	pos=ArrayFind(MLHOflag,1);
+	for(var i=0;i<pos.length;i++) for(var j=0;j<md.mesh.elements[0].length;j++) nodeonMLHO[md.mesh.elements[pos[i]][j]-1]=1;
 
 	nodeonFS=NewArrayFill(md.mesh.numberofvertices,0);
@@ -256,9 +262,10 @@
 	pos=ArrayFind(SSAflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=2;
 	pos=ArrayFind(L1L2flag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=3;
-	pos=ArrayFind(HOflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=4;
-	pos=ArrayFind(FSflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=5;
-	pos=ArrayFind(SSAHOflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=6;
-	pos=ArrayFind(SSAFSflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=7;
-	pos=ArrayFind(HOFSflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=8;
+	pos=ArrayFind(MLHOflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=4;
+	pos=ArrayFind(HOflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=5;
+	pos=ArrayFind(FSflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=6;
+	pos=ArrayFind(SSAHOflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=7;
+	pos=ArrayFind(SSAFSflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=8;
+	pos=ArrayFind(HOFSflag,1);for(var i=0;i<pos.length;i++)md.flowequation.element_equation[pos[i]]=9;
 
 
@@ -274,6 +281,7 @@
 	pos=ArrayFind(nodeonSSA,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=2;
 	pos=ArrayFind(nodeonL1L2,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=3;
-	pos=ArrayFind(nodeonHO,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=4;
-	pos=ArrayFind(nodeonFS,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=5;
+	pos=ArrayFind(nodeonMLHO,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=4;
+	pos=ArrayFind(nodeonHO,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=5;
+	pos=ArrayFind(nodeonFS,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=6;
 	//DO SIA LAST! Otherwise spcs might not be set up correctly (SIA should have priority)
 	pos=ArrayFind(nodeonSIA,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=1;
@@ -285,7 +293,7 @@
 	}
 
-	pos=ArrayFind(nodeonSSAHO,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=6;
-	pos=ArrayFind(nodeonHOFS,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=7;
-	pos=ArrayFind(nodeonSSAFS,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=8;
+	pos=ArrayFind(nodeonSSAHO,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=7;
+	pos=ArrayFind(nodeonHOFS,1);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=8;
+	pos=ArrayFind(nodeonSSAFS,2);for(var i=0;i<pos.length;i++)md.flowequation.vertex_equation[pos[i]]=9;
 
 	//figure out solution types
@@ -293,6 +301,7 @@
 	md.flowequation.isSSA  = ArrayAnyEqual(md.flowequation.element_equation,2);
 	md.flowequation.isL1L2 = ArrayAnyEqual(md.flowequation.element_equation,3);
-	md.flowequation.isHO   = ArrayAnyEqual(md.flowequation.element_equation,4);
-	md.flowequation.isFS   = ArrayAnyEqual(md.flowequation.element_equation,5);
+	md.flowequation.isMLHO = ArrayAnyEqual(md.flowequation.element_equation,4);
+	md.flowequation.isHO   = ArrayAnyEqual(md.flowequation.element_equation,5);
+	md.flowequation.isFS   = ArrayAnyEqual(md.flowequation.element_equation,6);
 	return
 
Index: /issm/trunk-jpl/src/m/parameterization/setflowequation.m
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/setflowequation.m	(revision 25608)
+++ /issm/trunk-jpl/src/m/parameterization/setflowequation.m	(revision 25609)
@@ -3,5 +3,5 @@
 %
 %   This routine works like plotmodel: it works with an even number of inputs
-%   'SIA','SSA','L1L2','HO','FS' and 'fill' are the possible options
+%   'SIA','SSA','L1L2','MLHO','HO','FS' and 'fill' are the possible options
 %   that must be followed by the corresponding exp file or flags list
 %   It can either be a domain file (argus type, .exp extension), or an array of element flags. 
@@ -11,5 +11,5 @@
 %   a string 'all' will be considered as the entire domain
 %   You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling'
-%   NB: L1L2 cannot currently be coupled to any other ice flow model
+%   NB: L1L2 and MLHO cannot currently be coupled to any other ice flow model
 %
 %   Usage:
@@ -39,4 +39,5 @@
 HOflag   = FlagElements(md,getfieldvalue(options,'HO',''));
 L1L2flag = FlagElements(md,getfieldvalue(options,'L1L2',''));
+MLHOflag = FlagElements(md,getfieldvalue(options,'MLHO',''));
 FSflag   = FlagElements(md,getfieldvalue(options,'FS',''));
 filltype = getfieldvalue(options,'fill','none');
@@ -53,10 +54,10 @@
 
 %check that each element has at least one flag
-if any(SIAflag+SSAflag+HOflag+L1L2flag+FSflag==0),
+if any(SIAflag+SSAflag+HOflag+L1L2flag+MLHOflag+FSflag==0),
 	error('elements type not assigned, supported models are ''SIA'',''SSA'',''HO'' and ''FS''')
 end
 
 %check that each element has only one flag
-if any(SIAflag+SSAflag+HOflag+L1L2flag+FSflag>1),
+if any(SIAflag+SSAflag+HOflag+L1L2flag+MLHOflag+FSflag>1),
 	disp('setflowequation warning message: some elements have several types, higher order type is used for them')
 	SIAflag(find(SIAflag & SSAflag))=0;
@@ -69,4 +70,7 @@
 	error('L1L2 cannot be coupled to any other model');
 end
+if any(MLHOflag) & any(SIAflag | SSAflag | HOflag | FSflag)
+	error('MLHO cannot be coupled to any other model');
+end
 
 %Check that no HO or FS for 2d mesh
@@ -83,12 +87,9 @@
 
 %Initialize node fields
-nodeonSIA=zeros(md.mesh.numberofvertices,1);
-nodeonSIA(md.mesh.elements(find(SIAflag),:))=1;
-nodeonSSA=zeros(md.mesh.numberofvertices,1);
-nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
-nodeonHO=zeros(md.mesh.numberofvertices,1);
-nodeonHO(md.mesh.elements(find(HOflag),:))=1;
-nodeonL1L2=zeros(md.mesh.numberofvertices,1);
-nodeonL1L2(md.mesh.elements(find(L1L2flag),:))=1;
+nodeonSIA=zeros(md.mesh.numberofvertices,1);  nodeonSIA(md.mesh.elements(find(SIAflag),:))=1;
+nodeonSSA=zeros(md.mesh.numberofvertices,1);  nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
+nodeonHO=zeros(md.mesh.numberofvertices,1);   nodeonHO(md.mesh.elements(find(HOflag),:))=1;
+nodeonL1L2=zeros(md.mesh.numberofvertices,1); nodeonL1L2(md.mesh.elements(find(L1L2flag),:))=1;
+nodeonMLHO=zeros(md.mesh.numberofvertices,1); nodeonMLHO(md.mesh.elements(find(MLHOflag),:))=1;
 nodeonFS=zeros(md.mesh.numberofvertices,1);
 noneflag=zeros(md.mesh.numberofelements,1);
@@ -247,9 +248,10 @@
 md.flowequation.element_equation(find(SSAflag))=2;
 md.flowequation.element_equation(find(L1L2flag))=3;
-md.flowequation.element_equation(find(HOflag))=4;
-md.flowequation.element_equation(find(FSflag))=5;
-md.flowequation.element_equation(find(SSAHOflag))=6;
-md.flowequation.element_equation(find(SSAFSflag))=7;
-md.flowequation.element_equation(find(HOFSflag))=8;
+md.flowequation.element_equation(find(MLHOflag))=4;
+md.flowequation.element_equation(find(HOflag))=5;
+md.flowequation.element_equation(find(FSflag))=6;
+md.flowequation.element_equation(find(SSAHOflag))=7;
+md.flowequation.element_equation(find(SSAFSflag))=8;
+md.flowequation.element_equation(find(HOFSflag))=9;
 
 %border
@@ -260,12 +262,9 @@
 %Create vertices_type
 md.flowequation.vertex_equation=zeros(md.mesh.numberofvertices,1);
-pos=find(nodeonSSA);
-md.flowequation.vertex_equation(pos)=2;
-pos=find(nodeonL1L2);
-md.flowequation.vertex_equation(pos)=3;
-pos=find(nodeonHO);
-md.flowequation.vertex_equation(pos)=4;
-pos=find(nodeonFS);
-md.flowequation.vertex_equation(pos)=5;
+pos=find(nodeonSSA);  md.flowequation.vertex_equation(pos)=2;
+pos=find(nodeonL1L2); md.flowequation.vertex_equation(pos)=3;
+pos=find(nodeonMLHO); md.flowequation.vertex_equation(pos)=4;
+pos=find(nodeonHO);   md.flowequation.vertex_equation(pos)=5;
+pos=find(nodeonFS);   md.flowequation.vertex_equation(pos)=6;
 %DO SIA LAST! Otherwise spcs might not be set up correctly (SIA should have priority)
 pos=find(nodeonSIA);
@@ -278,9 +277,9 @@
 end
 pos=find(nodeonSSAHO);
-md.flowequation.vertex_equation(pos)=6;
+md.flowequation.vertex_equation(pos)=7;
 pos=find(nodeonHOFS);
-md.flowequation.vertex_equation(pos)=7;
+md.flowequation.vertex_equation(pos)=8;
 pos=find(nodeonSSAFS);
-md.flowequation.vertex_equation(pos)=8;
+md.flowequation.vertex_equation(pos)=9;
 
 %figure out solution types
@@ -288,6 +287,7 @@
 md.flowequation.isSSA  = double(any(md.flowequation.element_equation == 2));
 md.flowequation.isL1L2 = double(any(md.flowequation.element_equation == 3));
-md.flowequation.isHO   = double(any(md.flowequation.element_equation == 4));
-md.flowequation.isFS   = double(any(md.flowequation.element_equation == 5));
+md.flowequation.isMLHO = double(any(md.flowequation.element_equation == 4));
+md.flowequation.isHO   = double(any(md.flowequation.element_equation == 5));
+md.flowequation.isFS   = double(any(md.flowequation.element_equation == 6));
 
 return
Index: /issm/trunk-jpl/src/m/parameterization/setflowequation.py
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/setflowequation.py	(revision 25608)
+++ /issm/trunk-jpl/src/m/parameterization/setflowequation.py	(revision 25609)
@@ -10,5 +10,5 @@
 
        This routine works like plotmodel: it works with an even number of inputs
-       'SIA', 'SSA', 'HO', 'L1L2', 'FS' and 'fill' are the possible options
+       'SIA', 'SSA', 'HO', 'L1L2', 'MLHO', 'FS' and 'fill' are the possible options
        that must be followed by the corresponding exp file or flags list
        It can either be a domain file (argus type, .exp extension), or an array of element flags.
@@ -44,4 +44,5 @@
     HOflag = FlagElements(md, options.getfieldvalue('HO', ''))
     L1L2flag = FlagElements(md, options.getfieldvalue('L1L2', ''))
+    MLHOflag = FlagElements(md, options.getfieldvalue('MLHO', ''))
     FSflag = FlagElements(md, options.getfieldvalue('FS', ''))
     filltype = options.getfieldvalue('fill', 'none')
@@ -55,9 +56,9 @@
         HOflag = ~SIAflag & ~SSAflag & ~FSflag
     #check that each element has at least one flag
-    if not any(SIAflag + SSAflag + L1L2flag + HOflag + FSflag):
+    if not any(SIAflag + SSAflag + L1L2flag + MLHOflag + HOflag + FSflag):
         raise TypeError("elements type not assigned, supported models are 'SIA', 'SSA', 'HO' and 'FS'")
 
     #check that each element has only one flag
-    if any(SIAflag + SSAflag + L1L2flag + HOflag + FSflag > 1):
+    if any(SIAflag + SSAflag + L1L2flag + MLHOflag + HOflag + FSflag > 1):
         print("setflowequation warning message: some elements have several types, higher order type is used for them")
         SIAflag[np.where(np.logical_and(SIAflag, SSAflag))] = False
@@ -65,7 +66,9 @@
         SSAflag[np.where(np.logical_and(SSAflag, HOflag))] = False
 
-        #check that L1L2 is not coupled to any other model for now
+        #check that L1L2 and MLHO is not coupled to any other model for now
         if any(L1L2flag) and any(SIAflag + SSAflag + HOflag + FSflag):
             raise TypeError('L1L2 cannot be coupled to any other model')
+        if any(MLHOflag) and any(SIAflag + SSAflag + HOflag + FSflag):
+            raise TypeError('MLHO cannot be coupled to any other model')
 
         #Check that no HO or FS for 2d mesh
@@ -85,4 +88,6 @@
     nodeonL1L2 = np.zeros(md.mesh.numberofvertices, bool)
     nodeonL1L2[md.mesh.elements[np.where(L1L2flag), :] - 1] = True
+    nodeonMLHO = np.zeros(md.mesh.numberofvertices, bool)
+    nodeonMLHO[md.mesh.elements[np.where(MLHOflag), :] - 1] = True
     nodeonHO = np.zeros(md.mesh.numberofvertices, bool)
     nodeonHO[md.mesh.elements[np.where(HOflag), :] - 1] = True
@@ -235,9 +240,10 @@
     md.flowequation.element_equation[np.where(SSAflag)] = 2
     md.flowequation.element_equation[np.where(L1L2flag)] = 3
-    md.flowequation.element_equation[np.where(HOflag)] = 4
-    md.flowequation.element_equation[np.where(FSflag)] = 5
-    md.flowequation.element_equation[np.where(SSAHOflag)] = 6
-    md.flowequation.element_equation[np.where(SSAFSflag)] = 7
-    md.flowequation.element_equation[np.where(HOFSflag)] = 8
+    md.flowequation.element_equation[np.where(MLHOflag)] = 4
+    md.flowequation.element_equation[np.where(HOflag)] = 5
+    md.flowequation.element_equation[np.where(FSflag)] = 6
+    md.flowequation.element_equation[np.where(SSAHOflag)] = 7
+    md.flowequation.element_equation[np.where(SSAFSflag)] = 8
+    md.flowequation.element_equation[np.where(HOFSflag)] = 9
 
     #border
@@ -252,8 +258,10 @@
     pos = np.where(nodeonL1L2)
     md.flowequation.vertex_equation[pos] = 3
+    pos = np.where(nodeonMLHO)
+    md.flowequation.vertex_equation[pos] = 4
     pos = np.where(nodeonHO)
-    md.flowequation.vertex_equation[pos] = 4
+    md.flowequation.vertex_equation[pos] = 5
     pos = np.where(nodeonFS)
-    md.flowequation.vertex_equation[pos] = 5
+    md.flowequation.vertex_equation[pos] = 6
     #DO SIA LAST! Otherwise spcs might not be set up correctly (SIA should have priority)
     pos = np.where(nodeonSIA)
@@ -264,16 +272,17 @@
             md.flowequation.vertex_equation[pos] = 0
     pos = np.where(nodeonSSAHO)
-    md.flowequation.vertex_equation[pos] = 6
+    md.flowequation.vertex_equation[pos] = 7
     pos = np.where(nodeonHOFS)
-    md.flowequation.vertex_equation[pos] = 7
+    md.flowequation.vertex_equation[pos] = 8
     pos = np.where(nodeonSSAFS)
-    md.flowequation.vertex_equation[pos] = 8
+    md.flowequation.vertex_equation[pos] = 9
 
     #figure out solution types
     md.flowequation.isSIA = any(md.flowequation.element_equation == 1)
     md.flowequation.isSSA = any(md.flowequation.element_equation == 2)
-    md.flowequation.isL1L2 = any(md.flowequation.element_equation == 3)
-    md.flowequation.isHO = any(md.flowequation.element_equation == 4)
-    md.flowequation.isFS = any(md.flowequation.element_equation == 5)
+    md.flowequation.isL1L2= any(md.flowequation.element_equation == 3)
+    md.flowequation.isMLHO= any(md.flowequation.element_equation == 4)
+    md.flowequation.isHO = any(md.flowequation.element_equation == 5)
+    md.flowequation.isFS = any(md.flowequation.element_equation == 6)
 
     return md
Index: /issm/trunk-jpl/src/m/plot/plot_elementstype.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_elementstype.m	(revision 25608)
+++ /issm/trunk-jpl/src/m/plot/plot_elementstype.m	(revision 25609)
@@ -16,101 +16,26 @@
 %plot
 subplot(width,width,i);
+p = [];
 
 if is2d
-	pos=find(data==0);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
-	p1=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',0,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==1);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
-	p2=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',2,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==2);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
-	p3=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',3,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==3);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
-	p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',4,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==4);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
-	p5=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',5,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==5);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
-	p6=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',6,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==6);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
-	p7=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',7,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==7);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
-	p8=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',8,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==8);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
-	p9=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',9,'FaceColor','flat','EdgeColor',edgecolor);
+	for i=0:9
+		pos=find(data==i);
+		A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); 
+		pnew=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData',i,'FaceColor','flat','EdgeColor',edgecolor);
+		p = [p;pnew];
+	end
 else
-	pos=find(data==0);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
-	p1=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 0,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==1);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
-	p2=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 1,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 1,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 1,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 1,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 1,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==2);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
-	p3=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 2,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 2,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 2,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 2,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 2,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==3);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
-	p4=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==3);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
-	p5=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 3,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==4);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
-	p6=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 4,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 4,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 4,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 4,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 4,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==5);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
-	p7=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 5,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 5,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 5,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 5,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 5,'FaceColor','flat','EdgeColor',edgecolor);
-	%HOFS elements
-	pos=find(data==7);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
-	p8=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 7,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 7,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 7,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 7,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 7,'FaceColor','flat','EdgeColor',edgecolor);
-	pos=find(data==6);
-	A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
-	p9=patch( 'Faces', [A B C],'Vertices', [x y z],'CData', 6,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', 6,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', 6,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', 6,'FaceColor','flat','EdgeColor',edgecolor);
-	patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', 6,'FaceColor','flat','EdgeColor',edgecolor);
+	for i=0:9
+		pos=find(data==i);
+		A=elements(pos,1); B=elements(pos,2); C=elements(pos,3); D=elements(pos,4); E=elements(pos,5); F=elements(pos,6);
+		pnew = patch( 'Faces', [A B C],  'Vertices', [x y z],'CData', i,'FaceColor','flat','EdgeColor',edgecolor);
+		p = [p;pnew];
+		patch( 'Faces', [D E F],  'Vertices', [x y z],'CData', i,'FaceColor','flat','EdgeColor',edgecolor);
+		patch( 'Faces', [A B E D],'Vertices', [x y z],'CData', i,'FaceColor','flat','EdgeColor',edgecolor);
+		patch( 'Faces', [B E F C],'Vertices', [x y z],'CData', i,'FaceColor','flat','EdgeColor',edgecolor);
+		patch( 'Faces', [C A D F],'Vertices', [x y z],'CData', i,'FaceColor','flat','EdgeColor',edgecolor);
+	end
 end
-legend([p1 p2 p3 p4 p5 p6 p7 p8 p9],...
-		'None','SIA','SSA','L1L2','HO',...
+legend(p,'None','SIA','SSA','L1L2','MLHO','HO',...
 		'SSAHO','FS','SSAFS','HOFS');