Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26204)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26205)
@@ -2732,5 +2732,5 @@
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/
 	
-	_error_("Mono Layer Higher-Order called, not fully tested. If you are sure in using it, comment this line.");
+	_error_("Mono Layer Higher-Order called, not fully tested. If you are sure of using it, comment this line.");
 	
 	/* Check if ice in element */
@@ -2914,5 +2914,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss      = element->NewGauss(5);
+	Gauss* gauss      = element->NewGauss(2);
 	Gauss* gauss_base = basalelement->NewGauss();
 	while(gauss->next()){
@@ -3192,4 +3192,5 @@
 	switch(domaintype){
 		case Domain2DhorizontalEnum:
+			// itapopo FIXME use basalelement or element only here
 			basalelement = element;
 			break;
@@ -3216,6 +3217,9 @@
 	IssmDouble* vsx       = xNew<IssmDouble>(numnodes);
 	IssmDouble* vsy       = xNew<IssmDouble>(numnodes);
+	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
 	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
 	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+	IssmDouble* n			 = xNew<IssmDouble>(numnodes);
 	
 	/*Use the dof list to index into the solution vector: */
@@ -3252,13 +3256,8 @@
 	element->AddBasalInput(VxShearEnum,vshx,element->GetElementType());
 	element->AddBasalInput(VyShearEnum,vshy,element->GetElementType());
-
-	/*Get Vz and compute vel (base)*/
-	basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
-	for(i=0;i<numnodes;i++) vel[i]=sqrt(vbx[i]*vbx[i] + vby[i]*vby[i] + vz[i]*vz[i]); 
-
+	
 	/*Add vx and vy as inputs to the tria element (base velocities): */
 	element->AddBasalInput(VxBaseEnum,vbx,element->GetElementType());
 	element->AddBasalInput(VyBaseEnum,vby,element->GetElementType());
-	element->AddBasalInput(VelEnum,vel,element->GetElementType()); //itapopo check this
 	
 	/*Add vx and vy as inputs to the tria element (surface velocities): */
@@ -3266,6 +3265,24 @@
 	element->AddBasalInput(VySurfaceEnum,vsy,element->GetElementType());
 
+	/*Compute the vertically averaged velocities on each node*/
+	basalelement->GetInputListOnNodes(&n[0],MaterialsRheologyNEnum,0.); 
+   for(i=0;i<numnodes;i++){ //numnodes of the 2D mesh in which the MLHO is written
+		vx[i]=vbx[i]+vshx[i]*(n[i]+1)/(n[i]+2);
+		vy[i]=vby[i]+vshy[i]*(n[i]+1)/(n[i]+2);
+	}
+		
+	/*Get Vz and compute vel (vertically averaged vel)*/
+	basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.); 
+	for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]); 
+
+	/*Add vx and vy as inputs to the tria element (vertically averaged velocities): */
+	element->AddBasalInput(VxEnum,vx,element->GetElementType());
+	element->AddBasalInput(VyEnum,vy,element->GetElementType());
+	element->AddBasalInput(VelEnum,vel,element->GetElementType()); 
+
 	/*Free ressources:*/
 	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(vy);
 	xDelete<IssmDouble>(vz);
 	xDelete<IssmDouble>(vby);
@@ -3277,4 +3294,5 @@
 	xDelete<IssmDouble>(values);
 	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(n);
 	xDelete<int>(doflist);
 	if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};