Index: ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18475)
+++ ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18476)
@@ -44,108 +44,8 @@
 	return Ke;
 }/*}}}*/
 ElementVector* AdjointBalancethickness2Analysis::CreatePVector(Element* element){/*{{{*/
+	_error_("Not implemented yet");
 
-	/*Intermediaries */
-	int         num_responses,i;
-	IssmDouble  hobs,hu2,weight,Jdet;
-	IssmDouble  NUMxH2,NUMyH2,DENH2;
-	IssmDouble  NUMxUbar,NUMyUbar,DENUbar;
-	IssmDouble  vxobs,vyobs,vxobsbar,vyobsbar,vbarobs2,vbarobs;
-	IssmDouble  nu,phi,dphi[2];
-	int        *responses = NULL;
-	IssmDouble *xyz_list  = NULL;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Initialize Element vector and vectors*/
-	ElementVector* pe     = element->NewElementVector(SSAApproximationEnum);
-	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
-	element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
-	element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
-	Input* thicknessobs_input = element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
-	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);
-		element->NodalFunctions(basis,gauss);
-		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-		vxobs_input->GetInputValue(&vxobs,gauss);
-		vyobs_input->GetInputValue(&vyobs,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 = nu*vxobs;
-		vyobsbar = nu*vyobs;
-
-		vbarobs2 = (nu*nu*vxobs*vxobs + nu*nu*vyobs*vyobs);
-		vbarobs  = sqrt(vbarobs2);
-		hu2 = hobs*hobs*vbarobs2;
-
-		/*Loop over all requested responses*/
-		for(int resp=0;resp<num_responses;resp++){
-			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 = 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 - 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;
-					//}
-					_error_("Not implemented yet");
-					break;
-				default:
-					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
-			}
-		}
-	}
-
-	/*Clean up and return*/
-	xDelete<int>(responses);
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(basis);
-	xDelete<IssmDouble>(dbasis);
-	delete gauss;
-	return pe;
 }/*}}}*/
 void AdjointBalancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	_error_("not implemented yet");
@@ -173,34 +73,11 @@
 		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 BalancethicknessNuEnum: /*Nothing: state equations do not depend on k*/ break;
 		default: _error_("control type not supported yet: " << EnumToStringx(control_type));
 	}
 
@@ -250,120 +127,6 @@
 	xDelete<int>(vertexpidlist);
 	delete gauss;
 }/*}}}*/
-void AdjointBalancethickness2Analysis::GradientJ1Nu(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble vxobs,vyobs,thicknessobs,weight;
-	IssmDouble nu,Jdet,dphi[2],dphinorm2,nuvobs2;
-	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* 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,IrrotationalH2MisfitEnum);
-		vxobs_input->GetInputValue(&vxobs,gauss);
-		vyobs_input->GetInputValue(&vyobs,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   = 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*basis[i]*(
-						-2./pow(nu,3) * dphinorm2/(vxobs*vxobs + vyobs*vyobs) * (dphinorm2/nuvobs2 - thicknessobs*thicknessobs) 
-						);
-			_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::GradientJ2Nu(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble vxobs,vyobs,thicknessobs,weight;
-	IssmDouble nu,Jdet,dphi[2],dphinorm2,nuvobs2;
-	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* 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,IrrotationalDirectionMisfitEnum);
-		vxobs_input->GetInputValue(&vxobs,gauss);
-		vyobs_input->GetInputValue(&vyobs,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   = 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*nu*(vyobs*dphi[0] -vxobs*dphi[1])*(vyobs*dphi[0] -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: ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.h	(revision 18475)
+++ ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.h	(revision 18476)
@@ -25,11 +25,12 @@
 		ElementMatrix* CreateJacobianMatrix(Element* element);
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		ElementVector* CreatePVectorVolume(Element* element);
-		ElementVector* CreatePVectorBoundary(Element* element);
 		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
 		void UpdateConstraints(FemModel* femmodel);
+
+		/*Specifics*/
+		void CreateDiffusionCoefficient(Element* element);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h	(revision 18475)
+++ ../trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h	(revision 18476)
@@ -28,8 +28,6 @@
 		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 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 18475)
+++ ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 18476)
@@ -25,13 +25,13 @@
 		}
 	}
 
-	iomodel->FetchDataToInput(elements,BalancethicknessApparentMassbalanceEnum);
-	iomodel->FetchDataToInput(elements,BalancethicknessNuEnum);
-	iomodel->FetchDataToInput(elements,BalancethicknessVxObsEnum);
-	iomodel->FetchDataToInput(elements,BalancethicknessVyObsEnum);
-	iomodel->FetchDataToInput(elements,BalancethicknessThicknessObsEnum);
-	iomodel->FetchDataToInput(elements,MeshVertexonboundaryEnum);
+	iomodel->FetchDataToInput(elements,ThicknessEnum);
+	iomodel->FetchDataToInput(elements,SurfaceEnum);
+	iomodel->FetchDataToInput(elements,BaseEnum);
 	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
+	iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
+	iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
+	iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum);
 }/*}}}*/
 void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
@@ -42,7 +42,7 @@
 
 
 	int finiteelement = P1Enum;
-	IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcpotentialEnum,Balancethickness2AnalysisEnum,finiteelement);
+	IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,Balancethickness2AnalysisEnum,finiteelement);
 
 }/*}}}*/
 void Balancethickness2Analysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
@@ -63,7 +63,7 @@
 ElementMatrix* Balancethickness2Analysis::CreateKMatrix(Element* element){/*{{{*/
 
 	/*Intermediaries */
-	IssmDouble Jdet,D_scalar;
+	IssmDouble  Jdet,D;
 	IssmDouble* xyz_list = NULL;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -73,8 +73,12 @@
 	ElementMatrix* Ke     = element->NewElementMatrix();
 	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
 
+	/*Create input D*/
+	this->CreateDiffusionCoefficient(element);
+
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
+	Input* D_input = element->GetInput(BalancethicknessDiffusionCoefficientEnum); _assert_(D_input);
 
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(2);
@@ -82,10 +86,11 @@
 		gauss->GaussPoint(ig);
 		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		D_input->GetInputValue(&D,gauss);
 
 		for(int i=0;i<numnodes;i++){
 			for(int j=0;j<numnodes;j++){
-				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
+				Ke->values[i*numnodes+j] += D*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
 			}
 		}
 	}
@@ -98,20 +103,8 @@
 }/*}}}*/
 ElementVector* Balancethickness2Analysis::CreatePVector(Element* element){/*{{{*/
 
-	/*compute all load vectors for this element*/
-	ElementVector* pe1=CreatePVectorVolume(element);
-	ElementVector* pe2=CreatePVectorBoundary(element);
-	ElementVector* pe =new ElementVector(pe1,pe2);
-
-	/*clean-up and return*/
-	delete pe1;
-	delete pe2;
-	return pe;
-}/*}}}*/
-ElementVector* Balancethickness2Analysis::CreatePVectorVolume(Element* element){/*{{{*/
-
 	/*Intermediaries */
-	IssmDouble  adot,Jdet;
+	IssmDouble  dhdt,mb,ms,dD[2],db[2],Jdet;
 	IssmDouble* xyz_list = NULL;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -123,7 +116,11 @@
 
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* adot_input   = element->GetInput(BalancethicknessApparentMassbalanceEnum); _assert_(adot_input);
+	Input* ms_input   = element->GetInput(SurfaceforcingsMassBalanceEnum);                _assert_(ms_input);
+	Input* mb_input   = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);       _assert_(mb_input);
+	Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum);            _assert_(dhdt_input);
+	Input* D_input   = element->GetInput(BalancethicknessDiffusionCoefficientEnum);       _assert_(D_input);
+	Input* bed_input  = element->GetInput(BaseEnum);                                      _assert_(bed_input);
 
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(2);
@@ -132,61 +129,21 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctions(basis,gauss);
-		adot_input->GetInputValue(&adot,gauss);
 
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*adot*basis[i];
-	}
+		ms_input->GetInputValue(&ms,gauss);
+		mb_input->GetInputValue(&mb,gauss);
+		dhdt_input->GetInputValue(&dhdt,gauss);
+		bed_input->GetInputDerivativeValue(&db[0],xyz_list,gauss);
+		D_input->GetInputDerivativeValue(&dD[0],xyz_list,gauss);
+		db[0]=0.;
+		db[1]=0.;
 
-	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(basis);
-	delete gauss;
-	return pe;
-}/*}}}*/
-ElementVector* Balancethickness2Analysis::CreatePVectorBoundary(Element* element){/*{{{*/
-
-	/*If no front, return NULL*/
-	if(!element->IsFaceOnBoundary()) return NULL;
-
-	/*Intermediaries*/
-	IssmDouble  Jdet,thickness,vx,vy;
-	IssmDouble *xyz_list = NULL;
-	IssmDouble *xyz_list_front = NULL;
-	IssmDouble  normal[2];
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Initialize Element vector and other vectors*/
-	ElementVector* pe    = element->NewElementVector();
-	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-	/*Retrieve all inputs and parameters*/
-	Input* thickness_input = element->GetInput(BalancethicknessThicknessObsEnum); _assert_(thickness_input);
-	Input* vx_input        = element->GetInput(BalancethicknessVxObsEnum);        _assert_(vx_input);
-	Input* vy_input        = element->GetInput(BalancethicknessVyObsEnum);        _assert_(vy_input);
-
-	element->GetVerticesCoordinates(&xyz_list);
-	element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
-	element->NormalSection(&normal[0],xyz_list_front);
-
-	/*Start looping on Gaussian points*/
-	Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-		thickness_input->GetInputValue(&thickness,gauss);
-		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
-		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
-		element->NodalFunctions(basis,gauss);
-
-		for(int i=0;i<numnodes;i++) pe->values[i] += Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1])*basis[i];
+		/*Since grad(b) is constant div(D grad(b) ) = grad(D).grad(b)*/
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt +dD[0]*db[0]+dD[1]*db[1])*basis[i];
 	}
 
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(xyz_list_front);
 	xDelete<IssmDouble>(basis);
 	delete gauss;
 	return pe;
@@ -198,79 +155,55 @@
 	_error_("Not implemented yet");
 }/*}}}*/
 void Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+	element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
+}/*}}}*/
+void Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+	/*Default, do nothing*/
+	return;
+}/*}}}*/
 
+/*Specifics*/
+void Balancethickness2Analysis::CreateDiffusionCoefficient(Element* element){/*{{{*/
+
 	/*Intermediaries */
-	int         Hinterpolation;
-	IssmDouble  vx,vy,vbar,nu,normdphi,dphi[2],H;
-	IssmDouble* xyz_list = NULL;
-	int       * doflist  = NULL;
+	IssmDouble       omega,h,mu0,ds[2],Cmu,B;
+	const int        n = 3.;
+	const IssmDouble Hstar = 500.;
+	const IssmDouble Lstar = 500.e+3;
+	IssmDouble *xyz_list  = NULL;
 
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes    = element->GetNumberOfNodes();
-	int numvertices = element->GetNumberOfVertices();
+	/*Fetch number of vertices and allocate output*/
+	int  numnodes = element->GetNumberOfNodes();
+	IssmDouble* D      = xNew<IssmDouble>(numnodes);
+	IssmDouble* Dgradb = xNew<IssmDouble>(numnodes);
 
-	/*Fetch dof list and allocate solution vector*/
-	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	IssmDouble* values         = xNew<IssmDouble>(numnodes);
-	IssmDouble* thickness_list = xNew<IssmDouble>(numvertices);
-	IssmDouble* vx_list        = xNew<IssmDouble>(numvertices);
-	IssmDouble* vy_list        = xNew<IssmDouble>(numvertices);
-
-	/*Use the dof list to index into the solution vector: */
-	for(int i=0;i<numnodes;i++){
-		values[i]=solution[doflist[i]];
-
-		/*Check solution*/
-		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
-	}
-
-	element->AddInput(PotentialEnum,values,element->GetElementType());
-
-	/*Retrieve all inputs and parameters*/
+	/*retrieve what we need: */
 	element->GetVerticesCoordinates(&xyz_list);
-	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* nu_input        = element->GetInput(BalancethicknessNuEnum);    _assert_(nu_input);
-	Input* thickness_input = element->GetInput(BalancethicknessThicknessObsEnum);             _assert_(thickness_input);
+	Input* thickness_input  = element->GetInput(ThicknessEnum);          _assert_(thickness_input);
+	Input* surface_input    = element->GetInput(SurfaceEnum);            _assert_(surface_input);
+	Input* B_input          = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
+	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
 
-	switch(element->GetElementType()){
-		case P1Enum: Hinterpolation = P0Enum; break;
-		default:     _error_("not implemented");
-	}
-
+	/*Calculate damage evolution source term: */
 	Gauss* gauss=element->NewGauss();
-	for (int iv=0;iv<1;iv++){
-		gauss->GaussNode(Hinterpolation,iv);//P0 Only for now
+	for (int i=0;i<numnodes;i++){
+		gauss->GaussNode(element->GetElementType(),i);
+		
+		B_input->GetInputValue(&B,gauss);
+		thickness_input->GetInputValue(&h,gauss);
+		surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
 
-		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
-		nu_input->GetInputValue(&nu,gauss);
-		thickness_input->GetInputValue(&H,gauss);
-		potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+		mu0   = pow(2.,(1.-3*n)/(2.*n)) * B;
+		omega = pow(rhog,3) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
+		Cmu   = 0.;
 
-		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]);
-
-		thickness_list[iv] = normdphi/vbar;
-		vx_list[iv]        = -1./H * dphi[0];
-		vy_list[iv]        = -1./H * dphi[1];
+		D[i] = omega*(ds[0]*ds[0]+ds[1]*ds[1])*pow(h,4)*(1./5.*h - Cmu);
 	}
-	element->AddInput(ThicknessEnum,thickness_list,Hinterpolation);
-	element->AddInput(VxEnum,vx_list,Hinterpolation);
-	element->AddInput(VyEnum,vy_list,Hinterpolation);
 
+	/*Add input*/
+	element->AddInput(BalancethicknessDiffusionCoefficientEnum,D,element->GetElementType());
+	
 	/*Clean up and return*/
+	xDelete<IssmDouble>(D);
 	delete gauss;
-	xDelete<int>(doflist);
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(values);
-	xDelete<IssmDouble>(thickness_list);
-	xDelete<IssmDouble>(vx_list);
-	xDelete<IssmDouble>(vy_list);
 }/*}}}*/
-void Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
-	/*Default, do nothing*/
-	return;
-}/*}}}*/
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18475)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18476)
@@ -300,15 +300,7 @@
 	BalancethicknessSpcpotentialEnum,
 	BalancethicknessApparentMassbalanceEnum,
 	Balancethickness2MisfitEnum,
-	IrrotationalH2MisfitEnum,
-	IrrotationalDirectionMisfitEnum,
-	IrrotationalVelMisfitEnum,
-	IrrotationalAlongGradientNuEnum,
-	IrrotationalAcrossGradientNuEnum,
-	BalancethicknessNuEnum,
-	BalancethicknessVxObsEnum,
-	BalancethicknessVyObsEnum,
-	BalancethicknessThicknessObsEnum,
+	BalancethicknessDiffusionCoefficientEnum,
 	/*}}}*/
 	/*Surfaceforcings{{{*/
 	SurfaceforcingsEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18475)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18476)
@@ -308,15 +308,7 @@
 		case BalancethicknessSpcpotentialEnum : return "BalancethicknessSpcpotential";
 		case BalancethicknessApparentMassbalanceEnum : return "BalancethicknessApparentMassbalance";
 		case Balancethickness2MisfitEnum : return "Balancethickness2Misfit";
-		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";
+		case BalancethicknessDiffusionCoefficientEnum : return "BalancethicknessDiffusionCoefficient";
 		case SurfaceforcingsEnum : return "Surfaceforcings";
 		case SMBEnum : return "SMB";
 		case SurfaceforcingsMassBalanceEnum : return "SurfaceforcingsMassBalance";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18475)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18476)
@@ -314,15 +314,7 @@
 	      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,"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;
+	      else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
 	      else if (strcmp(name,"Surfaceforcings")==0) return SurfaceforcingsEnum;
 	      else if (strcmp(name,"SMB")==0) return SMBEnum;
 	      else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
@@ -382,10 +374,7 @@
 	      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 stage=4;
-   }
-   if(stage==4){
-	      if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
+	      else 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;
@@ -393,7 +382,10 @@
 	      else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
 	      else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
 	      else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum;
-	      else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
+         else stage=4;
+   }
+   if(stage==4){
+	      if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
 	      else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
 	      else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
 	      else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
@@ -505,10 +497,7 @@
 	      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 stage=5;
-   }
-   if(stage==5){
-	      if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
+	      else 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;
@@ -516,7 +505,10 @@
 	      else if (strcmp(name,"Misfit")==0) return MisfitEnum;
 	      else if (strcmp(name,"Pressure")==0) return PressureEnum;
 	      else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
-	      else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
+         else stage=5;
+   }
+   if(stage==5){
+	      if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
 	      else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
 	      else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum;
 	      else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
@@ -628,10 +620,7 @@
 	      else if (strcmp(name,"MisfitName")==0) return MisfitNameEnum;
 	      else if (strcmp(name,"MisfitModelEnum")==0) return MisfitModelEnumEnum;
 	      else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
-         else stage=6;
-   }
-   if(stage==6){
-	      if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
+	      else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
 	      else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
 	      else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
 	      else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
@@ -639,7 +628,10 @@
 	      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;
+         else stage=6;
+   }
+   if(stage==6){
+	      if (strcmp(name,"MinVx")==0) return MinVxEnum;
 	      else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
 	      else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
 	      else if (strcmp(name,"MinVy")==0) return MinVyEnum;
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp	(revision 18475)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp	(revision 18476)
@@ -40,7 +40,6 @@
 			case VyEnum:
 			case ThicknessEnum:
 			case FrictionCoefficientEnum:
-			case BalancethicknessNuEnum:
 			case BalancethicknessApparentMassbalanceEnum:
 				iomodel->FetchData(1,control); 
 				break;
Index: ../trunk-jpl/src/c/cores/balancethickness2_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/balancethickness2_core.cpp	(revision 18475)
+++ ../trunk-jpl/src/c/cores/balancethickness2_core.cpp	(revision 18476)
@@ -25,8 +25,8 @@
 
 	if(save_results){
 		if(VerboseSolution()) _printf0_("   saving results\n");
-		int outputs[4] = {ThicknessEnum,PotentialEnum,VxEnum,VyEnum};
-		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4);
+		int outputs[1] = {ThicknessEnum};
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
 	}
 
 }
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18475)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18476)
@@ -592,11 +592,6 @@
 				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 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:
@@ -1319,58 +1314,8 @@
 
 		/*If on water, return 0: */
 		if(!element->IsIceInElement()) continue;
+		 _error_("Not implemented");
 
-		/* 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);
-		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* vx_input          =element->GetInput(VxEnum); _assert_(vx_input);
-		Input* vy_input          =element->GetInput(VyEnum); _assert_(vy_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);
-
-			/* Get Jacobian determinant: */
-			element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-
-			/*Get all parameters at gaussian point*/
-			weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
-			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);
-			vy_input->GetInputValue(&vybar,gauss);
-			nu_input->GetInputValue(&nu,gauss);
-
-			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;
-			/*J = phi^2*/
-			//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 = 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;
-		}
-
-		/*clean up and Return: */
-		xDelete<IssmDouble>(xyz_list);
-		delete gauss;
 	}
 
 	/*Sum all J from all cpus of the cluster:*/
@@ -1382,126 +1327,6 @@
 	*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 18475)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18476)
@@ -1031,7 +1031,6 @@
 				/*No yts conversion*/
 				case ThicknessEnum:
 				case FrictionCoefficientEnum:
-				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 18475)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18476)
@@ -1052,8 +1052,6 @@
 				name==BaseEnum ||
 				name==BedEnum ||
 				name==BalancethicknessThickeningRateEnum ||
-				name==BalancethicknessApparentMassbalanceEnum ||
-				name==BalancethicknessNuEnum ||
 				name==SigmaNNEnum ||
 				name==SurfaceSlopeXEnum ||
 				name==SurfaceSlopeYEnum ||
Index: ../trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.h	(revision 18475)
+++ ../trunk-jpl/src/c/classes/FemModel.h	(revision 18476)
@@ -80,11 +80,6 @@
 		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/dox/issm.dox
===================================================================
--- ../trunk-jpl/src/dox/issm.dox	(revision 18475)
+++ ../trunk-jpl/src/dox/issm.dox	(revision 18476)
@@ -46,44 +46,29 @@
 <th  bgcolor=#7AA9DD style="text-align:left;">Language</th><th  bgcolor=#7AA9DD style="text-align:right;">files</th><th  bgcolor=#7AA9DD style="text-align:right;">blank</th><th  bgcolor=#7AA9DD style="text-align:right;">comment</th><th  bgcolor=#7AA9DD style="text-align:right;">code</th><th  bgcolor=#7AA9DD style="text-align:right;">Total</th>
 </th>
 <tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> C++ </th><td  bgcolor=#FFFFFF style="text-align:right;">453</td><td  bgcolor=#FFFFFF style="text-align:right;">15861</td><td  bgcolor=#FFFFFF style="text-align:right;">18188</td><td  bgcolor=#FFFFFF style="text-align:right;">68849</td><td  bgcolor=#FFFFFF style="text-align:right;">102898</td>
+<th  bgcolor=#FFFFFF style="text-align:left;"> C++ </th><td  bgcolor=#FFFFFF style="text-align:right;">445</td><td  bgcolor=#FFFFFF style="text-align:right;">15967</td><td  bgcolor=#FFFFFF style="text-align:right;">16145</td><td  bgcolor=#FFFFFF style="text-align:right;">70397</td><td  bgcolor=#FFFFFF style="text-align:right;">102509</td>
 </tr>
 <tr>
-<th  bgcolor=#C6E2FF style="text-align:left;"> MATLAB </th><td  bgcolor=#C6E2FF style="text-align:right;">1378</td><td  bgcolor=#C6E2FF style="text-align:right;">8341</td><td  bgcolor=#C6E2FF style="text-align:right;">16259</td><td  bgcolor=#C6E2FF style="text-align:right;">38645</td><td  bgcolor=#C6E2FF style="text-align:right;">63245</td>
+<th  bgcolor=#C6E2FF style="text-align:left;"> MATLAB </th><td  bgcolor=#C6E2FF style="text-align:right;">1422</td><td  bgcolor=#C6E2FF style="text-align:right;">8497</td><td  bgcolor=#C6E2FF style="text-align:right;">16638</td><td  bgcolor=#C6E2FF style="text-align:right;">39538</td><td  bgcolor=#C6E2FF style="text-align:right;">64673</td>
 </tr>
 <tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> C/C++  Header </th><td  bgcolor=#FFFFFF style="text-align:right;">421</td><td  bgcolor=#FFFFFF style="text-align:right;">3504</td><td  bgcolor=#FFFFFF style="text-align:right;">3753</td><td  bgcolor=#FFFFFF style="text-align:right;">15489</td><td  bgcolor=#FFFFFF style="text-align:right;">22746</td>
+<th  bgcolor=#FFFFFF style="text-align:left;"> C/C++  Header </th><td  bgcolor=#FFFFFF style="text-align:right;">411</td><td  bgcolor=#FFFFFF style="text-align:right;">3443</td><td  bgcolor=#FFFFFF style="text-align:right;">3528</td><td  bgcolor=#FFFFFF style="text-align:right;">15243</td><td  bgcolor=#FFFFFF style="text-align:right;">22214</td>
 </tr>
 <tr>
-<th  bgcolor=#C6E2FF style="text-align:left;"> m4 </th><td  bgcolor=#C6E2FF style="text-align:right;">9</td><td  bgcolor=#C6E2FF style="text-align:right;">1588</td><td  bgcolor=#C6E2FF style="text-align:right;">151</td><td  bgcolor=#C6E2FF style="text-align:right;">11565</td><td  bgcolor=#C6E2FF style="text-align:right;">13304</td>
+<th  bgcolor=#C6E2FF style="text-align:left;"> m4 </th><td  bgcolor=#C6E2FF style="text-align:right;">8</td><td  bgcolor=#C6E2FF style="text-align:right;">1036</td><td  bgcolor=#C6E2FF style="text-align:right;">149</td><td  bgcolor=#C6E2FF style="text-align:right;">9756</td><td  bgcolor=#C6E2FF style="text-align:right;">10941</td>
 </tr>
 <tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> Python </th><td  bgcolor=#FFFFFF style="text-align:right;">144</td><td  bgcolor=#FFFFFF style="text-align:right;">2330</td><td  bgcolor=#FFFFFF style="text-align:right;">2559</td><td  bgcolor=#FFFFFF style="text-align:right;">9348</td><td  bgcolor=#FFFFFF style="text-align:right;">14237</td>
+<th  bgcolor=#FFFFFF style="text-align:left;"> Python </th><td  bgcolor=#FFFFFF style="text-align:right;">151</td><td  bgcolor=#FFFFFF style="text-align:right;">2421</td><td  bgcolor=#FFFFFF style="text-align:right;">2623</td><td  bgcolor=#FFFFFF style="text-align:right;">9733</td><td  bgcolor=#FFFFFF style="text-align:right;">14777</td>
 </tr>
 <tr>
-<th  bgcolor=#C6E2FF style="text-align:left;"> XML </th><td  bgcolor=#C6E2FF style="text-align:right;">4</td><td  bgcolor=#C6E2FF style="text-align:right;">160</td><td  bgcolor=#C6E2FF style="text-align:right;">92</td><td  bgcolor=#C6E2FF style="text-align:right;">3075</td><td  bgcolor=#C6E2FF style="text-align:right;">3327</td>
-</tr>
-<tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> Java </th><td  bgcolor=#FFFFFF style="text-align:right;">18</td><td  bgcolor=#FFFFFF style="text-align:right;">719</td><td  bgcolor=#FFFFFF style="text-align:right;">891</td><td  bgcolor=#FFFFFF style="text-align:right;">2321</td><td  bgcolor=#FFFFFF style="text-align:right;">3931</td>
-</tr>
-<tr>
 <th  bgcolor=#C6E2FF style="text-align:left;"> Fortran  77 </th><td  bgcolor=#C6E2FF style="text-align:right;">7</td><td  bgcolor=#C6E2FF style="text-align:right;">4</td><td  bgcolor=#C6E2FF style="text-align:right;">302</td><td  bgcolor=#C6E2FF style="text-align:right;">365</td><td  bgcolor=#C6E2FF style="text-align:right;">671</td>
 </tr>
 <tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> Bourne  Shell </th><td  bgcolor=#FFFFFF style="text-align:right;">3</td><td  bgcolor=#FFFFFF style="text-align:right;">61</td><td  bgcolor=#FFFFFF style="text-align:right;">88</td><td  bgcolor=#FFFFFF style="text-align:right;">266</td><td  bgcolor=#FFFFFF style="text-align:right;">415</td>
+<th  bgcolor=#FFFFFF style="text-align:left;"> Bourne  Shell </th><td  bgcolor=#FFFFFF style="text-align:right;">2</td><td  bgcolor=#FFFFFF style="text-align:right;">59</td><td  bgcolor=#FFFFFF style="text-align:right;">84</td><td  bgcolor=#FFFFFF style="text-align:right;">262</td><td  bgcolor=#FFFFFF style="text-align:right;">405</td>
 </tr>
 <tr>
-<th  bgcolor=#C6E2FF style="text-align:left;"> XSD </th><td  bgcolor=#C6E2FF style="text-align:right;">2</td><td  bgcolor=#C6E2FF style="text-align:right;">30</td><td  bgcolor=#C6E2FF style="text-align:right;">36</td><td  bgcolor=#C6E2FF style="text-align:right;">229</td><td  bgcolor=#C6E2FF style="text-align:right;">295</td>
+<th  bgcolor=#C6E2FF style="text-align:left;"> SUM: </th><td  bgcolor=#C6E2FF style="text-align:right;">2446</td><td  bgcolor=#C6E2FF style="text-align:right;">31427</td><td  bgcolor=#C6E2FF style="text-align:right;">39469</td><td  bgcolor=#C6E2FF style="text-align:right;">145294</td><td  bgcolor=#C6E2FF style="text-align:right;">216190</td>
 </tr>
-<tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> Perl </th><td  bgcolor=#FFFFFF style="text-align:right;">1</td><td  bgcolor=#FFFFFF style="text-align:right;">6</td><td  bgcolor=#FFFFFF style="text-align:right;">9</td><td  bgcolor=#FFFFFF style="text-align:right;">196</td><td  bgcolor=#FFFFFF style="text-align:right;">211</td>
-</tr>
-<tr>
-<th  bgcolor=#C6E2FF style="text-align:left;"> Ant </th><td  bgcolor=#C6E2FF style="text-align:right;">1</td><td  bgcolor=#C6E2FF style="text-align:right;">16</td><td  bgcolor=#C6E2FF style="text-align:right;">7</td><td  bgcolor=#C6E2FF style="text-align:right;">103</td><td  bgcolor=#C6E2FF style="text-align:right;">126</td>
-</tr>
-<tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> SUM: </th><td  bgcolor=#FFFFFF style="text-align:right;">2441</td><td  bgcolor=#FFFFFF style="text-align:right;">32620</td><td  bgcolor=#FFFFFF style="text-align:right;">42335</td><td  bgcolor=#FFFFFF style="text-align:right;">150451</td><td  bgcolor=#FFFFFF style="text-align:right;">225406</td>
-</tr>
 </table>
 
 	<I> Copyright (C) 2014 </I>
Index: ../trunk-jpl/src/m/classes/balancethickness2.m
===================================================================
--- ../trunk-jpl/src/m/classes/balancethickness2.m	(revision 18475)
+++ ../trunk-jpl/src/m/classes/balancethickness2.m	(revision 18476)
@@ -1,63 +0,0 @@
-%BALANCETHICKNESS2 class definition
-%
-%   Usage:
-%      balancethickness2=balancethickness2();
-
-classdef balancethickness2
-	properties (SetAccess=public) 
-		spcpotential         = NaN;
-		apparent_massbalance = NaN;
-		nu                  = NaN;
-		vx_obs               = NaN;
-		vy_obs               = NaN;
-		thickness_obs        = NaN;
-	end
-	methods
-		function obj = balancethickness2(varargin) % {{{
-			switch nargin
-				case 0
-					obj=setdefaultparameters(obj);
-				otherwise
-					error('constructor not supported');
-			end
-		end % }}}
-		function obj = setdefaultparameters(obj) % {{{
-
-		end % }}}
-		function md = checkconsistency(obj,md,solution,analyses) % {{{
-			%Early return
-			if solution~=Balancethickness2SolutionEnum(), return; end
-
-			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.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);
-		end % }}}
-		function disp(obj) % {{{
-			disp(sprintf('   balance thickness solution parameters:'));
-
-			fielddisplay(obj,'spcpotential','potential constraints (NaN means no constraint)');
-			fielddisplay(obj,'apparent_massbalance','Apparent mass balance [m/yr]');
-			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');
-
-		end % }}}
-		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','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
-end
Index: ../trunk-jpl/src/m/classes/inversion.m
===================================================================
--- ../trunk-jpl/src/m/classes/inversion.m	(revision 18475)
+++ ../trunk-jpl/src/m/classes/inversion.m	(revision 18476)
@@ -238,11 +238,6 @@
 			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==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/m1qn3inversion.m
===================================================================
--- ../trunk-jpl/src/m/classes/m1qn3inversion.m	(revision 18475)
+++ ../trunk-jpl/src/m/classes/m1qn3inversion.m	(revision 18476)
@@ -169,11 +169,6 @@
 			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==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/inversionvalidation.m
===================================================================
--- ../trunk-jpl/src/m/classes/inversionvalidation.m	(revision 18475)
+++ ../trunk-jpl/src/m/classes/inversionvalidation.m	(revision 18476)
@@ -139,11 +139,6 @@
 			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==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/enum/IrrotationalDirectionMisfitEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/IrrotationalDirectionMisfitEnum.m	(revision 18475)
+++ ../trunk-jpl/src/m/enum/IrrotationalDirectionMisfitEnum.m	(revision 18476)
@@ -1,11 +0,0 @@
-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/BalancethicknessThicknessObsEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BalancethicknessThicknessObsEnum.m	(revision 18475)
+++ ../trunk-jpl/src/m/enum/BalancethicknessThicknessObsEnum.m	(revision 18476)
@@ -1,11 +0,0 @@
-function macro=BalancethicknessThicknessObsEnum()
-%BALANCETHICKNESSTHICKNESSOBSENUM - Enum of BalancethicknessThicknessObs
-%
-%   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=BalancethicknessThicknessObsEnum()
-
-macro=StringToEnum('BalancethicknessThicknessObs');
Index: ../trunk-jpl/src/m/enum/IrrotationalAcrossGradientNuEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/IrrotationalAcrossGradientNuEnum.m	(revision 18475)
+++ ../trunk-jpl/src/m/enum/IrrotationalAcrossGradientNuEnum.m	(revision 18476)
@@ -1,11 +0,0 @@
-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 18475)
+++ ../trunk-jpl/src/m/enum/IrrotationalAlongGradientNuEnum.m	(revision 18476)
@@ -1,11 +0,0 @@
-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/BalancethicknessVxObsEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BalancethicknessVxObsEnum.m	(revision 18475)
+++ ../trunk-jpl/src/m/enum/BalancethicknessVxObsEnum.m	(revision 18476)
@@ -1,11 +0,0 @@
-function macro=BalancethicknessVxObsEnum()
-%BALANCETHICKNESSVXOBSENUM - Enum of BalancethicknessVxObs
-%
-%   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=BalancethicknessVxObsEnum()
-
-macro=StringToEnum('BalancethicknessVxObs');
Index: ../trunk-jpl/src/m/enum/IrrotationalH2MisfitEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/IrrotationalH2MisfitEnum.m	(revision 18475)
+++ ../trunk-jpl/src/m/enum/IrrotationalH2MisfitEnum.m	(revision 18476)
@@ -1,11 +0,0 @@
-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 18475)
+++ ../trunk-jpl/src/m/enum/IrrotationalVelMisfitEnum.m	(revision 18476)
@@ -1,11 +0,0 @@
-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/BalancethicknessVyObsEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BalancethicknessVyObsEnum.m	(revision 18475)
+++ ../trunk-jpl/src/m/enum/BalancethicknessVyObsEnum.m	(revision 18476)
@@ -1,11 +0,0 @@
-function macro=BalancethicknessVyObsEnum()
-%BALANCETHICKNESSVYOBSENUM - Enum of BalancethicknessVyObs
-%
-%   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=BalancethicknessVyObsEnum()
-
-macro=StringToEnum('BalancethicknessVyObs');
Index: ../trunk-jpl/src/m/enum/BalancethicknessNuEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BalancethicknessNuEnum.m	(revision 18475)
+++ ../trunk-jpl/src/m/enum/BalancethicknessNuEnum.m	(revision 18476)
@@ -1,11 +0,0 @@
-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');
Index: ../trunk-jpl/src/m/enum/BalancethicknessDiffusionCoefficientEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BalancethicknessDiffusionCoefficientEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/BalancethicknessDiffusionCoefficientEnum.m	(revision 18476)
@@ -0,0 +1,11 @@
+function macro=BalancethicknessDiffusionCoefficientEnum()
+%BALANCETHICKNESSDIFFUSIONCOEFFICIENTENUM - Enum of BalancethicknessDiffusionCoefficient
+%
+%   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=BalancethicknessDiffusionCoefficientEnum()
+
+macro=StringToEnum('BalancethicknessDiffusionCoefficient');
Index: ../trunk-jpl/src/m/enum/LACrouzeixRaviartEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/LACrouzeixRaviartEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/LACrouzeixRaviartEnum.m	(revision 18476)
@@ -0,0 +1,11 @@
+function macro=LACrouzeixRaviartEnum()
+%LACROUZEIXRAVIARTENUM - Enum of LACrouzeixRaviart
+%
+%   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=LACrouzeixRaviartEnum()
+
+macro=StringToEnum('LACrouzeixRaviart');
Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18475)
+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18476)
@@ -300,15 +300,7 @@
 def BalancethicknessSpcpotentialEnum(): return StringToEnum("BalancethicknessSpcpotential")[0]
 def BalancethicknessApparentMassbalanceEnum(): return StringToEnum("BalancethicknessApparentMassbalance")[0]
 def Balancethickness2MisfitEnum(): return StringToEnum("Balancethickness2Misfit")[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]
+def BalancethicknessDiffusionCoefficientEnum(): return StringToEnum("BalancethicknessDiffusionCoefficient")[0]
 def SurfaceforcingsEnum(): return StringToEnum("Surfaceforcings")[0]
 def SMBEnum(): return StringToEnum("SMB")[0]
 def SurfaceforcingsMassBalanceEnum(): return StringToEnum("SurfaceforcingsMassBalance")[0]
@@ -589,6 +581,7 @@
 def XTaylorHoodEnum(): return StringToEnum("XTaylorHood")[0]
 def OneLayerP4zEnum(): return StringToEnum("OneLayerP4z")[0]
 def CrouzeixRaviartEnum(): return StringToEnum("CrouzeixRaviart")[0]
+def LACrouzeixRaviartEnum(): return StringToEnum("LACrouzeixRaviart")[0]
 def SaveResultsEnum(): return StringToEnum("SaveResults")[0]
 def BoolExternalResultEnum(): return StringToEnum("BoolExternalResult")[0]
 def DoubleExternalResultEnum(): return StringToEnum("DoubleExternalResult")[0]
