Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26132)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26133)
@@ -2886,5 +2886,5 @@
 
 	/*Intermediaries */
-	IssmDouble  thickness,Jdet,slope[2];
+	IssmDouble  thickness,Jdet,slope[2],n;
 	IssmDouble* xyz_list = NULL;
 
@@ -2900,4 +2900,5 @@
 	Input*     thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
 	Input*     surface_input  =element->GetInput(SurfaceEnum);   _assert_(surface_input);
+	Input*     n_input        =element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
 	IssmDouble rhog = element->FindParam(MaterialsRhoIceEnum)*element->FindParam(ConstantsGEnum);
 
@@ -2910,8 +2911,11 @@
 		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++){
-			pe->values[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i];
-			pe->values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i];
+			pe->values[i+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; //F1
+			pe->values[i+3]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F2
+			pe->values[i+6]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]; //F3
+			pe->values[i+9]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F4
 		}
 	}