Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 17968)
+++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 17969)
@@ -208,4 +208,5 @@
 	/*Intermediaries */
 	bool        incomplete_adjoint;
+	int         dim,epssize;
 	IssmDouble  Jdet,mu_prime;
 	IssmDouble  eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij,eps3dotdphii,eps3dotdphij;
@@ -213,8 +214,13 @@
 	IssmDouble *xyz_list = NULL;
 
+	/*Get problem dimension*/
+	element->FindParam(&dim,DomainDimensionEnum);
+	if(dim==2) epssize = 3;
+	else       epssize = 6;
+
 	/*Fetch number of nodes and dof for this finite element*/
 	int vnumnodes = element->NumberofNodesVelocity();
 	int pnumnodes = element->NumberofNodesPressure();
-	int numdof    = vnumnodes*3 + pnumnodes;
+	int numdof    = vnumnodes*dim + pnumnodes;
 
 	/*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
@@ -229,8 +235,14 @@
 	Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input);
 	Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input);
-	Input* vz_input = element->GetInput(VzEnum);_assert_(vz_input);
+	Input* vz_input = NULL;
+	if(dim==3){
+		vz_input = element->GetInput(VzEnum);
+	}
+	else{
+		_error_("Not implemented yet");
+	}
 
 	/*Allocate dbasis*/
-	IssmDouble* dbasis = xNew<IssmDouble>(3*vnumnodes);
+	IssmDouble* dbasis = xNew<IssmDouble>(dim*vnumnodes);
 
 	/* Start  looping on the number of gaussian points: */
@@ -376,7 +388,9 @@
 					for(i=0;i<vnumnodes;i++){
 						dux=vxobs-vx;
-						duy=vyobs-vy;
 						pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
-						pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+						if(dim==3){
+							duy=vyobs-vy;
+							pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+						}
 					}
 					break;
@@ -415,11 +429,20 @@
 					 */
 					for(i=0;i<vnumnodes;i++){
-						velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
-						obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
-						scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
-						dux=scale*vx;
-						duy=scale*vy;
-						pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
-						pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+						if(dim==3){
+							velocity_mag    =sqrt(vx*vx+vy*vy)+epsvel;
+							obs_velocity_mag=sqrt(vxobs*vxobs+vyobs*vyobs)+epsvel;
+							scale=-8.*meanvel*meanvel/(velocity_mag*velocity_mag)*log(velocity_mag/obs_velocity_mag);
+							dux=scale*vx;
+							duy=scale*vy;
+							pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
+							pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+						}
+						else{
+							velocity_mag    =fabs(vx)+epsvel;
+							obs_velocity_mag=fabs(vxobs)+epsvel;
+							scale=-8.*meanvel*meanvel/(velocity_mag*velocity_mag)*log(velocity_mag/obs_velocity_mag);
+							dux=scale*vx;
+							pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
+						}
 					}
 					break;
Index: /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp	(revision 17968)
+++ /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp	(revision 17969)
@@ -21,5 +21,5 @@
 	for (i=0;i<elements->Size();i++){
 		element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-		J+=element->SurfaceAbsVelMisfit();
+		J+=SurfaceAbsVelMisfit(element);
 	}
 
@@ -32,2 +32,90 @@
 	*pJ=J;
 }
+
+IssmDouble SurfaceAbsVelMisfit(Element* element){
+
+	int        domaintype,numcomponents;
+	bool       surfaceintegral;
+	IssmDouble Jelem=0.;
+	IssmDouble misfit,Jdet;
+	IssmDouble vx,vy,vxobs,vyobs,weight;
+	IssmDouble* xyz_list_top = NULL;
+
+	/*Get basal element*/
+	if(!element->IsOnSurface()) return 0.;
+
+	/*If on water, return 0: */
+	if(!element->IsIceInElement()) return 0.;
+
+	/*Get problem dimension*/
+	element->FindParam(&domaintype,DomainTypeEnum);
+	switch(domaintype){
+		case Domain2DverticalEnum:
+			surfaceintegral = true;
+			numcomponents   = 1;
+			break;
+		case Domain3DEnum:
+			surfaceintegral = true;
+			numcomponents   = 2;
+			break;
+		case Domain2DhorizontalEnum:
+			surfaceintegral = false;
+			numcomponents   = 2;
+			break;
+		default: _error_("not supported yet");
+	}
+
+	/* Get node coordinates*/
+	if(surfaceintegral) element->GetVerticesCoordinatesTop(&xyz_list_top);
+	else                element->GetVerticesCoordinates(&xyz_list_top);
+
+	/*Retrieve all inputs we will be needing: */
+	Input* weights_input=element->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
+	Input* vx_input     =element->GetInput(VxEnum);        _assert_(vx_input);
+	Input* vxobs_input  =element->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
+	Input* vy_input     = NULL;
+	Input* vyobs_input  = NULL;
+	if(numcomponents==2){
+		vy_input    =element->GetInput(VyEnum);              _assert_(vy_input);
+		vyobs_input =element->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=NULL;
+	if(surfaceintegral) gauss=element->NewGaussTop(2);
+	else                gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		/* Get Jacobian determinant: */
+		if(surfaceintegral) element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
+		else                element->JacobianDeterminant(&Jdet,xyz_list_top,gauss);
+
+		/*Get all parameters at gaussian point*/
+		weights_input->GetInputValue(&weight,gauss,SurfaceAbsVelMisfitEnum);
+		vx_input->GetInputValue(&vx,gauss);
+		vxobs_input->GetInputValue(&vxobs,gauss);
+		if(numcomponents==2){
+			vy_input->GetInputValue(&vy,gauss);
+			vyobs_input->GetInputValue(&vyobs,gauss);
+		}
+
+		/*Compute SurfaceAbsVelMisfitEnum:
+		 *
+		 *      1  [           2              2 ]
+		 * J = --- | (u - u   )  +  (v - v   )  |
+		 *      2  [       obs            obs   ]
+		 *
+		 */
+		misfit=0.5*(vx-vxobs)*(vx-vxobs);
+		if(numcomponents==2) misfit+=0.5*(vy-vyobs)*(vy-vyobs);
+
+		/*Add to cost function*/
+		Jelem+=misfit*weight*Jdet*gauss->weight;
+	}
+
+	/*clean up and Return: */
+	delete gauss;
+	return Jelem;
+}
Index: /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h	(revision 17968)
+++ /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h	(revision 17969)
@@ -10,4 +10,5 @@
 /* local prototypes: */
 void SurfaceAbsVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters);
+IssmDouble SurfaceAbsVelMisfit(Element* element);
 
 #endif