Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8610)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8611)
@@ -3849,5 +3849,14 @@
 double Penta::DragCoefficientAbsGradient(bool process_units,int weight_index){
 
-	_error_("Not implemented yet");
+	double J;
+	Tria*  tria=NULL;
+
+	/*If on water, on shelf or not on bed, skip: */
+	if(IsOnWater()|| IsOnShelf() || !IsOnBed()) return 0;
+
+	tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria
+	J=tria->DragCoefficientAbsGradient(process_units,weight_index);
+	delete tria->matice; delete tria;
+	return J;
 }
 /*}}}*/
@@ -6751,5 +6760,14 @@
 double Penta::RheologyBbarAbsGradient(bool process_units,int weight_index){
 
-	_error_("Not implemented yet");
+	double J;
+	Tria*  tria=NULL;
+
+	/*If on water, on shelf or not on bed, skip: */
+	if(IsOnWater() || !IsOnBed()) return 0;
+
+	tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria
+	J=tria->RheologyBbarAbsGradient(process_units,weight_index);
+	delete tria->matice; delete tria;
+	return J;
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8610)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8611)
@@ -1812,4 +1812,13 @@
 						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
 					}
+					break;
+				case DragCoefficientAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case ThicknessAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case RheologyBbarAbsGradientEnum:
+					/*Nothing in P vector*/
 					break;
 				default:
@@ -2396,5 +2405,43 @@
 double Tria::DragCoefficientAbsGradient(bool process_units,int weight_index){
 
-	_error_("Not implemented yet");
+	/* Intermediaries */
+	int        ig;
+	double     Jelem = 0;
+	double     weight;
+	double     Jdet;
+	double     xyz_list[NUMVERTICES][3];
+	double     dp[NDOF2];
+	GaussTria *gauss = NULL;
+
+	/*retrieve parameters and inputs*/
+
+	/*If on water, return 0: */
+	if(IsOnWater()) return 0;
+
+	/*Retrieve all inputs we will be needing: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	Input* weights_input=inputs->GetInput(WeightsEnum);         _assert_(weights_input);
+	Input* drag_input   =inputs->GetInput(DragCoefficientEnum); _assert_(drag_input);
+
+	/* Start looping on the number of gaussian points: */
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+
+		/*Get all parameters at gaussian point*/
+		weights_input->GetParameterValue(&weight,gauss,weight_index);
+		drag_input->GetParameterDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
+
+		/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+		//Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	return Jelem;
 }
 /*}}}*/
@@ -4679,5 +4726,43 @@
 double Tria::RheologyBbarAbsGradient(bool process_units,int weight_index){
 
-	_error_("Not implemented yet");
+	/* Intermediaries */
+	int        ig;
+	double     Jelem = 0;
+	double     weight;
+	double     Jdet;
+	double     xyz_list[NUMVERTICES][3];
+	double     dp[NDOF2];
+	GaussTria *gauss = NULL;
+
+	/*retrieve parameters and inputs*/
+
+	/*If on water, return 0: */
+	if(IsOnWater()) return 0;
+
+	/*Retrieve all inputs we will be needing: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	Input* weights_input  =inputs->GetInput(WeightsEnum);              _assert_(weights_input);
+	Input* rheologyb_input=matice->inputs->GetInput(RheologyBbarEnum); _assert_(rheologyb_input);
+
+	/* Start looping on the number of gaussian points: */
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+
+		/*Get all parameters at gaussian point*/
+		weights_input->GetParameterValue(&weight,gauss,weight_index);
+		rheologyb_input->GetParameterDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
+
+		/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+		//Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	return Jelem;
 }
 /*}}}*/
@@ -5211,12 +5296,47 @@
 double Tria::ThicknessAbsGradient(bool process_units,int weight_index){
 
-	_error_("Not implemented yet");
+	/* Intermediaries */
+	int        ig;
+	double     Jelem = 0;
+	double     weight;
+	double     Jdet;
+	double     xyz_list[NUMVERTICES][3];
+	double     dp[NDOF2];
+	GaussTria *gauss = NULL;
+
+	/*retrieve parameters and inputs*/
+
+	/*If on water, return 0: */
+	if(IsOnWater()) return 0;
+
+	/*Retrieve all inputs we will be needing: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	Input* weights_input  =inputs->GetInput(WeightsEnum);   _assert_(weights_input);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+
+	/* Start looping on the number of gaussian points: */
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+
+		/*Get all parameters at gaussian point*/
+		weights_input->GetParameterValue(&weight,gauss,weight_index);
+		thickness_input->GetParameterDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
+
+		/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
+		//Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	return Jelem;
 }
 /*}}}*/
 /*FUNCTION Tria::ThicknessAbsMisfit {{{1*/
 double Tria::ThicknessAbsMisfit(bool process_units,int weight_index){
-
-	/* Constants */
-	const int    numdof=1*NUMVERTICES;
 
 	/*Intermediaries*/
@@ -5233,8 +5353,6 @@
 	if(IsOnWater())return 0;
 
-	/* Get node coordinates and dof list: */
+	/*Retrieve all inputs we will be needing: */
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-
-	/*Retrieve all inputs we will be needing: */
 	Input* thickness_input   =inputs->GetInput(ThicknessEnum);   _assert_(thickness_input);
 	Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);_assert_(thicknessobs_input);
Index: /issm/trunk/src/m/model/ismodelselfconsistent.m
===================================================================
--- /issm/trunk/src/m/model/ismodelselfconsistent.m	(revision 8610)
+++ /issm/trunk/src/m/model/ismodelselfconsistent.m	(revision 8611)
@@ -199,5 +199,5 @@
 	checkvalues(md,{'cm_responses'},...
 		[SurfaceAbsVelMisfitEnum SurfaceRelVelMisfitEnum SurfaceLogVelMisfitEnum SurfaceLogVxVyMisfitEnum SurfaceAverageVelMisfitEnum...
-		ThicknessAbsMisfitEnum]);
+		ThicknessAbsMisfitEnum DragCoefficientAbsGradientEnum RheologyBbarAbsGradientEnum]);
 
 	%WEIGHTS