Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18389)
+++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18390)
@@ -108,28 +108,27 @@
 					/*J = (H^2 - Hobs^2)^2*/
 					//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
-					/*J = grad J ^2*/
+					/*J = grad phi ^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++){
+					//	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]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight;
+					//}
+					/*J = 1/2 (vbar ^ gard(phi))^2*/
 					for(i=0;i<numnodes;i++){
-						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]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight;
+						pe->values[i]+= weight*Jdet*gauss->weight*
+						  (nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*
+						  (nuy*vyobs*dbasis[0*numnodes+i] - nux*vxobs*dbasis[1*numnodes+i]);
 					}
+
 					break;
 				default:
@@ -168,9 +167,4 @@
 	element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
 
-	/*Check that control_type is supported*/
-	if(control_type!=BalancethicknessApparentMassbalanceEnum){
-		_error_("Control "<<EnumToStringx(control_type)<<" not supported");
-	}
-
 	/*Deal with first part (partial derivative a J with respect to k)*/
 	for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
@@ -182,4 +176,6 @@
 	switch(control_type){
 		case BalancethicknessApparentMassbalanceEnum: GradientJAdot(element,gradient,control_index); break;
+		case BalancethicknessNuxEnum: GradientJNux(element,gradient,control_index); break;
+		case BalancethicknessNuyEnum: GradientJNuy(element,gradient,control_index); break;
 		default: _error_("control type not supported yet: " << EnumToStringx(control_type));
 	}
@@ -234,4 +230,114 @@
 	delete gauss;
 }/*}}}*/
+void AdjointBalancethickness2Analysis::GradientJNux(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble Jdet,weight;
+	IssmDouble vxobs,vyobs;
+	IssmDouble nux,nuy,dphi[2];
+	IssmDouble *xyz_list= NULL;
+
+	/*Fetch number of vertices for this finite element*/
+	int numvertices = element->GetNumberOfVertices();
+
+	/*Initialize some vectors*/
+	IssmDouble* basis         = xNew<IssmDouble>(numvertices);
+	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
+	int*        vertexpidlist = xNew<int>(numvertices);
+
+	/*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);
+
+
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctionsP1(basis,gauss);
+		weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
+		vxobs_input->GetInputValue(&vxobs,gauss);
+		vyobs_input->GetInputValue(&vyobs,gauss);
+		nux_input->GetInputValue(&nux,gauss);
+		nuy_input->GetInputValue(&nuy,gauss);
+		potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+
+		/*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];
+			_assert_(!xIsNan<IssmDouble>(ge[i]));
+		}
+	}
+	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(ge);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<int>(vertexpidlist);
+	delete gauss;
+}/*}}}*/
+void AdjointBalancethickness2Analysis::GradientJNuy(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble Jdet,weight;
+	IssmDouble vxobs,vyobs;
+	IssmDouble nux,nuy,dphi[2];
+	IssmDouble *xyz_list= NULL;
+
+	/*Fetch number of vertices for this finite element*/
+	int numvertices = element->GetNumberOfVertices();
+
+	/*Initialize some vectors*/
+	IssmDouble* basis         = xNew<IssmDouble>(numvertices);
+	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
+	int*        vertexpidlist = xNew<int>(numvertices);
+
+	/*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);
+
+
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctionsP1(basis,gauss);
+		weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
+		vxobs_input->GetInputValue(&vxobs,gauss);
+		vyobs_input->GetInputValue(&vyobs,gauss);
+		nux_input->GetInputValue(&nux,gauss);
+		nuy_input->GetInputValue(&nuy,gauss);
+		potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+
+		/*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];
+			_assert_(!xIsNan<IssmDouble>(ge[i]));
+		}
+	}
+	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(ge);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<int>(vertexpidlist);
+	delete gauss;
+}/*}}}*/
 void AdjointBalancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 	element->InputUpdateFromSolutionOneDof(solution,AdjointEnum);
Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h	(revision 18389)
+++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h	(revision 18390)
@@ -29,4 +29,6 @@
 		void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
 		void GradientJAdot(Element* element,Vector<IssmDouble>* gradient,int control_index);
+		void GradientJNux(Element* element,Vector<IssmDouble>* gradient,int control_index);
+		void GradientJNuy(Element* element,Vector<IssmDouble>* gradient,int control_index);
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
 		void UpdateConstraints(FemModel* femmodel);
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 18389)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 18390)
@@ -1361,5 +1361,7 @@
 			//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;
+			/*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;
 		}