Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26147)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26148)
@@ -2707,6 +2707,5 @@
 /*MLHO*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/
-	_error_("not implemented yet");
-
+_error_("not implemented yet");
 	/* Check if ice in element */
 	if(!element->IsIceInElement()) return NULL;
@@ -2723,5 +2722,4 @@
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHOFriction(Element* element){/*{{{*/
-	_error_("not implemented yet");
 
 	if(!element->IsOnBase() || element->IsFloating()) return NULL;
@@ -2746,5 +2744,4 @@
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHOViscous(Element* element){/*{{{*/
-	_error_("not implemented yet");
 
 	/*Intermediaries*/
@@ -2861,5 +2858,4 @@
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorMLHO(Element* element){/*{{{*/
-	_error_("not implemented yet");
 
 	/* Check if ice in element */
@@ -2895,5 +2891,4 @@
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorMLHODrivingStress(Element* element){/*{{{*/
-	_error_("not implemented yet");
 
 	/*Intermediaries */
@@ -2943,5 +2938,4 @@
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorMLHOFront(Element* element){/*{{{*/
-	_error_("not implemented yet");
 
 	/*If no front, return NULL*/
@@ -2992,6 +2986,6 @@
 		pressure = ice_pressure + water_pressure;
 		/*Vertically integrated pressure - HO type*/
-		s=max(0.,thickness+base-sealevel);
-		b=min(0.,base-sealevel); 
+		b=min(0.,base-sealevel); // 0 if the bottom of the glacier is above water level
+		s=thickness+b; // ice surface regardless of whether the top of the glacier is above water level or not
 		water_pressure_sh = gravity*rho_water*(-b*b/2 + pow(-b,n+2)*( -s/(n+2) -b/(n+3) )/pow(thickness,n+1));
 		ice_pressure_sh   = gravity*rho_ice*thickness*thickness*(n+1)/(2*(n+3));
@@ -3017,6 +3011,4 @@
 }/*}}}*/
 void           StressbalanceAnalysis::InputUpdateFromSolutionMLHO(IssmDouble* solution,Element* element){/*{{{*/
-
-	_error_("not implemented yet");
 
 	int         i,dim,domaintype;
@@ -3056,7 +3048,8 @@
 			break;
 		case Domain3DEnum:
-			if(!element->IsOnBase()){xDelete<IssmDouble>(xyz_list); return;}
-			basalelement=element->SpawnBasalElement();
-			break;
+			_error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+			//if(!element->IsOnBase()){xDelete<IssmDouble>(xyz_list); return;}
+			//basalelement=element->SpawnBasalElement();
+			//break;
 		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 	}
@@ -3064,14 +3057,17 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = basalelement->GetNumberOfNodes();
-	int numdof   = numnodes*2;
+	int numdof   = numnodes*2*2; //4 DOFs per node
 
 	/*Fetch dof list and allocate solution vectors*/
-	basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
+	basalelement->GetDofListLocal(&doflist,MLHOApproximationEnum,GsetEnum); //itapopo check this
+	//basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
 	IssmDouble* values    = xNew<IssmDouble>(numdof);
 	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
 	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+   IssmDouble* vshx      = xNew<IssmDouble>(numnodes); 
+   IssmDouble* vshy      = xNew<IssmDouble>(numnodes);
 	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
 	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
-
+	
 	/*Use the dof list to index into the solution vector: */
 	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
@@ -3079,21 +3075,33 @@
 	/*Transform solution in Cartesian Space*/
 	basalelement->TransformSolutionCoord(&values[0],XYEnum);
-	basalelement->FindParam(&domaintype,DomainTypeEnum);
-
-	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numnodes;i++){
-		vx[i]=values[i*2+0];
-		vy[i]=values[i*2+1];
-
-		/*Check solution*/
-		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
-		if(xIsInf<IssmDouble>(vx[i])) _error_("Inf found in solution vector");
-		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
-		if(xIsInf<IssmDouble>(vy[i])) _error_("Inf found in solution vector");
+	//basalelement->FindParam(&domaintype,DomainTypeEnum);
+
+   /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+   for(i=0;i<numnodes;i++){ //numnodes of the 2D mesh in which the MLHO is written
+      vx[i]  =values[i+0*numnodes]; //basal vx
+      vshx[i]=values[i+1*numnodes];
+		if(xIsNan<IssmDouble>(vx[i]))		_error_("NaN found in solution vector");
+      if(xIsInf<IssmDouble>(vx[i]))		_error_("Inf found in solution vector");
+		if(xIsNan<IssmDouble>(vshx[i]))	_error_("NaN found in solution vector");
+      if(xIsInf<IssmDouble>(vshx[i]))	_error_("Inf found in solution vector");
+      //if(dim==3){
+         vy[i] =values[i+2*numnodes];
+         vshy[i]=values[i+3*numnodes];
+         if(xIsNan<IssmDouble>(vy[i]))		_error_("NaN found in solution vector");
+         if(xIsInf<IssmDouble>(vy[i]))		_error_("Inf found in solution vector");
+         if(xIsNan<IssmDouble>(vshy[i]))	_error_("NaN found in solution vector");
+         if(xIsInf<IssmDouble>(vshy[i]))	_error_("Inf found in solution vector");
+      //}
+   }
+
+	/*Compute suface velocities vx and vy*/
+	if(domaintype==Domain2DhorizontalEnum) {
+		for(i=0;i<numnodes;i++) vx[i]=vx[i]+vshx[i];
+		for(i=0;i<numnodes;i++) vy[i]=vy[i]+vshy[i];
 	}
 
 	/*Get Vz and compute 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]);
+	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: */
@@ -3107,4 +3115,6 @@
 	xDelete<IssmDouble>(vy);
 	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(vshy);
+	xDelete<IssmDouble>(vshx);
 	xDelete<IssmDouble>(values);
 	xDelete<IssmDouble>(xyz_list);