Index: /issm/trunk/src/m/solutions/ice/VelocityExtrude.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/VelocityExtrude.m	(revision 647)
+++ /issm/trunk/src/m/solutions/ice/VelocityExtrude.m	(revision 648)
@@ -14,5 +14,5 @@
 	%And remove grids on the border.
 	%These grids are on macayeal but not on pattyn nor hutter nor dirichlet
-	grids2d=md.gridonmacayeal(1:md.numberofgrids2d) & ~(md.gridonhutter(1:md.numberofgrids2d) | md.gridonpattyn(1:md.numberofgrids2d) | md.gridondirichlet_diag(1:md.numberofgrids2d));
+	grids2d=md.gridonmacayeal(1:md.numberofgrids2d);%  & ~(md.gridonhutter(1:md.numberofgrids2d) | md.gridonpattyn(1:md.numberofgrids2d) | md.gridondirichlet_diag(1:md.numberofgrids2d));
 	
 	vx2d=u_g(1:6:6*md.numberofgrids2d).*grids2d;
@@ -20,14 +20,10 @@
 
 	%Extrude across the 3d mesh
-	vx3d=project3d(md,vx2d,'node');
-	vy3d=project3d(md,vy2d,'node');
-
-	%Remove velocity of the first layer
-	vx3d(1:md.numberofgrids2d)=0;
-	vy3d(1:md.numberofgrids2d)=0;
+	vx3d=project3d(md,vx2d,'node')+u_g(1:6:end).*double(~md.gridonmacayeal);
+	vy3d=project3d(md,vy2d,'node')+u_g(2:6:end).*double(~md.gridonmacayeal);
 
 	%Plug back into u_g
-	u_g(1:6:end)=u_g(1:6:end)+vx3d;
-	u_g(2:6:end)=u_g(2:6:end)+vy3d;
+	u_g(1:6:end)=vx3d;
+	u_g(2:6:end)=vy3d;
 
 end
Index: /issm/trunk/src/m/solutions/ice/diagnostic3d.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/diagnostic3d.m	(revision 647)
+++ /issm/trunk/src/m/solutions/ice/diagnostic3d.m	(revision 648)
@@ -39,7 +39,15 @@
 	disp(sprintf('\n%s',['computing hutter velocities...']));
 	u_g=diagnostic_core_linear(m_dhu,'diagnostic_hutter',inputs);
+
+	if fem.ismacayealpattyn,
+		gridset=fem.m_dh.gridset;
+		fem.m_dh.ys=Reducevector_g(u_g);
+	end
+
 end
 
 if fem.ismacayealpattyn,
+
+	%Get field of fem
 	m_dh=fem.m_dh;