Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26136)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26137)
@@ -2555,5 +2555,5 @@
 
 	/*Intermediaries*/
-	IssmDouble  Jdet,thickness,bed,sealevel,water_pressure,ice_pressure;
+	IssmDouble  Jdet,thickness,base,sealevel,water_pressure,ice_pressure;
 	IssmDouble  surface_under_water,base_under_water,pressure;
 	IssmDouble *xyz_list = NULL;
@@ -2583,11 +2583,11 @@
 	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
-		base_input->GetInputValue(&bed,gauss);
+		base_input->GetInputValue(&base,gauss);
 		sealevel_input->GetInputValue(&sealevel,gauss);
 		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
 		element->NodalFunctions(basis,gauss);
 
-		surface_under_water = min(0.,thickness+bed-sealevel); // 0 if the top of the glacier is above water level
-		base_under_water    = min(0.,bed-sealevel);           // 0 if the bottom of the glacier is above water level
+		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
 		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;
@@ -2949,5 +2949,6 @@
 
 	/*Intermediaries*/
-	IssmDouble  Jdet,thickness,bed,sealevel,water_pressure,ice_pressure;
+	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 *xyz_list = NULL;
@@ -2964,6 +2965,7 @@
 	/*Retrieve all inputs and parameters*/
 	Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
-	Input* base_input       = element->GetInput(BaseEnum);       _assert_(base_input);
-	Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
+	Input* base_input      = element->GetInput(BaseEnum);       _assert_(base_input);
+	Input* sealevel_input  = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
+	Input* n_input         = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
 	IssmDouble rho_water   = element->FindParam(MaterialsRhoSeawaterEnum);
 	IssmDouble rho_ice     = element->FindParam(MaterialsRhoIceEnum);
@@ -2977,18 +2979,28 @@
 	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
-		base_input->GetInputValue(&bed,gauss);
+		base_input->GetInputValue(&base,gauss);
+		n_input->GetInputValue(&n,gauss);
 		sealevel_input->GetInputValue(&sealevel,gauss);
 		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
 		element->NodalFunctions(basis,gauss);
 
-		surface_under_water = min(0.,thickness+bed-sealevel); // 0 if the top of the glacier is above water level
-		base_under_water    = min(0.,bed-sealevel);           // 0 if the bottom of the glacier is above water level
+		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
+		/*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;
 		pressure = ice_pressure + water_pressure;
+		/*Vertically integrated pressure - HO type*/
+		s=max(0.,thickness+base-sealevel);
+		b=min(0.,base-sealevel); 
+		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));
+		pressure_sh = ice_pressure_sh + water_pressure_sh;
 
 		for (int i=0;i<numnodes;i++){
-			pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
-			pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
+			pe->values[i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; // F1
+			pe->values[i+3]+= pressure_sh*Jdet*gauss->weight*normal[0]*basis[i]; // F2
+			pe->values[i+6]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; // F3
+			pe->values[i+9]+= pressure_sh*Jdet*gauss->weight*normal[1]*basis[i]; // F4
 		}
 	}