Index: ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18368)
+++ ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18369)
@@ -99,16 +99,6 @@
 		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]);
@@ -116,11 +106,26 @@
 			switch(responses[resp]){
 				case Balancethickness2MisfitEnum:
 					/*J = (H^2 - Hobs^2)^2*/
-					//for(i=0;i<numnodes;i++) pe->values[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight; NOT WORKING
+					for(i=0;i<numnodes;i++){
+
+						/*H^2 - Hobs^2*/
+						NUMxH2 = 2.*dbasis[0*numnodes+i]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
+						NUMyH2 = 2.*dbasis[1*numnodes+i]*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*numnodes+i];
+						NUMyUbar = (vyobsbar*dphi[0]*dphi[0] - vxobsbar*dphi[0]*dphi[1])*vbarobs*dbasis[1*numnodes+i];
+						DENUbar  = pow(dphi[0]*dphi[0] + dphi[1]*dphi[1],3./2.)+1.e-14;
+
+						pe->values[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight;
+					}
 					/*J = phi^2*/
-					//for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight; OK
+					//for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight; //OK
+					/*J = grad J ^2*/
+					//for(i=0;i<numnodes;i++) pe->values[i]+= (dphi[0]*dbasis[0*numnodes+i] + dphi[1]*dbasis[1*numnodes+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;
+					//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");
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18368)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18369)
@@ -1303,7 +1303,7 @@
 	IssmDouble J=0.;
 	IssmDouble J_sum;
 
-	IssmDouble  weight,thicknessobs,thickness,potential;
+	IssmDouble  weight,thicknessobs,thickness,potential,dpotential[2];
 	IssmDouble  vx,vy,vxbar,vybar,vxobs,vyobs,vxbarobs,vybarobs,nux,nuy;
 	IssmDouble  Jdet;
 	IssmDouble* xyz_list = NULL;
@@ -1342,6 +1342,7 @@
 			thickness_input->GetInputValue(&thickness,gauss);
 			thicknessobs_input->GetInputValue(&thicknessobs,gauss);
 			potential_input->GetInputValue(&potential,gauss);
+			potential_input->GetInputDerivativeValue(&dpotential[0],xyz_list,gauss);
 			vxobs_input->GetInputValue(&vxobs,gauss);
 			vyobs_input->GetInputValue(&vyobs,gauss);
 			vx_input->GetInputValue(&vxbar,gauss);
@@ -1353,11 +1354,13 @@
 			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 +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
 			/*J = phi^2*/
-			//J +=.5*potential*potential*weight*Jdet*gauss->weight; OK
+			//J +=.5*potential*potential*weight*Jdet*gauss->weight;// OK
+			/*J = grad phi^2*/
+			//J +=.5*(dpotential[0]*dpotential[0] + dpotential[1]*dpotential[1])*weight*Jdet*gauss->weight;
 			/*J = (ubar - nux*uobs)^2*/
-			J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
+			//J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
 		}
 
 		/*clean up and Return: */
