Index: ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp	(revision 18402)
+++ ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp	(revision 18403)
@@ -198,8 +198,7 @@
 void AdjointBalancethicknessAnalysis::GradientJVx(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 
 	/*Intermediaries*/
-	IssmDouble Jdet,weight;
-	IssmDouble thickness,Dlambda[3],dp[3];
+	IssmDouble thickness,Jdet,Dlambda[3],dp[3];
 	IssmDouble *xyz_list= NULL;
 
 	/*Fetch number of vertices for this finite element*/
@@ -246,8 +245,7 @@
 void AdjointBalancethicknessAnalysis::GradientJVy(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 
 	/*Intermediaries*/
-	IssmDouble Jdet,weight;
-	IssmDouble thickness,Dlambda[3],dp[3];
+	IssmDouble thickness,Jdet,Dlambda[3],dp[3];
 	IssmDouble *xyz_list= NULL;
 
 	/*Fetch number of vertices for this finite element*/
Index: ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18402)
+++ ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18403)
@@ -51,7 +51,7 @@
 	IssmDouble  NUMxH2,NUMyH2,DENH2;
 	IssmDouble  NUMxUbar,NUMyUbar,DENUbar;
 	IssmDouble  vxobs,vyobs,vxobsbar,vyobsbar,vbarobs2,vbarobs;
-	IssmDouble  nux,nuy,phi,dphi[2];
+	IssmDouble  nu,phi,dphi[2];
 	int        *responses = NULL;
 	IssmDouble *xyz_list  = NULL;
 
@@ -72,8 +72,7 @@
 	Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_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* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_input);
 
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(2);
@@ -86,16 +85,15 @@
 
 		vxobs_input->GetInputValue(&vxobs,gauss);
 		vyobs_input->GetInputValue(&vyobs,gauss);
-		nux_input->GetInputValue(&nux,gauss);
-		nuy_input->GetInputValue(&nuy,gauss);
+		nu_input->GetInputValue(&nu,gauss);
 		potential_input->GetInputValue(&phi,gauss);
 		potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
 		thicknessobs_input->GetInputValue(&hobs,gauss);
 
-		vxobsbar = nux*vxobs;
-		vyobsbar = nuy*vyobs;
+		vxobsbar = nu*vxobs;
+		vyobsbar = nu*vyobs;
 
-		vbarobs2 = (nux*nux*vxobs*vxobs + nuy*nuy*vyobs*vyobs);
+		vbarobs2 = (nu*nu*vxobs*vxobs + nu*nu*vyobs*vyobs);
 		vbarobs  = sqrt(vbarobs2);
 		hu2 = hobs*hobs*vbarobs2;
 
@@ -104,32 +102,36 @@
 			weights_input->GetInputValue(&weight,gauss,responses[resp]);
 
 			switch(responses[resp]){
+				case IrrotationalH2MisfitEnum:
+					/*J = (H^2 - Hobs^2)^2*/
+					for(i=0;i<numnodes;i++){
+						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;
+						pe->values[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight;
+					}
+					break;
+				case IrrotationalDirectionMisfitEnum:
+					/*J = 1/2 (vbar ^ gard(phi))^2*/
+					for(i=0;i<numnodes;i++){
+						pe->values[i]+= weight*Jdet*gauss->weight*
+						  nu*nu*(vyobs*dphi[0] - vxobs*dphi[1])*
+						  (vyobs*dbasis[0*numnodes+i] - vxobs*dbasis[1*numnodes+i]);
+					}
+					break;
 				case Balancethickness2MisfitEnum:
-					/*J = (H^2 - Hobs^2)^2*/
-					//for(i=0;i<numnodes;i++){
-					//	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;
-					//	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 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*/
+					/*J = (ubar - nu*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++){
-						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]);
-					}
-
+					_error_("Not implemented yet");
 					break;
 				default:
 					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
@@ -168,15 +170,37 @@
 
 	/*Deal with first part (partial derivative a J with respect to k)*/
 	for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
-		case Balancethickness2MisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
+		case Balancethickness2MisfitEnum:
+			/*Nothing, \partial J/\partial k = 0*/
+			break;
+		case IrrotationalH2MisfitEnum:
+			switch(control_type){
+				case BalancethicknessApparentMassbalanceEnum:
+					/*Nothing, \partial J/\partial k = 0*/
+					break;
+				case BalancethicknessNuEnum:
+					GradientJ1Nu(element,gradient,control_index);
+					break;
+				default: _error_("control type not supported yet: " << EnumToStringx(control_type));
+			} 
+			break;
+		case IrrotationalDirectionMisfitEnum:
+			switch(control_type){
+				case BalancethicknessApparentMassbalanceEnum:
+					/*Nothing, \partial J/\partial k = 0*/
+					break;
+				case BalancethicknessNuEnum:
+					GradientJ2Nu(element,gradient,control_index); break; //Might need to be renamed ?
+				default: _error_("control type not supported yet: " << EnumToStringx(control_type));
+			} 
+			break;
 		default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
 	}
 
 	/*Deal with second term*/
 	switch(control_type){
 		case BalancethicknessApparentMassbalanceEnum: GradientJAdot(element,gradient,control_index); break;
-		case BalancethicknessNuxEnum: GradientJNux(element,gradient,control_index); break;//should be in first part....
-		case BalancethicknessNuyEnum: GradientJNuy(element,gradient,control_index); break;//should be in first part....
+		case BalancethicknessNuEnum: /*Nothing: state equations do not depend on k*/ break;
 		default: _error_("control type not supported yet: " << EnumToStringx(control_type));
 	}
 
@@ -187,8 +211,7 @@
 void AdjointBalancethickness2Analysis::GradientJAdot(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 
 	/*Intermediaries*/
-	IssmDouble Jdet,weight;
-	IssmDouble lambda; 
+	IssmDouble lambda,Jdet; 
 	IssmDouble *xyz_list= NULL;
 
 	/*Fetch number of vertices for this finite element*/
@@ -202,8 +225,7 @@
 	/*Retrieve all inputs we will be needing: */
 	element->GetVerticesCoordinates(&xyz_list);
 	element->GradientIndexing(&vertexpidlist[0],control_index);
-	Input* adjoint_input = element->GetInput(AdjointEnum);                            _assert_(adjoint_input);
-	Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+	Input* adjoint_input = element->GetInput(AdjointEnum); _assert_(adjoint_input);
 
 	Gauss* gauss=element->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
@@ -211,12 +233,11 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctionsP1(basis,gauss);
-		weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
 		adjoint_input->GetInputValue(&lambda,gauss);
 
 		/*Build gradient vector (actually -dJ/da): */
 		for(int i=0;i<numvertices;i++){
-			ge[i]+= - weight*Jdet*gauss->weight*basis[i]*lambda;
+			ge[i]+= - Jdet*gauss->weight*basis[i]*lambda;
 			_assert_(!xIsNan<IssmDouble>(ge[i]));
 		}
 	}
@@ -229,12 +250,11 @@
 	xDelete<int>(vertexpidlist);
 	delete gauss;
 }/*}}}*/
-void AdjointBalancethickness2Analysis::GradientJNux(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+void AdjointBalancethickness2Analysis::GradientJ1Nu(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 
 	/*Intermediaries*/
-	IssmDouble Jdet,weight;
-	IssmDouble vxobs,vyobs,thicknessobs;
-	IssmDouble nux,nuy,dphi[2],dphinorm2,nuvobs2;
+	IssmDouble vxobs,vyobs,thicknessobs,weight;
+	IssmDouble nu,Jdet,dphi[2],dphinorm2,nuvobs2;
 	IssmDouble *xyz_list= NULL;
 
 	/*Fetch number of vertices for this finite element*/
@@ -251,8 +271,7 @@
 	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* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_input);
 	Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
 	Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
 
@@ -262,24 +281,20 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctionsP1(basis,gauss);
-		weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
+		weights_input->GetInputValue(&weight,gauss,IrrotationalH2MisfitEnum);
 		vxobs_input->GetInputValue(&vxobs,gauss);
 		vyobs_input->GetInputValue(&vyobs,gauss);
-		nux_input->GetInputValue(&nux,gauss);
-		nuy_input->GetInputValue(&nuy,gauss);
+		nu_input->GetInputValue(&nu,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;
+		nuvobs2   = nu*nu*(vxobs*vxobs + vyobs*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]; //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)
+			ge[i]+= -weight*Jdet*gauss->weight*basis[i]*(
+						-2./pow(nu,3) * dphinorm2/(vxobs*vxobs + vyobs*vyobs) * (dphinorm2/nuvobs2 - thicknessobs*thicknessobs) 
 						);
 			_assert_(!xIsNan<IssmDouble>(ge[i]));
 		}
@@ -293,12 +308,11 @@
 	xDelete<int>(vertexpidlist);
 	delete gauss;
 }/*}}}*/
-void AdjointBalancethickness2Analysis::GradientJNuy(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+void AdjointBalancethickness2Analysis::GradientJ2Nu(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 
 	/*Intermediaries*/
-	IssmDouble Jdet,weight;
-	IssmDouble vxobs,vyobs,thicknessobs;
-	IssmDouble nux,nuy,dphi[2],dphinorm2,nuvobs2;
+	IssmDouble vxobs,vyobs,thicknessobs,weight;
+	IssmDouble nu,Jdet,dphi[2],dphinorm2,nuvobs2;
 	IssmDouble *xyz_list= NULL;
 
 	/*Fetch number of vertices for this finite element*/
@@ -312,41 +326,32 @@
 	/*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* thicknessobs_input = element->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_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* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_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++){
 		gauss->GaussPoint(ig);
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctionsP1(basis,gauss);
-		weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
+		weights_input->GetInputValue(&weight,gauss,IrrotationalDirectionMisfitEnum);
 		vxobs_input->GetInputValue(&vxobs,gauss);
 		vyobs_input->GetInputValue(&vyobs,gauss);
-		nux_input->GetInputValue(&nux,gauss);
-		nuy_input->GetInputValue(&nuy,gauss);
+		nu_input->GetInputValue(&nu,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;
+		nuvobs2   = nu*vxobs*nu*vxobs + nu*vyobs*nu*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*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)
-						);
+			ge[i]+= - weight*Jdet*gauss->weight*nu*(vyobs*dphi[0] -vxobs*dphi[1])*(vyobs*dphi[0] -vxobs*dphi[1])*basis[i];
 			_assert_(!xIsNan<IssmDouble>(ge[i]));
 		}
 	}
Index: ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h	(revision 18402)
+++ ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h	(revision 18403)
@@ -28,8 +28,8 @@
 		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		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 GradientJ1Nu(Element* element,Vector<IssmDouble>* gradient,int control_index);
+		void GradientJ2Nu(Element* element,Vector<IssmDouble>* gradient,int control_index);
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
 		void UpdateConstraints(FemModel* femmodel);
 };
Index: ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 18402)
+++ ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 18403)
@@ -26,8 +26,7 @@
 	}
 
 	iomodel->FetchDataToInput(elements,BalancethicknessApparentMassbalanceEnum);
-	iomodel->FetchDataToInput(elements,BalancethicknessNuxEnum);
-	iomodel->FetchDataToInput(elements,BalancethicknessNuyEnum);
+	iomodel->FetchDataToInput(elements,BalancethicknessNuEnum);
 	iomodel->FetchDataToInput(elements,BalancethicknessVxObsEnum);
 	iomodel->FetchDataToInput(elements,BalancethicknessVyObsEnum);
 	iomodel->FetchDataToInput(elements,BalancethicknessThicknessObsEnum);
@@ -202,7 +201,7 @@
 
 	/*Intermediaries */
 	int         Hinterpolation;
-	IssmDouble  vx,vy,vbar,nux,nuy,normdphi,dphi[2];
+	IssmDouble  vx,vy,vbar,nu,normdphi,dphi[2];
 	IssmDouble* xyz_list = NULL;
 	int       * doflist  = NULL;
 
@@ -232,8 +231,7 @@
 	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* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
-	Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
+	Input* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_input);
 
 	switch(element->GetElementType()){
 		case P1Enum: Hinterpolation = P0Enum; break;
@@ -246,11 +244,10 @@
 
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
-		nux_input->GetInputValue(&nux,gauss);
-		nuy_input->GetInputValue(&nuy,gauss);
+		nu_input->GetInputValue(&nu,gauss);
 		potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
 
-		vx = vx*nux; vy = vy*nuy;
+		vx = vx*nu; vy = vy*nu;
 		vbar = sqrt(vx*vx + vy*vy) + 1.e-10;
 		normdphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]);
 
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18402)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18403)
@@ -300,8 +300,12 @@
 	BalancethicknessSpcpotentialEnum,
 	BalancethicknessApparentMassbalanceEnum,
 	Balancethickness2MisfitEnum,
-	BalancethicknessNuxEnum,
-	BalancethicknessNuyEnum,
+	IrrotationalH2MisfitEnum,
+	IrrotationalDirectionMisfitEnum,
+	IrrotationalVelMisfitEnum,
+	IrrotationalAlongGradientNuEnum,
+	IrrotationalAcrossGradientNuEnum,
+	BalancethicknessNuEnum,
 	BalancethicknessVxObsEnum,
 	BalancethicknessVyObsEnum,
 	BalancethicknessThicknessObsEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18402)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18403)
@@ -308,8 +308,12 @@
 		case BalancethicknessSpcpotentialEnum : return "BalancethicknessSpcpotential";
 		case BalancethicknessApparentMassbalanceEnum : return "BalancethicknessApparentMassbalance";
 		case Balancethickness2MisfitEnum : return "Balancethickness2Misfit";
-		case BalancethicknessNuxEnum : return "BalancethicknessNux";
-		case BalancethicknessNuyEnum : return "BalancethicknessNuy";
+		case IrrotationalH2MisfitEnum : return "IrrotationalH2Misfit";
+		case IrrotationalDirectionMisfitEnum : return "IrrotationalDirectionMisfit";
+		case IrrotationalVelMisfitEnum : return "IrrotationalVelMisfit";
+		case IrrotationalAlongGradientNuEnum : return "IrrotationalAlongGradientNu";
+		case IrrotationalAcrossGradientNuEnum : return "IrrotationalAcrossGradientNu";
+		case BalancethicknessNuEnum : return "BalancethicknessNu";
 		case BalancethicknessVxObsEnum : return "BalancethicknessVxObs";
 		case BalancethicknessVyObsEnum : return "BalancethicknessVyObs";
 		case BalancethicknessThicknessObsEnum : return "BalancethicknessThicknessObs";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18402)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18403)
@@ -314,8 +314,12 @@
 	      else if (strcmp(name,"BalancethicknessSpcpotential")==0) return BalancethicknessSpcpotentialEnum;
 	      else if (strcmp(name,"BalancethicknessApparentMassbalance")==0) return BalancethicknessApparentMassbalanceEnum;
 	      else if (strcmp(name,"Balancethickness2Misfit")==0) return Balancethickness2MisfitEnum;
-	      else if (strcmp(name,"BalancethicknessNux")==0) return BalancethicknessNuxEnum;
-	      else if (strcmp(name,"BalancethicknessNuy")==0) return BalancethicknessNuyEnum;
+	      else if (strcmp(name,"IrrotationalH2Misfit")==0) return IrrotationalH2MisfitEnum;
+	      else if (strcmp(name,"IrrotationalDirectionMisfit")==0) return IrrotationalDirectionMisfitEnum;
+	      else if (strcmp(name,"IrrotationalVelMisfit")==0) return IrrotationalVelMisfitEnum;
+	      else if (strcmp(name,"IrrotationalAlongGradientNu")==0) return IrrotationalAlongGradientNuEnum;
+	      else if (strcmp(name,"IrrotationalAcrossGradientNu")==0) return IrrotationalAcrossGradientNuEnum;
+	      else if (strcmp(name,"BalancethicknessNu")==0) return BalancethicknessNuEnum;
 	      else if (strcmp(name,"BalancethicknessVxObs")==0) return BalancethicknessVxObsEnum;
 	      else if (strcmp(name,"BalancethicknessVyObs")==0) return BalancethicknessVyObsEnum;
 	      else if (strcmp(name,"BalancethicknessThicknessObs")==0) return BalancethicknessThicknessObsEnum;
@@ -378,14 +382,14 @@
 	      else if (strcmp(name,"FlaimSolution")==0) return FlaimSolutionEnum;
 	      else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
 	      else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
-	      else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
+         else stage=4;
+   }
+   if(stage==4){
+	      if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
 	      else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
 	      else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
 	      else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
-         else stage=4;
-   }
-   if(stage==4){
-	      if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
+	      else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
 	      else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
 	      else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
 	      else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum;
@@ -501,14 +505,14 @@
 	      else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
 	      else if (strcmp(name,"Converged")==0) return ConvergedEnum;
 	      else if (strcmp(name,"Fill")==0) return FillEnum;
-	      else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
+         else stage=5;
+   }
+   if(stage==5){
+	      if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
 	      else if (strcmp(name,"Friction")==0) return FrictionEnum;
 	      else if (strcmp(name,"Internal")==0) return InternalEnum;
 	      else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
-         else stage=5;
-   }
-   if(stage==5){
-	      if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
+	      else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
 	      else if (strcmp(name,"Misfit")==0) return MisfitEnum;
 	      else if (strcmp(name,"Pressure")==0) return PressureEnum;
 	      else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
@@ -624,14 +628,14 @@
 	      else if (strcmp(name,"MisfitModelEnum")==0) return MisfitModelEnumEnum;
 	      else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
 	      else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
-	      else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
+         else stage=6;
+   }
+   if(stage==6){
+	      if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
 	      else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
 	      else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
 	      else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
-         else stage=6;
-   }
-   if(stage==6){
-	      if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
+	      else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
 	      else if (strcmp(name,"MinVel")==0) return MinVelEnum;
 	      else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
 	      else if (strcmp(name,"MinVx")==0) return MinVxEnum;
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp	(revision 18402)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp	(revision 18403)
@@ -40,8 +40,7 @@
 			case VyEnum:
 			case ThicknessEnum:
 			case FrictionCoefficientEnum:
-			case BalancethicknessNuxEnum:
-			case BalancethicknessNuyEnum:
+			case BalancethicknessNuEnum:
 			case BalancethicknessApparentMassbalanceEnum:
 				iomodel->FetchData(1,control); 
 				break;
Index: ../trunk-jpl/src/c/cores/controlvalidation_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/controlvalidation_core.cpp	(revision 18402)
+++ ../trunk-jpl/src/c/cores/controlvalidation_core.cpp	(revision 18403)
@@ -11,8 +11,10 @@
 void controlvalidation_core(FemModel* femmodel){
 
 	int         solution_type,n;
+	int         num_responses;
 	IssmDouble  j0,j;
 	IssmDouble  Ialpha,exponent,alpha;
+	IssmDouble* jlist = NULL;
 	IssmDouble *G = NULL;
 	IssmDouble *X = NULL;
 	IssmDouble *X0= NULL;
@@ -24,6 +26,7 @@
 	/*Recover parameters used throughout the solution*/
 	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
 	femmodel->parameters->SetParam(false,SaveResultsEnum);
+	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
 
 	/*Get initial guess*/
 	Vector<IssmDouble> *Xpetsc = NULL;
@@ -44,8 +47,16 @@
 	solutioncore(femmodel);
 	if(VerboseControl()) _printf0_("   Compute Adjoint\n");
 	adjointcore(femmodel);
+
 	if(VerboseControl()) _printf0_("   Compute Initial cost function\n");
-	femmodel->CostFunctionx(&j0,NULL,NULL);
+	femmodel->CostFunctionx(&j0,&jlist,NULL);
+	_printf0_("Initial J(x+dk)      |  List of contributions\n");
+	_printf0_("_____________________________________________\n");
+	_printf0_("J(x) = "<<setw(12)<<setprecision(7)<<j0<<"  |  ");
+	for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<jlist[i]);
+	_printf0_("\n");
+	xDelete<IssmDouble>(jlist);
+
 	if(VerboseControl()) _printf0_("   Compute Gradient\n");
 	Gradjx(&G,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 	for(int i=0;i<n;i++) G[i] = -G[i];
@@ -94,6 +105,8 @@
 	#else
 	femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,output,num,2,1,0));
 	#endif
+	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
+	femmodel->OutputControlsx(&femmodel->results);
 
 	/*Clean up and return*/
 	xDelete<IssmDouble>(G);
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18402)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18403)
@@ -591,7 +591,12 @@
 				case RheologyBbarAbsGradientEnum:   RheologyBbarAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters);    break;
 				case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break;
 				case BalancethicknessMisfitEnum:    BalancethicknessMisfitx(&double_result);                                                        break;
-				case Balancethickness2MisfitEnum:  Balancethickness2Misfitx(&double_result); break;
+				case Balancethickness2MisfitEnum:   Balancethickness2Misfitx(&double_result); break;
+				case IrrotationalH2MisfitEnum:      IrrotationalH2Misfitx(&double_result); break;
+				case IrrotationalDirectionMisfitEnum:  IrrotationalDirectionMisfitx(&double_result); break;
+				case IrrotationalVelMisfitEnum:        IrrotationalVelMisfitx(&double_result); break;
+				case IrrotationalAlongGradientNuEnum:  IrrotationalAlongGradientNux(&double_result); break;
+				case IrrotationalAcrossGradientNuEnum: IrrotationalAcrossGradientNux(&double_result); break;
 
 			   /*Vector */
 				default:
@@ -1304,7 +1309,7 @@
 	IssmDouble J_sum;
 
 	IssmDouble  weight,thicknessobs,thickness,potential,dpotential[2];
-	IssmDouble  vx,vy,vxbar,vybar,vxobs,vyobs,vxbarobs,vybarobs,nux,nuy;
+	IssmDouble  vx,vy,vxbar,vybar,vxobs,vyobs,vxbarobs,vybarobs,nu;
 	IssmDouble  Jdet;
 	IssmDouble* xyz_list = NULL;
 
@@ -1325,8 +1330,7 @@
 		Input* vyobs_input       =element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
 		Input* vx_input          =element->GetInput(VxEnum); _assert_(vx_input);
 		Input* vy_input          =element->GetInput(VyEnum); _assert_(vy_input);
-		Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
-		Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
+		Input* nu_input       = element->GetInput(BalancethicknessNuEnum);   _assert_(nu_input);
 
 		/* Start  looping on the number of gaussian points: */
 		Gauss* gauss=element->NewGauss(2);
@@ -1347,11 +1351,10 @@
 			vyobs_input->GetInputValue(&vyobs,gauss);
 			vx_input->GetInputValue(&vxbar,gauss);
 			vy_input->GetInputValue(&vybar,gauss);
-			nux_input->GetInputValue(&nux,gauss);
-			nuy_input->GetInputValue(&nuy,gauss);
+			nu_input->GetInputValue(&nu,gauss);
 
-			vxbarobs = nux*vxobs;
-			vybarobs = nuy*vyobs;
+			vxbarobs = nu*vxobs;
+			vybarobs = nu*vyobs;
 
 			/*J = (H^2 - Hobs^2)^2*/
 			J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
@@ -1379,6 +1382,126 @@
 	*presponse=J;
 
 }/*}}}*/
+void FemModel::IrrotationalH2Misfitx(IssmDouble* presponse){/*{{{*/
+
+	/*output: */
+	IssmDouble J=0.;
+	IssmDouble J_sum;
+
+	IssmDouble  weight,thicknessobs,thickness;
+	IssmDouble  Jdet;
+	IssmDouble* xyz_list = NULL;
+
+	/*Compute Misfit: */
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+
+		/*If on water, return 0: */
+		if(!element->IsIceInElement()) continue;
+
+		/* Get node coordinates*/
+		element->GetVerticesCoordinates(&xyz_list);
+		Input* weights_input     =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+		Input* thickness_input   =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
+		Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
+
+		/* Start  looping on the number of gaussian points: */
+		Gauss* gauss=element->NewGauss(2);
+		for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+			gauss->GaussPoint(ig);
+			element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+
+			/*Get all parameters at gaussian point*/
+			weights_input->GetInputValue(&weight,gauss,IrrotationalH2MisfitEnum);
+			thickness_input->GetInputValue(&thickness,gauss);
+			thicknessobs_input->GetInputValue(&thicknessobs,gauss);
+
+			/*J = (H^2 - Hobs^2)^2*/
+			J +=0.5*(thickness*thickness - thicknessobs*thicknessobs)*(thickness*thickness - thicknessobs*thicknessobs)*weight*Jdet*gauss->weight;
+		}
+
+		/*clean up and Return: */
+		xDelete<IssmDouble>(xyz_list);
+		delete gauss;
+	}
+
+	/*Sum all J from all cpus of the cluster:*/
+	ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	J=J_sum;
+
+	/*Assign output pointers: */
+	*presponse=J;
+
+}/*}}}*/
+void FemModel::IrrotationalDirectionMisfitx(IssmDouble* presponse){/*{{{*/
+
+	/*output: */
+	IssmDouble J=0.;
+	IssmDouble J_sum;
+
+	IssmDouble  weight,thicknessobs,thickness,dpotential[2];
+	IssmDouble  vx,vy,vxobs,vyobs,nu;
+	IssmDouble  Jdet;
+	IssmDouble* xyz_list = NULL;
+
+	/*Compute Misfit: */
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+
+		/*If on water, return 0: */
+		if(!element->IsIceInElement()) continue;
+
+		/* Get node coordinates*/
+		element->GetVerticesCoordinates(&xyz_list);
+		Input* weights_input     =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+		Input* potential_input   =element->GetInput(PotentialEnum);                          _assert_(potential_input);
+		Input* vxobs_input       =element->GetInput(BalancethicknessVxObsEnum);              _assert_(vxobs_input);
+		Input* vyobs_input       =element->GetInput(BalancethicknessVyObsEnum);              _assert_(vyobs_input);
+		Input* nu_input          = element->GetInput(BalancethicknessNuEnum);                _assert_(nu_input);
+
+		/* Start  looping on the number of gaussian points: */
+		Gauss* gauss=element->NewGauss(2);
+		for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+			gauss->GaussPoint(ig);
+			element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+
+			/*Get all parameters at gaussian point*/
+			weights_input->GetInputValue(&weight,gauss,IrrotationalDirectionMisfitEnum);
+			potential_input->GetInputDerivativeValue(&dpotential[0],xyz_list,gauss);
+			vxobs_input->GetInputValue(&vxobs,gauss);
+			vyobs_input->GetInputValue(&vyobs,gauss);
+			nu_input->GetInputValue(&nu,gauss);
+
+			/*J = 1/2 (vbar ^ gard(phi))^2*/
+			J += 0.5*(nu*vyobs*dpotential[0] - nu*vxobs*dpotential[1])*(nu*vyobs*dpotential[0] - nu*vxobs*dpotential[1])*weight*Jdet*gauss->weight;
+		}
+
+		/*clean up and Return: */
+		xDelete<IssmDouble>(xyz_list);
+		delete gauss;
+	}
+
+	/*Sum all J from all cpus of the cluster:*/
+	ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	J=J_sum;
+
+	/*Assign output pointers: */
+	*presponse=J;
+
+}/*}}}*/
+void FemModel::IrrotationalVelMisfitx(IssmDouble* presponse){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
+void FemModel::IrrotationalAlongGradientNux(IssmDouble* presponse){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
+void FemModel::IrrotationalAcrossGradientNux(IssmDouble* presponse){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
 void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/
 
 	/*output: */
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18402)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18403)
@@ -1031,8 +1031,7 @@
 				/*No yts conversion*/
 				case ThicknessEnum:
 				case FrictionCoefficientEnum:
-				case BalancethicknessNuxEnum:
-				case BalancethicknessNuyEnum:
+				case BalancethicknessNuEnum:
 					if(iomodel->Data(control)){
 						for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(control)[tria_vertex_ids[j]-1];
 						for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18402)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18403)
@@ -1053,8 +1053,7 @@
 				name==BedEnum ||
 				name==BalancethicknessThickeningRateEnum ||
 				name==BalancethicknessApparentMassbalanceEnum ||
-				name==BalancethicknessNuxEnum ||
-				name==BalancethicknessNuyEnum ||
+				name==BalancethicknessNuEnum ||
 				name==SigmaNNEnum ||
 				name==SurfaceSlopeXEnum ||
 				name==SurfaceSlopeYEnum ||
Index: ../trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.h	(revision 18402)
+++ ../trunk-jpl/src/c/classes/FemModel.h	(revision 18403)
@@ -80,6 +80,11 @@
 		void ElementResponsex(IssmDouble* presponse,int response_enum);
 		void BalancethicknessMisfitx(IssmDouble* pV);
 		void Balancethickness2Misfitx(IssmDouble* pV);
+		void IrrotationalH2Misfitx(IssmDouble* pV);
+		void IrrotationalDirectionMisfitx(IssmDouble* pV);
+		void IrrotationalVelMisfitx(IssmDouble* pV);
+		void IrrotationalAlongGradientNux(IssmDouble* pV);
+		void IrrotationalAcrossGradientNux(IssmDouble* pV);
 		#ifdef  _HAVE_DAKOTA_
 		void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);
 		#endif
Index: ../trunk-jpl/src/m/classes/inversion.m
===================================================================
--- ../trunk-jpl/src/m/classes/inversion.m	(revision 18402)
+++ ../trunk-jpl/src/m/classes/inversion.m	(revision 18403)
@@ -131,7 +131,8 @@
 			md = checkfield(md,'fieldname','inversion.iscontrol','values',[0 1]);
 			md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0 1]);
 			md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',...
-				{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar' 'Vx' 'Vy' 'Thickness' 'BalancethicknessNux' 'BalancethicknessNuy' 'BalancethicknessApparentMassbalance'});
+				{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar' 'Vx' 'Vy' 'Thickness',...
+				'BalancethicknessNux' 'BalancethicknessNuy' 'BalancethicknessApparentMassbalance'});
 			md = checkfield(md,'fieldname','inversion.nsteps','numel',1,'>=',0);
 			md = checkfield(md,'fieldname','inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0);
 			md = checkfield(md,'fieldname','inversion.step_threshold','size',[md.inversion.nsteps 1]);
Index: ../trunk-jpl/src/m/classes/inversionvalidation.m
===================================================================
--- ../trunk-jpl/src/m/classes/inversionvalidation.m	(revision 18402)
+++ ../trunk-jpl/src/m/classes/inversionvalidation.m	(revision 18403)
@@ -52,8 +52,9 @@
 			md = checkfield(md,'fieldname','inversion.iscontrol','values',[0 1]);
 			md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0 1]);
 			md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',...
-				{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness' 'BalancethicknessApparentMassbalance'});
-			md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:507]);
+				{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar' 'Vx' 'Vy' 'Thickness',...
+				'BalancethicknessNu' 'BalancethicknessApparentMassbalance'});
+			md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506 601:604]);
 			md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
 			md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
 			md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);
@@ -138,7 +139,11 @@
 			pos=find(obj.cost_functions==504); data(pos)=ThicknessAlongGradientEnum();
 			pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum();
 			pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum();
-			pos=find(obj.cost_functions==507); data(pos)=Balancethickness2MisfitEnum();
+			pos=find(obj.cost_functions==601); data(pos)=IrrotationalH2MisfitEnum();
+			pos=find(obj.cost_functions==602); data(pos)=IrrotationalDirectionMisfitEnum();
+			pos=find(obj.cost_functions==603); data(pos)=IrrotationalVelMisfitEnum();
+			pos=find(obj.cost_functions==604); data(pos)=IrrotationalAlongGradientNuEnum();
+			pos=find(obj.cost_functions==605); data(pos)=IrrotationalAcrossGradientNuEnum();
 			WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
 			WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');
 		end % }}}
Index: ../trunk-jpl/src/m/classes/balancethickness2.m
===================================================================
--- ../trunk-jpl/src/m/classes/balancethickness2.m	(revision 18402)
+++ ../trunk-jpl/src/m/classes/balancethickness2.m	(revision 18403)
@@ -7,8 +7,7 @@
 	properties (SetAccess=public) 
 		spcpotential         = NaN;
 		apparent_massbalance = NaN;
-		nux                  = NaN;
-		nuy                  = NaN;
+		nu                  = NaN;
 		vx_obs               = NaN;
 		vy_obs               = NaN;
 		thickness_obs        = NaN;
@@ -31,8 +30,7 @@
 
 			md = checkfield(md,'fieldname','balancethickness.spcpotential','size',[md.mesh.numberofvertices 1]);
 			md = checkfield(md,'fieldname','balancethickness.apparent_massbalance','size',[md.mesh.numberofvertices 1],'NaN',1);
-			md = checkfield(md,'fieldname','balancethickness.nux','size',[md.mesh.numberofvertices 1],'NaN',1,'>',0,'<=',1);
-			md = checkfield(md,'fieldname','balancethickness.nuy','size',[md.mesh.numberofvertices 1],'NaN',1,'>',0,'<=',1);
+			md = checkfield(md,'fieldname','balancethickness.nu','size',[md.mesh.numberofvertices 1],'NaN',1,'>',0,'<=',1);
 			md = checkfield(md,'fieldname','balancethickness.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
 			md = checkfield(md,'fieldname','balancethickness.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
 			md = checkfield(md,'fieldname','balancethickness.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'>=',0);
@@ -42,8 +40,7 @@
 
 			fielddisplay(obj,'spcpotential','potential constraints (NaN means no constraint)');
 			fielddisplay(obj,'apparent_massbalance','Apparent mass balance [m/yr]');
-			fielddisplay(obj,'nux','vx_bar = nux vx_s (in ]0 1])');
-			fielddisplay(obj,'nuy','vy_bar = nuy vy_s (in ]0 1])');
+			fielddisplay(obj,'nu','v_bar = nu v_s (in ]0 1])');
 			fielddisplay(obj,'vx_obs','observed vx');
 			fielddisplay(obj,'vy_obs','observed vy');
 			fielddisplay(obj,'thickness_obs','observed H');
@@ -52,10 +49,14 @@
 		function marshall(obj,md,fid) % {{{
 			WriteData(fid,'object',obj,'class','Balancethickness','fieldname','spcpotential','format','DoubleMat','mattype',1);
 			WriteData(fid,'object',obj,'class','Balancethickness','fieldname','apparent_massbalance','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'class','Balancethickness','fieldname','nux','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'class','Balancethickness','fieldname','nuy','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'class','Balancethickness','fieldname','vx_obs','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'class','Balancethickness','fieldname','vy_obs','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'class','Balancethickness','fieldname','nu','format','DoubleMat','mattype',1);
+			if(numel(obj.vx_obs)==md.mesh.numberofelements),
+				mattype=2;
+			else
+				mattype=1;
+			end
+			WriteData(fid,'object',obj,'class','Balancethickness','fieldname','vx_obs','format','DoubleMat','mattype',mattype);
+			WriteData(fid,'object',obj,'class','Balancethickness','fieldname','vy_obs','format','DoubleMat','mattype',mattype);
 			WriteData(fid,'object',obj,'class','Balancethickness','fieldname','thickness_obs','format','DoubleMat','mattype',1);
 		end % }}}
 	end
Index: ../trunk-jpl/src/m/enum/BalancethicknessNuxEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BalancethicknessNuxEnum.m	(revision 18402)
+++ ../trunk-jpl/src/m/enum/BalancethicknessNuxEnum.m	(revision 18403)
@@ -1,11 +0,0 @@
-function macro=BalancethicknessNuxEnum()
-%BALANCETHICKNESSNUXENUM - Enum of BalancethicknessNux
-%
-%   WARNING: DO NOT MODIFY THIS FILE
-%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
-%            Please read src/c/shared/Enum/README for more information
-%
-%   Usage:
-%      macro=BalancethicknessNuxEnum()
-
-macro=StringToEnum('BalancethicknessNux');
Index: ../trunk-jpl/src/m/enum/BalancethicknessNuyEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BalancethicknessNuyEnum.m	(revision 18402)
+++ ../trunk-jpl/src/m/enum/BalancethicknessNuyEnum.m	(revision 18403)
@@ -1,11 +0,0 @@
-function macro=BalancethicknessNuyEnum()
-%BALANCETHICKNESSNUYENUM - Enum of BalancethicknessNuy
-%
-%   WARNING: DO NOT MODIFY THIS FILE
-%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
-%            Please read src/c/shared/Enum/README for more information
-%
-%   Usage:
-%      macro=BalancethicknessNuyEnum()
-
-macro=StringToEnum('BalancethicknessNuy');
Index: ../trunk-jpl/src/m/enum/IrrotationalDirectionMisfitEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/IrrotationalDirectionMisfitEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/IrrotationalDirectionMisfitEnum.m	(revision 18403)
@@ -0,0 +1,11 @@
+function macro=IrrotationalDirectionMisfitEnum()
+%IRROTATIONALDIRECTIONMISFITENUM - Enum of IrrotationalDirectionMisfit
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=IrrotationalDirectionMisfitEnum()
+
+macro=StringToEnum('IrrotationalDirectionMisfit');
Index: ../trunk-jpl/src/m/enum/IrrotationalAcrossGradientNuEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/IrrotationalAcrossGradientNuEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/IrrotationalAcrossGradientNuEnum.m	(revision 18403)
@@ -0,0 +1,11 @@
+function macro=IrrotationalAcrossGradientNuEnum()
+%IRROTATIONALACROSSGRADIENTNUENUM - Enum of IrrotationalAcrossGradientNu
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=IrrotationalAcrossGradientNuEnum()
+
+macro=StringToEnum('IrrotationalAcrossGradientNu');
Index: ../trunk-jpl/src/m/enum/IrrotationalAlongGradientNuEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/IrrotationalAlongGradientNuEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/IrrotationalAlongGradientNuEnum.m	(revision 18403)
@@ -0,0 +1,11 @@
+function macro=IrrotationalAlongGradientNuEnum()
+%IRROTATIONALALONGGRADIENTNUENUM - Enum of IrrotationalAlongGradientNu
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=IrrotationalAlongGradientNuEnum()
+
+macro=StringToEnum('IrrotationalAlongGradientNu');
Index: ../trunk-jpl/src/m/enum/IrrotationalH2MisfitEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/IrrotationalH2MisfitEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/IrrotationalH2MisfitEnum.m	(revision 18403)
@@ -0,0 +1,11 @@
+function macro=IrrotationalH2MisfitEnum()
+%IRROTATIONALH2MISFITENUM - Enum of IrrotationalH2Misfit
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=IrrotationalH2MisfitEnum()
+
+macro=StringToEnum('IrrotationalH2Misfit');
Index: ../trunk-jpl/src/m/enum/IrrotationalVelMisfitEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/IrrotationalVelMisfitEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/IrrotationalVelMisfitEnum.m	(revision 18403)
@@ -0,0 +1,11 @@
+function macro=IrrotationalVelMisfitEnum()
+%IRROTATIONALVELMISFITENUM - Enum of IrrotationalVelMisfit
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=IrrotationalVelMisfitEnum()
+
+macro=StringToEnum('IrrotationalVelMisfit');
Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18402)
+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18403)
@@ -300,8 +300,12 @@
 def BalancethicknessSpcpotentialEnum(): return StringToEnum("BalancethicknessSpcpotential")[0]
 def BalancethicknessApparentMassbalanceEnum(): return StringToEnum("BalancethicknessApparentMassbalance")[0]
 def Balancethickness2MisfitEnum(): return StringToEnum("Balancethickness2Misfit")[0]
-def BalancethicknessNuxEnum(): return StringToEnum("BalancethicknessNux")[0]
-def BalancethicknessNuyEnum(): return StringToEnum("BalancethicknessNuy")[0]
+def IrrotationalH2MisfitEnum(): return StringToEnum("IrrotationalH2Misfit")[0]
+def IrrotationalDirectionMisfitEnum(): return StringToEnum("IrrotationalDirectionMisfit")[0]
+def IrrotationalVelMisfitEnum(): return StringToEnum("IrrotationalVelMisfit")[0]
+def IrrotationalAlongGradientNuEnum(): return StringToEnum("IrrotationalAlongGradientNu")[0]
+def IrrotationalAcrossGradientNuEnum(): return StringToEnum("IrrotationalAcrossGradientNu")[0]
+def BalancethicknessNuEnum(): return StringToEnum("BalancethicknessNu")[0]
 def BalancethicknessVxObsEnum(): return StringToEnum("BalancethicknessVxObs")[0]
 def BalancethicknessVyObsEnum(): return StringToEnum("BalancethicknessVyObs")[0]
 def BalancethicknessThicknessObsEnum(): return StringToEnum("BalancethicknessThicknessObs")[0]
Index: ../trunk-jpl/src/m/enum/BalancethicknessNuEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BalancethicknessNuEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/BalancethicknessNuEnum.m	(revision 18403)
@@ -0,0 +1,11 @@
+function macro=BalancethicknessNuEnum()
+%BALANCETHICKNESSNUENUM - Enum of BalancethicknessNu
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=BalancethicknessNuEnum()
+
+macro=StringToEnum('BalancethicknessNu');
