Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18392)
+++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18393)
@@ -234,6 +234,6 @@
 	/*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;
 
@@ -255,4 +255,5 @@
 	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);
@@ -268,8 +269,17 @@
 		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]));
 		}
@@ -291,6 +301,6 @@
 	/*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;
 
@@ -306,10 +316,11 @@
 	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);
 
 
@@ -326,8 +337,18 @@
 		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: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 18392)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 18393)
@@ -1355,5 +1355,5 @@
 
 			/*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
@@ -1363,5 +1363,5 @@
 			//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;
 		}