Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18060)
+++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 18061)
@@ -124,5 +124,86 @@
 }/*}}}*/
 void AdjointBalancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
-	_error_("Not implemented yet");
+	/*The gradient of the cost function is calculated in 2 parts.
+	 *
+	 * dJ    \partial J   \partial lambda^T(KU-F)
+	 * --  = ---------- + ------------------------
+	 * dk    \partial k   \parial k                  
+	 *
+	 * */
+
+	/*If on water, grad = 0: */
+	if(!element->IsIceInElement()) return;
+
+	/*Get list of cost functions*/
+	int *responses = NULL;
+	int num_responses,resp;
+	element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+	element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
+
+	/*Check that control_type is supported*/
+	if(control_type!=BalancethicknessApparentMassbalanceEnum){
+		_error_("Control "<<EnumToStringx(control_type)<<" not supported");
+	}
+
+	/*Deal with first part (partial derivative a J with respect to k)*/
+	for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
+		case Balancethickness2MisfitEnum: /*Nothing, \partial J/\partial k = 0*/ 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;
+		default: _error_("control type not supported yet: " << EnumToStringx(control_type));
+	}
+
+	/*Clean up and return*/
+	xDelete<int>(responses);
+
+}/*}}}*/
+void AdjointBalancethickness2Analysis::GradientJAdot(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble Jdet,weight;
+	IssmDouble lambda; 
+	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* adjoint_input = element->GetInput(AdjointEnum);                            _assert_(adjoint_input);
+	Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_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);
+		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;
+			_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){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h	(revision 18060)
+++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h	(revision 18061)
@@ -28,4 +28,5 @@
 		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 InputUpdateFromSolution(IssmDouble* solution,Element* element);
 		void UpdateConstraints(FemModel* femmodel);
Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 18060)
+++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 18061)
@@ -1035,4 +1035,83 @@
 		dragcoefficient_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
 
+		/*Build gradient vector (actually -dJ/ddrag): */
+		for(int i=0;i<numvertices;i++){
+			if(dim==2){
+				ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
+			}
+			else{
+				ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0];
+			}
+			_assert_(!xIsNan<IssmDouble>(ge[i]));
+		}
+	}
+	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(dbasis);
+	xDelete<IssmDouble>(ge);
+	xDelete<int>(vertexpidlist);
+	delete gauss;
+	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+
+}/*}}}*/
+void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+	/*Intermediaries*/
+	int      domaintype,dim;
+	Element* basalelement;
+
+	/*Get basal element*/
+	element->FindParam(&domaintype,DomainTypeEnum);
+	switch(domaintype){
+		case Domain2DhorizontalEnum:
+			basalelement = element;
+			dim          = 2;
+			break;
+		case Domain2DverticalEnum:
+			if(!element->IsOnBase()) return;
+			basalelement = element->SpawnBasalElement();
+			dim          = 1;
+			break;
+		case Domain3DEnum:
+			if(!element->IsOnBase()) return;
+			basalelement = element->SpawnBasalElement();
+			dim          = 2;
+			break;
+		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+	}
+
+	/*Intermediaries*/
+	IssmDouble Jdet,weight;
+	IssmDouble dk[3]; 
+	IssmDouble *xyz_list= NULL;
+
+	/*Fetch number of vertices for this finite element*/
+	int numvertices = basalelement->GetNumberOfVertices();
+
+	/*Initialize some vectors*/
+	IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
+	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
+	int*        vertexpidlist = xNew<int>(numvertices);
+
+	/*Retrieve all inputs we will be needing: */
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	basalelement->GradientIndexing(&vertexpidlist[0],control_index);
+	Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum);              _assert_(rheologyb_input);
+	Input* weights_input   = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=basalelement->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
+		weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum);
+
+		/*Build alpha_complement_list: */
+		rheologyb_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
+
 		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
 		for(int i=0;i<numvertices;i++){
@@ -1055,83 +1134,4 @@
 	delete gauss;
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-
-}/*}}}*/
-void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
-
-	/*Intermediaries*/
-	int      domaintype,dim;
-	Element* basalelement;
-
-	/*Get basal element*/
-	element->FindParam(&domaintype,DomainTypeEnum);
-	switch(domaintype){
-		case Domain2DhorizontalEnum:
-			basalelement = element;
-			dim          = 2;
-			break;
-		case Domain2DverticalEnum:
-			if(!element->IsOnBase()) return;
-			basalelement = element->SpawnBasalElement();
-			dim          = 1;
-			break;
-		case Domain3DEnum:
-			if(!element->IsOnBase()) return;
-			basalelement = element->SpawnBasalElement();
-			dim          = 2;
-			break;
-		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-	}
-
-	/*Intermediaries*/
-	IssmDouble Jdet,weight;
-	IssmDouble dk[3]; 
-	IssmDouble *xyz_list= NULL;
-
-	/*Fetch number of vertices for this finite element*/
-	int numvertices = basalelement->GetNumberOfVertices();
-
-	/*Initialize some vectors*/
-	IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
-	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
-	int*        vertexpidlist = xNew<int>(numvertices);
-
-	/*Retrieve all inputs we will be needing: */
-	basalelement->GetVerticesCoordinates(&xyz_list);
-	basalelement->GradientIndexing(&vertexpidlist[0],control_index);
-	Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum);              _assert_(rheologyb_input);
-	Input* weights_input   = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
-		weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum);
-
-		/*Build alpha_complement_list: */
-		rheologyb_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
-
-		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
-		for(int i=0;i<numvertices;i++){
-			if(dim==2){
-				ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
-			}
-			else{
-				ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0];
-			}
-			_assert_(!xIsNan<IssmDouble>(ge[i]));
-		}
-	}
-	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
-
-	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(dbasis);
-	xDelete<IssmDouble>(ge);
-	xDelete<int>(vertexpidlist);
-	delete gauss;
-	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
 void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/