Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 18221)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 18222)
@@ -1207,4 +1207,5 @@
 
 	IssmDouble  weight,thicknessobs,thickness,potential;
+	IssmDouble  vx,vy,vxbar,vybar,vxobs,vyobs,vxbarobs,vybarobs,nux,nuy;
 	IssmDouble  Jdet;
 	IssmDouble* xyz_list = NULL;
@@ -1223,4 +1224,10 @@
 		Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
 		Input* potential_input   =element->GetInput(PotentialEnum);                          _assert_(potential_input);
+		Input* vxobs_input       =element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
+		Input* vyobs_input       =element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
+		Input* vx_input          =element->GetInput(VxEnum); _assert_(vx_input);
+		Input* vy_input          =element->GetInput(VyEnum); _assert_(vy_input);
+		Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
+		Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
 
 		/* Start  looping on the number of gaussian points: */
@@ -1238,8 +1245,20 @@
 			thicknessobs_input->GetInputValue(&thicknessobs,gauss);
 			potential_input->GetInputValue(&potential,gauss);
-			gauss->GaussPoint(ig);
-
+			vxobs_input->GetInputValue(&vxobs,gauss);
+			vyobs_input->GetInputValue(&vyobs,gauss);
+			vx_input->GetInputValue(&vxbar,gauss);
+			vy_input->GetInputValue(&vybar,gauss);
+			nux_input->GetInputValue(&nux,gauss);
+			nuy_input->GetInputValue(&nuy,gauss);
+
+			vxbarobs = nux*vxobs;
+			vybarobs = nuy*vyobs;
+
+			/*J = (H^2 - Hobs^2)^2*/
 			//J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
-			J +=.5*potential*potential*weight*Jdet*gauss->weight;
+			/*J = phi^2*/
+			//J +=.5*potential*potential*weight*Jdet*gauss->weight; OK
+			/*J = (ubar - nux*uobs)^2*/
+			J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
 		}