Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 17347)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 17348)
@@ -257,5 +257,6 @@
 
 	/*Intermediaries */
-	IssmDouble ub,vb,slope2,drag,thickness,connectivity;
+	int        frictionlaw = 1;
+	IssmDouble ub,vb,slope2,drag,thickness,surface,connectivity;
 	IssmDouble slope[2];
 
@@ -274,4 +275,5 @@
 	Input* slopey_input    = element->GetInput(SurfaceSlopeYEnum);        _assert_(slopey_input);
 	Input* thickness_input = element->GetInput(ThicknessEnum);            _assert_(thickness_input);
+	Input* surface_input   = element->GetInput(SurfaceEnum);              _assert_(surface_input);
 	Input* drag_input      = element->GetInput(FrictionCoefficientEnum);  _assert_(drag_input);
 
@@ -283,4 +285,5 @@
 
 		thickness_input->GetInputValue(&thickness,gauss);
+		surface_input->GetInputValue(&surface,gauss);
 		drag_input->GetInputValue(&drag,gauss);
 		slopex_input->GetInputValue(&slope[0],gauss);
@@ -288,13 +291,29 @@
 		slope2=slope[0]*slope[0]+slope[1]*slope[1];
 
-		/*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
-		ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
-		vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
-		///*Ritz et al. 1996*/
-		//ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
-		//vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
-		///*Rutt et al. 2009*/
-		//ub=-drag*rho_ice*gravity*thickness*slope[0];
-		//vb=-drag*rho_ice*gravity*thickness*slope[1];
+		switch(frictionlaw){
+			case 1:
+				/*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
+				ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
+				vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
+				break;
+			case 2:
+				/*Ritz et al. 1996*/
+				ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
+				vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
+				break;
+			case 3:
+				/*Rutt et al. 2009*/
+				ub=-drag*rho_ice*gravity*thickness*slope[0];
+				vb=-drag*rho_ice*gravity*thickness*slope[1];
+				break;
+			case 4:
+				/*Henning Akesson*/
+				drag = -4e-15 * surface + 8.6e-12;
+				ub=-drag*rho_ice*gravity*thickness*slope[0];
+				vb=-drag*rho_ice*gravity*thickness*slope[1];
+				break;
+			default:
+				_error_("Not supported yet");
+		}
 
 		pe->values[2*iv+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/connectivity;
@@ -309,4 +328,5 @@
 
 	/*Intermediaries */
+	int         frictionlaw = 1;
 	int         nodeup,nodedown,numsegments;
 	IssmDouble  ub,vb,slope2,drag,surface,thickness,constant_part,z,Jdet;
@@ -374,13 +394,29 @@
 			drag_input->GetInputValue(&drag,gauss);
 
-			/*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
-			ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
-			vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
-			///*Ritz et al. 1996*/
-			//ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
-			//vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
-			///*Rutt et al. 2009*/
-			//ub=-drag*rho_ice*gravity*thickness*slope[0];
-			//vb=-drag*rho_ice*gravity*thickness*slope[1];
+			switch(frictionlaw){
+				case 1:
+					/*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
+					ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
+					vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
+					break;
+				case 2:
+					/*Ritz et al. 1996*/
+					ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
+					vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
+					break;
+				case 3:
+					/*Rutt et al. 2009*/
+					ub=-drag*rho_ice*gravity*thickness*slope[0];
+					vb=-drag*rho_ice*gravity*thickness*slope[1];
+					break;
+				case 4:
+					/*Henning Akesson*/
+					drag = -4e-15 * surface + 8.6e-12;
+					ub=-drag*rho_ice*gravity*thickness*slope[0];
+					vb=-drag*rho_ice*gravity*thickness*slope[1];
+					break;
+				default:
+					_error_("Not supported yet");
+			}
 
 			pe->values[2*nodedown+0]+=ub/connectivity[0];