Index: ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18392)
+++ ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18393)
@@ -233,8 +233,8 @@
 
 	/*Intermediaries*/
 	IssmDouble Jdet,weight;
-	IssmDouble vxobs,vyobs;
-	IssmDouble nux,nuy,dphi[2];
+	IssmDouble vxobs,vyobs,thicknessobs;
+	IssmDouble nux,nuy,dphi[2],dphinorm2,nuvobs2;
 	IssmDouble *xyz_list= NULL;
 
 	/*Fetch number of vertices for this finite element*/
@@ -254,6 +254,7 @@
 	Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
 	Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
 	Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
+	Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
 
 	Gauss* gauss=element->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
@@ -267,10 +268,19 @@
 		nux_input->GetInputValue(&nux,gauss);
 		nuy_input->GetInputValue(&nuy,gauss);
 		potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+		thicknessobs_input->GetInputValue(&thicknessobs,gauss);
 
+		dphinorm2 = dphi[0]*dphi[0] + dphi[1]*dphi[1];
+		nuvobs2   = nux*vxobs*nux*vxobs + nuy*vyobs*nuy*vyobs;
+
 		/*Build gradient vector (actually -dJ/da): */
 		for(int i=0;i<numvertices;i++){
-			ge[i]+= - weight*Jdet*gauss->weight*(-vxobs)*dphi[1]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i];
+			//ge[i]+= - weight*Jdet*gauss->weight*(-vxobs)*dphi[1]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i]; //THIS IS FOR Q//V
+			ge[i]+= - weight*Jdet*gauss->weight*basis[i]*(
+						-2.*nux*vxobs*vxobs*(
+							dphinorm2*dphinorm2 - thicknessobs*thicknessobs*dphinorm2*nuvobs2
+							)/ ( (pow(nux,6)+pow(nuy,6))*pow(vxobs,6) + 3*(nux*vxobs*nuy*vyobs)*(nux*vxobs*nuy*vyobs)*nuvobs2)
+						);
 			_assert_(!xIsNan<IssmDouble>(ge[i]));
 		}
 	}
@@ -290,8 +300,8 @@
 
 	/*Intermediaries*/
 	IssmDouble Jdet,weight;
-	IssmDouble vxobs,vyobs;
-	IssmDouble nux,nuy,dphi[2];
+	IssmDouble vxobs,vyobs,thicknessobs;
+	IssmDouble nux,nuy,dphi[2],dphinorm2,nuvobs2;
 	IssmDouble *xyz_list= NULL;
 
 	/*Fetch number of vertices for this finite element*/
@@ -305,12 +315,13 @@
 	/*Retrieve all inputs we will be needing: */
 	element->GetVerticesCoordinates(&xyz_list);
 	element->GradientIndexing(&vertexpidlist[0],control_index);
-	Input* weights_input   = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_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);
-	Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
+	Input* weights_input      = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_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);
+	Input* potential_input    = element->GetInput(PotentialEnum);             _assert_(potential_input);
+	Input* thicknessobs_input = element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
 
 
 	Gauss* gauss=element->NewGauss(2);
@@ -325,10 +336,20 @@
 		nux_input->GetInputValue(&nux,gauss);
 		nuy_input->GetInputValue(&nuy,gauss);
 		potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+		thicknessobs_input->GetInputValue(&thicknessobs,gauss);
 
+		dphinorm2 = dphi[0]*dphi[0] + dphi[1]*dphi[1];
+		nuvobs2   = nux*vxobs*nux*vxobs + nuy*vyobs*nuy*vyobs;
+
+
 		/*Build gradient vector (actually -dJ/da): */
 		for(int i=0;i<numvertices;i++){
-			ge[i]+= - weight*Jdet*gauss->weight*(+vyobs)*dphi[0]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i];
+			//ge[i]+= - weight*Jdet*gauss->weight*(+vyobs)*dphi[0]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i];
+			ge[i]+= - weight*Jdet*gauss->weight*basis[i]*(
+						-2.*nuy*vyobs*vyobs*(
+							dphinorm2*dphinorm2 - thicknessobs*thicknessobs*dphinorm2*nuvobs2
+							)/ ( (pow(nux,6)+pow(nuy,6))*pow(vyobs,6) + 3*(nux*vxobs*nuy*vyobs)*(nux*vxobs*nuy*vyobs)*nuvobs2)
+						);
 			_assert_(!xIsNan<IssmDouble>(ge[i]));
 		}
 	}
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18392)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18393)
@@ -1354,7 +1354,7 @@
 			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 = grad phi^2*/
@@ -1362,7 +1362,7 @@
 			/*J = (ubar - nux*uobs)^2*/
 			//J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
 			/*J = 1/2 (vbar ^ gard(phi))^2*/
-			J += 0.5*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*weight*Jdet*gauss->weight;
+			//J += 0.5*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*weight*Jdet*gauss->weight;
 		}
 
 		/*clean up and Return: */
