Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26289)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26290)
@@ -193,13 +193,11 @@
 					IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,1);
 				}
 			}
-			else{//itapopo FIXME only base velocities are constrained here
+			else{//MLHO itapopo FIXME base and shear velocities are constrained
 				IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0);
-				//IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1);
-				if(iomodel->domaintype!=Domain2DverticalEnum){
-					IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2);
-					//IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3);
-				}
+				IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1);
+				IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2);
+				IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3);
 			}
 		}
 
@@ -3065,7 +3063,7 @@
 		thickness_input->GetInputValue(&thickness,gauss);
 		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
 		n_input->GetInputValue(&n,gauss);
-
+		
 		for(int i=0;i<numnodes;i++){// per node: vx (basal vx), vshx, vy (basal vy), vshy
 			pe->values[i*4+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; //F1
 			pe->values[i*4+1]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F2
@@ -3090,8 +3088,7 @@
 
 	/*Intermediaries*/
 	IssmDouble  Jdet,thickness,base,sealevel,water_pressure,ice_pressure;
-	IssmDouble  water_pressure_sh,ice_pressure_sh,pressure_sh,n,s,b;
-	IssmDouble  surface_under_water,base_under_water,pressure;
+	IssmDouble  water_pressure_sh,ice_pressure_sh,pressure_sh,pressure,n,s,b;
 	IssmDouble *xyz_list = NULL;
 	IssmDouble *xyz_list_front = NULL;
 	IssmDouble  normal[2];
@@ -3124,17 +3121,15 @@
 		sealevel_input->GetInputValue(&sealevel,gauss);
 		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
 		element->NodalFunctions(basis,gauss);
-
-		surface_under_water = min(0.,thickness+base-sealevel); // 0 if the top of the glacier is above water level
-		base_under_water    = min(0.,base-sealevel);           // 0 if the bottom of the glacier is above water level
+		
+		b = base-sealevel; // ice base shifted by the sea level
+		s = thickness+b;   // ice surface shifted by the sea level
 		/*Vertically integrated pressure - SSA type*/
-		water_pressure = (1.0/2.0)*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
-		ice_pressure   = (1.0/2.0)*gravity*rho_ice*thickness*thickness;
+		water_pressure = -(1.0/2.0)*gravity*rho_water*(b*b);
+		ice_pressure   =  (1.0/2.0)*gravity*rho_ice*thickness*thickness;
 		pressure = ice_pressure + water_pressure;
 		/*Vertically integrated pressure - HO type*/
-		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)/pow(thickness,n+1))*(-s/(n+2)-b/(n+3)) );
+		water_pressure_sh = gravity*rho_water*( -b*b/2 - (s*thickness/(n+2))*(1-pow(s/thickness,n+2)) + (thickness*thickness/(n+3))*(1-pow(s/thickness,n+3)) );
 		ice_pressure_sh   = gravity*rho_ice*thickness*thickness*(n+1)/(2*(n+3));
 		pressure_sh = ice_pressure_sh + water_pressure_sh;