Index: ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18223)
+++ ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18224)
@@ -47,8 +47,11 @@
 
 	/*Intermediaries */
 	int         num_responses,i;
-	IssmDouble  hobs,hu2,weight,NUMx,NUMy,DEN,Jdet;
-	IssmDouble  vx,vy,vbar2,nux,nuy,phi,dphi[2];
+	IssmDouble  hobs,hu2,weight,Jdet;
+	IssmDouble  NUMxH2,NUMyH2,DENH2;
+	IssmDouble  NUMxUbar,NUMyUbar,DENUbar;
+	IssmDouble  vxobs,vyobs,vxobsbar,vyobsbar,vbarobs2,vbarobs;
+	IssmDouble  nux,nuy,phi,dphi[2];
 	int        *responses = NULL;
 	IssmDouble *xyz_list  = NULL;
 
@@ -67,8 +70,8 @@
 	Input* thicknessobs_input = element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
 	Input* weights_input      = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
 	Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
-	Input* vx_input        = element->GetInput(BalancethicknessVxObsEnum); _assert_(vx_input);
-	Input* vy_input        = element->GetInput(BalancethicknessVyObsEnum); _assert_(vy_input);
+	Input* vxobs_input     = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
+	Input* vyobs_input     = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
 	Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
 	Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
 
@@ -81,29 +84,43 @@
 		element->NodalFunctions(basis,gauss);
 		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
-		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
+		vxobs_input->GetInputValue(&vxobs,gauss);
+		vyobs_input->GetInputValue(&vyobs,gauss);
 		nux_input->GetInputValue(&nux,gauss);
 		nuy_input->GetInputValue(&nuy,gauss);
 		potential_input->GetInputValue(&phi,gauss);
 		potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
 		thicknessobs_input->GetInputValue(&hobs,gauss);
 
-		vbar2 = (nux*nux*vx*vx + nuy*nuy*vy*vy);
-		hu2 = hobs*hobs*vbar2;
+		vxobsbar = nux*vxobs;
+		vyobsbar = nuy*vyobs;
 
-		NUMx = 2.*dbasis[0]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
-		NUMy = 2.*dbasis[1]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
-		DEN = vbar2*vbar2+1.e-14;
+		vbarobs2 = (nux*nux*vxobs*vxobs + nuy*nuy*vyobs*vyobs);
+		vbarobs  = sqrt(vbarobs2);
+		hu2 = hobs*hobs*vbarobs2;
 
+		/*H^2 - Hobs^2*/
+		NUMxH2 = 2.*dbasis[0]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
+		NUMyH2 = 2.*dbasis[1]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
+		DENH2 = vbarobs2*vbarobs2+1.e-14;
+
+		/*Ubar-Ubar_obs*/
+		NUMxUbar = (vyobsbar*dphi[0]*dphi[1] - vxobsbar*dphi[1]*dphi[1])*vbarobs*dbasis[0];
+		NUMyUbar = (vyobsbar*dphi[0]*dphi[0] - vxobsbar*dphi[0]*dphi[1])*vbarobs*dbasis[1];
+		DENUbar  = pow(dphi[0]*dphi[0] + dphi[1]*dphi[1],3./2.)+1.e-14;
+
 		/*Loop over all requested responses*/
 		for(int resp=0;resp<num_responses;resp++){
 			weights_input->GetInputValue(&weight,gauss,responses[resp]);
 
 			switch(responses[resp]){
 				case Balancethickness2MisfitEnum:
-					//for(i=0;i<numnodes;i++) pe->values[i]+=(NUMx+NUMy)/DEN *weight*Jdet*gauss->weight;
-					for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight;
+					/*J = (H^2 - Hobs^2)^2*/
+					//for(i=0;i<numnodes;i++) pe->values[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight; NOT WORKING
+					/*J = phi^2*/
+					//for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight; OK
+					/*J = (ubar - nux*uobs)^2*/
+					for(i=0;i<numnodes;i++) pe->values[i]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight;
 					break;
 				default:
 					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");