Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 26530)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 26531)
@@ -179,38 +179,4 @@
 		else _error_("not implemented yet");
 	}
-
-}
-/*}}}*/
-void       Penta::Recover3DMLHOInput(int targetVel_enum, int numnodes, IssmDouble* vb,  IssmDouble* vsh, IssmDouble* n, IssmDouble* H, IssmDouble* s){/*{{{*/
-   /* Recover the velocity acording to v=vb+(1-\zeta^{n+1})vsh, where \zeta=(s-z)/H
-    * The variables vb, vsh, n, H and s are all from the 2D horizontal mesh(Tria), with "numnodes" DOFs
-    * To project to penta the DOFs are doubled in size
-    *
-    */
-   _assert_(this->inputs);
-   if(!IsOnBase()) return;
-   else{
-      if(targetVel_enum==VxEnum || targetVel_enum==VyEnum){
-         IssmDouble vel[NUMVERTICES];
-         IssmDouble* xyz_list = NULL;
-         Penta* penta = this;
-         _assert_(NUMVERTICES-2*numnodes==0);
-
-         for(;;){
-            penta->GetVerticesCoordinates(&xyz_list);
-            for(int i=0;i<numnodes;i++) {
-               vel[i] = vb[i] + vsh[i]*(1.0-pow((s[i]-xyz_list[i*3+2])/H[i], (n[i]+1)));
-               vel[i+NUMVERTICES2D] = vb[i] + vsh[i]*(1-pow((s[i]-xyz_list[(i+NUMVERTICES2D)*3+2])/H[i], (n[i]+1)));
-            }
-				xDelete<IssmDouble>(xyz_list);
-
-            /*Add to the bottom side of the element*/
-            penta->AddInput(targetVel_enum,&vel[0],P1Enum);
-            if (penta->IsOnSurface()) break;
-            penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
-         }
-      }
-      else _error_("not implemented yet");
-   }
 
 }
@@ -3298,4 +3264,45 @@
 int        Penta::PressureInterpolation(void){/*{{{*/
 	return PentaRef::PressureInterpolation(this->element_type);
+}
+/*}}}*/
+void       Penta::Recover3DMLHOInput(int targetVel_enum, int numnodes, IssmDouble* vb,  IssmDouble* vsh, IssmDouble* n, IssmDouble* H, IssmDouble* s){/*{{{*/
+   /* Recover the velocity acording to v=vb+(1-\zeta^{n+1})vsh, where \zeta=(s-z)/H
+    * The variables vb, vsh, n, H and s are all from the 2D horizontal mesh(Tria), with "numnodes" DOFs
+    * To project to penta the DOFs are doubled in size
+    *
+    */
+   _assert_(this->inputs);
+   if(!IsOnBase()) return;
+   else{
+      if(targetVel_enum==VxEnum || targetVel_enum==VyEnum){
+         IssmDouble vel[NUMVERTICES2D*5];
+         IssmDouble* xyz_list = NULL;
+			IssmDouble zi;
+         Penta* penta = this;
+			numnodes = penta->NumberofNodes(P1xP4Enum);
+         _assert_(NUMVERTICES2D*5-numnodes==0);
+			GaussPenta gauss;
+
+         for(;;){
+            penta->GetVerticesCoordinates(&xyz_list);
+            for(int i=0;i<NUMVERTICES2D;i++) {
+					for (int j=0;j<5;j++){
+						/* Get z-coordinate of the current node */
+						gauss.GaussNode(P1xP4Enum, i+j*NUMVERTICES2D);
+						zi = this->GetZcoord(xyz_list, &gauss);
+	               vel[i+j*NUMVERTICES2D] = vb[i] + vsh[i]*(1.0-pow((s[i]-zi)/H[i], (n[i]+1)));
+					}
+            }
+				xDelete<IssmDouble>(xyz_list);
+
+            /*Add to the bottom side of the element*/
+            penta->AddInput(targetVel_enum,&vel[0],P1xP4Enum);
+            if (penta->IsOnSurface()) break;
+            penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
+         }
+      }
+      else _error_("not implemented yet");
+   }
+
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26530)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26531)
@@ -2604,4 +2604,5 @@
 				input->Serve(numindices,&indices[0]);
 				break;
+			case P1xP4Enum:
 			case P1DGEnum:
 			case P1bubbleEnum: