Index: ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 18833)
+++ ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 18834)
@@ -336,15 +336,15 @@
 	if(!element->IsOnSurface()) return NULL;
 
 	/*Intermediaries */
-	int        num_responses,i,dim;
+	int        num_responses,i,domaintype;
 	IssmDouble Jdet,obs_velocity_mag,velocity_mag;
 	IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight;
 	IssmDouble scalex,scaley,scale,S;
 	int        *responses    = NULL;
 	IssmDouble *xyz_list_top = NULL;
 
-	/*Get problem dimension*/
-	element->FindParam(&dim,DomainDimensionEnum);
+	/* Get domaintype*/
+	element->FindParam(&domaintype,DomainTypeEnum);
 
 	/*Fetch number of nodes and dof for this finite element*/
 	int vnumnodes = element->NumberofNodesVelocity();
@@ -352,7 +352,7 @@
 
 	/*Prepare coordinate system list*/
 	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
-	if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
+	if(domaintype==Domain2DverticalEnum) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
 	else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
 	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
 
@@ -366,9 +366,13 @@
 	element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
 	Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
 	Input* vx_input      = element->GetInput(VxEnum);                                 _assert_(vx_input);
-	Input* vy_input      = element->GetInput(VyEnum);                                 _assert_(vy_input);
 	Input* vxobs_input   = element->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
-	Input* vyobs_input   = element->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
+	Input* vy_input=NULL;
+	Input* vyobs_input=NULL;
+	if (domaintype!=Domain2DverticalEnum) {
+		Input* vy_input      = element->GetInput(VyEnum);                                 _assert_(vy_input);
+		Input* vyobs_input   = element->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
+	}
 	IssmDouble epsvel  = 2.220446049250313e-16;
 	IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/
 
@@ -388,9 +392,11 @@
 		element->NodalFunctionsVelocity(vbasis,gauss);
 
 		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
 		vxobs_input->GetInputValue(&vxobs,gauss);
-		vyobs_input->GetInputValue(&vyobs,gauss);
+		if (domaintype!=Domain2DverticalEnum) {
+			vy_input->GetInputValue(&vy,gauss);
+			vyobs_input->GetInputValue(&vyobs,gauss);
+		}
 
 		/*Loop over all requested responses*/
 		for(int resp=0;resp<num_responses;resp++){
@@ -408,12 +414,16 @@
 					 *        du     obs
 					 */
 					for(i=0;i<vnumnodes;i++){
-						dux=vxobs-vx;
-						pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
-						if(dim==3){
+						if(domaintype!=Domain2DverticalEnum){
+							dux=vxobs-vx;
+							pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
 							duy=vyobs-vy;
-							pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+							pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
 						}
+						else{
+							dux=vxobs-vx;
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
+						}
 					}
 					break;
 				case SurfaceRelVelMisfitEnum:
@@ -429,12 +439,20 @@
 					 *               obs
 					 */
 					for(i=0;i<vnumnodes;i++){
-						scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
-						scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
-						dux=scalex*(vxobs-vx);
-						duy=scaley*(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(domaintype!=Domain2DverticalEnum){
+							scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
+							scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
+							dux=scalex*(vxobs-vx);
+							duy=scaley*(vyobs-vy);
+							pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
+							pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+						}
+						else{
+							scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
+							dux=scalex*(vxobs-vx);
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
+							pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+						}
 					}
 					break;
 				case SurfaceLogVelMisfitEnum:
@@ -450,21 +468,21 @@
 					 *            
 					 */
 					for(i=0;i<vnumnodes;i++){
-						if(dim==3){
+						if(domaintype!=Domain2DverticalEnum){
 							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]; 
+							pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
+							pe->values[i*3+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]; 
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
 						}
 					}
 					break;
@@ -479,11 +497,18 @@
 					 *        du      S  2 sqrt(...)           obs
 					 */
 					for(i=0;i<vnumnodes;i++){
-						scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
-						dux=scale*(vxobs-vx);
-						duy=scale*(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 (domaintype!=Domain2DverticalEnum){
+							scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
+							dux=scale*(vxobs-vx);
+							duy=scale*(vyobs-vy);
+							pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
+							pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+						}
+						else{
+							scale=1./(S*2*fabs(vx-vxobs)+epsvel);
+							dux=scale*(vxobs-vx);
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
+						}
 					}
 					break;
 				case SurfaceLogVxVyMisfitEnum:
@@ -497,10 +522,16 @@
 					 *        du                         |u| + eps  |u|                           u + eps
 					 */
 					for(i=0;i<vnumnodes;i++){
-						dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
-						duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
-						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(domaintype!=Domain2DverticalEnum){
+							dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
+							duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
+							pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
+							pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+						}
+						else{
+							dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
+						}
 					}
 					break;
 				case DragCoefficientAbsGradientEnum:
@@ -541,13 +572,16 @@
 	if(!element->IsOnSurface()) return NULL;
 
 	/*Intermediaries */
-	int        num_responses,i;
+	int        num_responses,i,domaintype;
 	IssmDouble Jdet,obs_velocity_mag,velocity_mag;
 	IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight;
 	IssmDouble scalex,scaley,scale,S;
 	int        *responses    = NULL;
 	IssmDouble *xyz_list_top = NULL;
 
+	/* Get domaintype*/
+	element->FindParam(&domaintype,DomainTypeEnum);
+
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
 
@@ -561,9 +595,13 @@
 	element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
 	Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
 	Input* vx_input      = element->GetInput(VxEnum);                                 _assert_(vx_input);
-	Input* vy_input      = element->GetInput(VyEnum);                                 _assert_(vy_input);
 	Input* vxobs_input   = element->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
-	Input* vyobs_input   = element->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
+	Input* vy_input=NULL;
+	Input* vyobs_input=NULL;
+	if(domaintype!=Domain2DverticalEnum){
+		vy_input      = element->GetInput(VyEnum);                                 _assert_(vy_input);
+		vyobs_input   = element->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
+	}
 	IssmDouble epsvel  = 2.220446049250313e-16;
 	IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/
 
@@ -583,10 +621,11 @@
 		element->NodalFunctions(basis, gauss);
 
 		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
 		vxobs_input->GetInputValue(&vxobs,gauss);
-		vyobs_input->GetInputValue(&vyobs,gauss);
-
+		if(domaintype!=Domain2DverticalEnum){
+			vy_input->GetInputValue(&vy,gauss);
+			vyobs_input->GetInputValue(&vyobs,gauss);
+		}
 		/*Loop over all requested responses*/
 		for(int resp=0;resp<num_responses;resp++){
 			weights_input->GetInputValue(&weight,gauss,responses[resp]);
@@ -603,10 +642,16 @@
 					 *        du     obs
 					 */
 					for(i=0;i<numnodes;i++){
-						dux=vxobs-vx;
-						duy=vyobs-vy;
-						pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
-						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						if(domaintype!=Domain2DverticalEnum){
+							dux=vxobs-vx;
+							duy=vyobs-vy;
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+							pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						}
+						else{ 
+							dux=vxobs-vx;
+							pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						}
 					}
 					break;
 				case SurfaceRelVelMisfitEnum:
@@ -622,12 +667,19 @@
 					 *               obs
 					 */
 					for(i=0;i<numnodes;i++){
-						scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
-						scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
-						dux=scalex*(vxobs-vx);
-						duy=scaley*(vyobs-vy);
-						pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
-						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						if(domaintype!=Domain2DverticalEnum){
+							scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
+							scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
+							dux=scalex*(vxobs-vx);
+							duy=scaley*(vyobs-vy);
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+							pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						}
+						else{
+							scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
+							dux=scalex*(vxobs-vx);
+							pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						}
 					}
 					break;
 				case SurfaceLogVelMisfitEnum:
@@ -643,13 +695,22 @@
 					 *            
 					 */
 					for(i=0;i<numnodes;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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
-						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						if(domaintype!=Domain2DverticalEnum){
+							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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+							pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						}
+						else{
+							velocity_mag    =fabs(vx)+epsvel;
+							obs_velocity_mag=fabs(vxobs)+epsvel;
+							scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
+							dux=scale*vx;
+							pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						}
 					}
 					break;
 				case SurfaceAverageVelMisfitEnum:
@@ -663,11 +724,18 @@
 					 *        du      S  2 sqrt(...)           obs
 					 */
 					for(i=0;i<numnodes;i++){
-						scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
-						dux=scale*(vxobs-vx);
-						duy=scale*(vyobs-vy);
-						pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
-						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						if(domaintype!=Domain2DverticalEnum){
+							scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
+							dux=scale*(vxobs-vx);
+							duy=scale*(vyobs-vy);
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+							pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						}
+						else{
+							scale=1./(S*2*fabs(vx-vxobs)+epsvel);
+							dux=scale*(vxobs-vx);
+							pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						}
 					}
 					break;
 				case SurfaceLogVxVyMisfitEnum:
@@ -681,10 +749,16 @@
 					 *        du                         |u| + eps  |u|                           u + eps
 					 */
 					for(i=0;i<numnodes;i++){
-						dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
-						duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
-						pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
-						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						if(domaintype!=Domain2DverticalEnum){
+							dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
+							duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+							pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						}
+						else{
+							dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
+							pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						}
 					}
 					break;
 				case DragCoefficientAbsGradientEnum:
@@ -702,6 +776,9 @@
 				case RheologyBbarAbsGradientEnum:
 					/*Nothing in P vector*/
 					break;
+				case RheologyBAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
 				default:
 					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
 			}
@@ -709,7 +786,7 @@
 	}
 
 	/*Transform coordinate system*/
-	element->TransformLoadVectorCoord(pe,XYEnum);
+	if(domaintype!=Domain2DverticalEnum) element->TransformLoadVectorCoord(pe,XYEnum);
 
 	/*Clean up and return*/
 	xDelete<int>(responses);
@@ -731,6 +808,10 @@
 		case Domain2DhorizontalEnum:
 			basalelement = element;
 			break;
+		case Domain2DverticalEnum:
+			if(!element->IsOnBase()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
 		case Domain3DEnum:
 			if(!element->IsOnBase()) return NULL;
 			basalelement = element->SpawnBasalElement();
@@ -759,9 +840,13 @@
 	basalelement->FindParam(&responses,NULL,InversionCostFunctionsEnum);
 	Input* weights_input = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
 	Input* vx_input      = basalelement->GetInput(VxEnum);                                 _assert_(vx_input);
-	Input* vy_input      = basalelement->GetInput(VyEnum);                                 _assert_(vy_input);
 	Input* vxobs_input   = basalelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
-	Input* vyobs_input   = basalelement->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
+	Input* vy_input=NULL;
+	Input* vyobs_input=NULL;
+	if(domaintype!=Domain2DverticalEnum){
+		vy_input      = basalelement->GetInput(VyEnum);                                 _assert_(vy_input);
+		vyobs_input   = basalelement->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
+	}
 	IssmDouble epsvel  = 2.220446049250313e-16;
 	IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/
 
@@ -781,10 +866,11 @@
 		basalelement->NodalFunctions(basis, gauss);
 
 		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
 		vxobs_input->GetInputValue(&vxobs,gauss);
-		vyobs_input->GetInputValue(&vyobs,gauss);
-
+		if(domaintype!=Domain2DverticalEnum){
+			vy_input->GetInputValue(&vy,gauss);
+			vyobs_input->GetInputValue(&vyobs,gauss);
+		}
 		/*Loop over all requested responses*/
 		for(int resp=0;resp<num_responses;resp++){
 			weights_input->GetInputValue(&weight,gauss,responses[resp]);
@@ -801,10 +887,16 @@
 					 *        du     obs
 					 */
 					for(i=0;i<numnodes;i++){
-						dux=vxobs-vx;
-						duy=vyobs-vy;
-						pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
-						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						if(domaintype!=Domain2DverticalEnum){
+							dux=vxobs-vx;
+							duy=vyobs-vy;
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+							pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						}
+						else {
+							dux=vxobs-vx;
+							pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						}
 					}
 					break;
 				case SurfaceRelVelMisfitEnum:
@@ -820,12 +912,19 @@
 					 *               obs
 					 */
 					for(i=0;i<numnodes;i++){
-						scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
-						scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
-						dux=scalex*(vxobs-vx);
-						duy=scaley*(vyobs-vy);
-						pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
-						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						if(domaintype!=Domain2DverticalEnum){
+							scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
+							scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
+							dux=scalex*(vxobs-vx);
+							duy=scaley*(vyobs-vy);
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+							pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						}
+						else{
+							scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
+							dux=scalex*(vxobs-vx);
+							pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						}
 					}
 					break;
 				case SurfaceLogVelMisfitEnum:
@@ -841,13 +940,22 @@
 					 *            
 					 */
 					for(i=0;i<numnodes;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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
-						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						if(domaintype!=Domain2DverticalEnum){
+							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*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+							pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						}
+						else{
+							velocity_mag    =fabs(vx)+epsvel;
+							obs_velocity_mag=fabs(vxobs)+epsvel;
+							scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
+							dux=scale*vx;
+							pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						}
 					}
 					break;
 				case SurfaceAverageVelMisfitEnum:
@@ -861,11 +969,18 @@
 					 *        du      S  2 sqrt(...)           obs
 					 */
 					for(i=0;i<numnodes;i++){
-						scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
-						dux=scale*(vxobs-vx);
-						duy=scale*(vyobs-vy);
-						pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
-						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						if(domaintype!=Domain2DverticalEnum){
+							scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
+							dux=scale*(vxobs-vx);
+							duy=scale*(vyobs-vy);
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+							pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						}
+						else{
+							scale=1./(S*2*fabs(vx-vxobs)+epsvel);
+							dux=scale*(vxobs-vx);
+							pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						}
 					}
 					break;
 				case SurfaceLogVxVyMisfitEnum:
@@ -879,10 +994,16 @@
 					 *        du                         |u| + eps  |u|                           u + eps
 					 */
 					for(i=0;i<numnodes;i++){
-						dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
-						duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
-						pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
-						pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						if(domaintype!=Domain2DverticalEnum){
+							dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
+							duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
+							pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+							pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
+						}
+						else{
+							dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
+							pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						}
 					}
 					break;
 				case DragCoefficientAbsGradientEnum:
@@ -900,6 +1021,9 @@
 				case RheologyBbarAbsGradientEnum:
 					/*Nothing in P vector*/
 					break;
+				case RheologyBAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
 				default:
 					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
 			}
@@ -907,7 +1031,7 @@
 	}
 
 	/*Transform coordinate system*/
-	basalelement->TransformLoadVectorCoord(pe,XYEnum);
+	if(domaintype!=Domain2DverticalEnum)	basalelement->TransformLoadVectorCoord(pe,XYEnum);
 
 	/*Clean up and return*/
 	xDelete<int>(responses);
@@ -950,7 +1074,8 @@
 	/*Check that control_type is supported*/
 	if(control_type!=MaterialsRheologyBbarEnum && 
 		control_type!=FrictionCoefficientEnum   && 
-		control_type!=DamageDbarEnum){
+		control_type!=DamageDbarEnum            &&
+		control_type!=MaterialsRheologyBEnum){
 		_error_("Control "<<EnumToStringx(control_type)<<" not supported");
 	}
 
@@ -962,7 +1087,8 @@
 		case SurfaceLogVxVyMisfitEnum:    /*Nothing, \partial J/\partial k = 0*/ break;
 		case SurfaceAverageVelMisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
 		case DragCoefficientAbsGradientEnum: GradientJDragGradient(element,gradient,control_index); break;
-		case RheologyBbarAbsGradientEnum:    GradientJBGradient(element,gradient,control_index);    break;
+		case RheologyBbarAbsGradientEnum:    GradientJBbarGradient(element,gradient,control_index); break;
+		case RheologyBAbsGradientEnum:       GradientJBGradient(element,gradient,control_index);    break;
 		default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
 	}
 
@@ -988,6 +1114,15 @@
 				default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
 			}
 			break;
+		case MaterialsRheologyBEnum:
+			switch(approximation){
+				case SSAApproximationEnum: GradientJBSSA(element,gradient,control_index); break;
+				case HOApproximationEnum:  GradientJBHO( element,gradient,control_index); break;
+				case FSApproximationEnum:  GradientJBFS( element,gradient,control_index); break;
+				case NoneApproximationEnum: /*Gradient is 0*/                    break;
+				default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
+			}
+			break;
 		case DamageDbarEnum:
 			switch(approximation){
 				case SSAApproximationEnum: GradientJDSSA(element,gradient,control_index); break;
@@ -1084,10 +1219,10 @@
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 
 }/*}}}*/
-void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+void AdjointHorizAnalysis::GradientJBbarGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 
 	/*Intermediaries*/
-	int      domaintype,dim;
+	int      domaintype;
 	Element* basalelement;
 
 	/*Get basal element*/
@@ -1095,17 +1230,14 @@
 	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");
 	}
@@ -1143,7 +1275,7 @@
 
 		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
 		for(int i=0;i<numvertices;i++){
-			if(dim==2){
+			if(domaintype!=Domain2DverticalEnum){
 				ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
 			}
 			else{
@@ -1162,6 +1294,72 @@
 	delete gauss;
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
+void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+
+	/*Intermediaries*/
+	int      domaintype;
+
+	/*Get basal element*/
+	element->FindParam(&domaintype,DomainTypeEnum);
+	switch(domaintype){
+		case Domain2DhorizontalEnum:
+			break;
+		case Domain2DverticalEnum:
+			break;
+		case Domain3DEnum:
+			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 = element->GetNumberOfVertices();
+
+	/*Initialize some vectors*/
+	IssmDouble* dbasis        = xNew<IssmDouble>(3*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* rheology_input = element->GetInput(MaterialsRheologyBEnum);              _assert_(rheology_input);
+	Input* weights_input   = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_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->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
+		weights_input->GetInputValue(&weight,gauss,RheologyBAbsGradientEnum);
+
+		/*Build alpha_complement_list: */
+		rheology_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
+
+		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
+		for(int i=0;i<numvertices;i++){
+			if(domaintype!=Domain2DverticalEnum){
+				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]+dbasis[1*numvertices+i]*dk[1]+dbasis[2*numvertices+i]*dk[2]);
+			}
+			_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;
+}/*}}}*/
 void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 
 	/*return if floating (gradient is 0)*/
@@ -1262,12 +1460,13 @@
 	if(!element->IsOnBase()) return;
 
 	/*Intermediaries*/
-	int        dim=3;
 	IssmDouble Jdet,weight;
 	IssmDouble drag,dalpha2dk;
 	IssmDouble vx,vy,lambda,mu;
 	IssmDouble *xyz_list_base= NULL;
 
+	int      domaintype,dim;
+	element->FindParam(&domaintype,DomainTypeEnum);
 	/*Fetch number of vertices for this finite element*/
 	int numvertices = element->GetNumberOfVertices();
 
@@ -1277,26 +1476,33 @@
 	int*        vertexpidlist = xNew<int>(numvertices);
 
 	/*Build friction element, needed later: */
+	if(domaintype!=Domain2DverticalEnum) dim=3;
+	else dim=2;
 	Friction* friction=new Friction(element,dim);
 
 	/*Retrieve all inputs we will be needing: */
 	element->GetVerticesCoordinatesBase(&xyz_list_base);
 	element->GradientIndexing(&vertexpidlist[0],control_index);
 	Input* vx_input        = element->GetInput(VxEnum);                   _assert_(vx_input);
-	Input* vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
+	Input* vy_input        = NULL;
 	Input* adjointx_input  = element->GetInput(AdjointxEnum);             _assert_(adjointx_input);
-	Input* adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
+	Input* adjointy_input  = NULL;
 	Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum);  _assert_(dragcoeff_input);
-
+	if(domaintype!=Domain2DverticalEnum){
+		vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
+		adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
+	}
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGaussBase(4);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
 
 		adjointx_input->GetInputValue(&lambda, gauss);
-		adjointy_input->GetInputValue(&mu, gauss);
 		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
+		if(domaintype!=Domain2DverticalEnum){
+			adjointy_input->GetInputValue(&mu, gauss);
+			vy_input->GetInputValue(&vy,gauss);
+		}
 		dragcoeff_input->GetInputValue(&drag, gauss);
 
 		friction->GetAlphaComplement(&dalpha2dk,gauss);
@@ -1306,7 +1512,8 @@
 
 		/*Build gradient vector (actually -dJ/dD): */
 		for(int i=0;i<numvertices;i++){
-			ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
+			if(domaintype!=Domain2DverticalEnum) ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
+			else ge[i]+=-2.*drag*dalpha2dk*(lambda*vx)*Jdet*gauss->weight*basis[i];
 			_assert_(!xIsNan<IssmDouble>(ge[i]));
 		}
 	}
@@ -1327,7 +1534,7 @@
 	if(!element->IsOnBase()) return;
 
 	/*Intermediaries*/
-	int        dim=3;
+	int        domaintype,dim;
 	IssmDouble Jdet,weight;
 	IssmDouble drag,dalpha2dk,normal[3];
 	IssmDouble vx,vy,vz,lambda,mu,xi;
@@ -1341,7 +1548,12 @@
 	IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
 	int*        vertexpidlist = xNew<int>(numvertices);
 
+	/* get domaintype */
+	element->FindParam(&domaintype,DomainTypeEnum);
+
 	/*Build friction element, needed later: */
+	if(domaintype!=Domain2DverticalEnum) dim=3;
+	else dim=2;
 	Friction* friction=new Friction(element,dim);
 
 	/*Retrieve all inputs we will be needing: */
@@ -1349,10 +1561,14 @@
 	element->GradientIndexing(&vertexpidlist[0],control_index);
 	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_(vy_input);
 	Input* adjointx_input  = element->GetInput(AdjointxEnum);             _assert_(adjointx_input);
 	Input* adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
-	Input* adjointz_input  = element->GetInput(AdjointzEnum);             _assert_(adjointz_input);
+	Input* vz_input        = NULL;
+	Input* adjointz_input  = NULL;
+	if(domaintype!=Domain2DverticalEnum){
+		Input* vz_input        = element->GetInput(VzEnum);                   _assert_(vy_input);
+		Input* adjointz_input  = element->GetInput(AdjointzEnum);             _assert_(adjointz_input);
+	}
 	Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum);  _assert_(dragcoeff_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -1362,10 +1578,12 @@
 
 		adjointx_input->GetInputValue(&lambda, gauss);
 		adjointy_input->GetInputValue(&mu, gauss);
-		adjointz_input->GetInputValue(&xi    ,gauss);
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
-		vz_input->GetInputValue(&vz,gauss);
+		if(domaintype!=Domain2DverticalEnum){
+			adjointz_input->GetInputValue(&xi    ,gauss);
+			vz_input->GetInputValue(&vz,gauss);
+		}
 		dragcoeff_input->GetInputValue(&drag, gauss);
 
 		friction->GetAlphaComplement(&dalpha2dk,gauss);
@@ -1375,14 +1593,25 @@
 		element->NodalFunctionsP1(basis,gauss);
 
 		/*Build gradient vector (actually -dJ/dk): */
-		for(int i=0;i<numvertices;i++){
-			ge[i]+=(
-						-lambda*(2*drag*dalpha2dk*(vx - vz*normal[0]*normal[2]))
-						-mu    *(2*drag*dalpha2dk*(vy - vz*normal[1]*normal[2]))
-						-xi    *(2*drag*dalpha2dk*(-vx*normal[0]*normal[2]-vy*normal[1]*normal[2]))
-						)*Jdet*gauss->weight*basis[i];
-			_assert_(!xIsNan<IssmDouble>(ge[i]));
+		if(domaintype!=Domain2DverticalEnum){
+			for(int i=0;i<numvertices;i++){
+				ge[i]+=(
+							-lambda*(2*drag*dalpha2dk*(vx - vz*normal[0]*normal[2]))
+							-mu    *(2*drag*dalpha2dk*(vy - vz*normal[1]*normal[2]))
+							-xi    *(2*drag*dalpha2dk*(-vx*normal[0]*normal[2]-vy*normal[1]*normal[2]))
+						 )*Jdet*gauss->weight*basis[i];
+				_assert_(!xIsNan<IssmDouble>(ge[i]));
+			}
 		}
+		else{
+			for(int i=0;i<numvertices;i++){
+				ge[i]+=(
+							-lambda*2*drag*dalpha2dk*vx
+							-mu    *2*drag*dalpha2dk*vy
+						 )*Jdet*gauss->weight*basis[i];
+				_assert_(!xIsNan<IssmDouble>(ge[i]));
+			}
+		}
 	}
 	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
 
@@ -1492,6 +1721,170 @@
 	/*WARNING: We use SSA as an estimate for now*/
 	this->GradientJBbarSSA(element,gradient,control_index);
 }/*}}}*/
+void AdjointHorizAnalysis::GradientJBSSA(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 thickness,dmudB;
+	IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3]; 
+	IssmDouble *xyz_list= NULL;
+
+	/*Fetch number of vertices for this finite element*/
+	int numvertices = basalelement->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: */
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	basalelement->GradientIndexing(&vertexpidlist[0],control_index);
+	Input* thickness_input = basalelement->GetInput(ThicknessEnum);             _assert_(thickness_input);
+	Input* vx_input        = basalelement->GetInput(VxEnum);                    _assert_(vx_input);
+	Input* vy_input        = basalelement->GetInput(VyEnum);                    _assert_(vy_input);
+	Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);              _assert_(adjointx_input);
+	Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);              _assert_(adjointy_input);
+	Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBEnum); _assert_(rheologyb_input);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=basalelement->NewGauss(4);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		thickness_input->GetInputValue(&thickness,gauss);
+		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+		adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
+		adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss);
+
+		basalelement->dViscositydBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input);
+
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctionsP1(basis,gauss);
+
+		/*Build gradient vector (actually -dJ/dB): */
+		for(int i=0;i<numvertices;i++){
+			ge[i]+=-dmudB*thickness*(
+						(2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
+						)*Jdet*gauss->weight*basis[i];
+			_assert_(!xIsNan<IssmDouble>(ge[i]));
+		}
+	}
+	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(ge);
+	xDelete<int>(vertexpidlist);
+	delete gauss;
+	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+}/*}}}*/
+void AdjointHorizAnalysis::GradientJBHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+	/*Intermediaries*/
+	int      domaintype,dim;
+
+	/*Get domaintype*/
+	element->FindParam(&domaintype,DomainTypeEnum);
+
+	/*Intermediaries*/
+	IssmDouble Jdet,weight;
+	IssmDouble thickness,dmudB;
+	IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3]; 
+	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* thickness_input = element->GetInput(ThicknessEnum);             _assert_(thickness_input);
+	Input* vx_input        = element->GetInput(VxEnum);                    _assert_(vx_input);
+	Input* vy_input        = NULL;
+	Input* adjointx_input  = element->GetInput(AdjointxEnum);              _assert_(adjointx_input);
+	Input* adjointy_input  = NULL;
+	Input* rheologyb_input = element->GetInput(MaterialsRheologyBEnum); _assert_(rheologyb_input);
+	if(domaintype!=Domain2DverticalEnum){
+		vy_input        = element->GetInput(VyEnum);                   _assert_(vy_input);
+		adjointy_input  = element->GetInput(AdjointyEnum);             _assert_(adjointy_input);
+	}
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(4);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		thickness_input->GetInputValue(&thickness,gauss);
+		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+		adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
+		dim=2;
+		if(domaintype!=Domain2DverticalEnum){
+			adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list, gauss);
+			vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+			dim=3;
+		}
+
+		element->dViscositydBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctionsP1(basis,gauss);
+
+		/*Build gradient vector (actually -dJ/dB): */
+		for(int i=0;i<numvertices;i++){
+			if(domaintype!=Domain2DverticalEnum){
+				ge[i]+=-dmudB*thickness*(
+							(2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
+							)*Jdet*gauss->weight*basis[i];
+			}
+			else{
+				ge[i]+=-dmudB*thickness*4*dvx[0]*dadjx[0]*Jdet*gauss->weight*basis[i];
+			}
+			_assert_(!xIsNan<IssmDouble>(ge[i]));
+		}
+	}
+	gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(ge);
+	xDelete<int>(vertexpidlist);
+	delete gauss;
+}/*}}}*/
+void AdjointHorizAnalysis::GradientJBFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
+	/*WARNING: We use HO as an estimate for now*/
+	this->GradientJBHO(element,gradient,control_index);
+}/*}}}*/
 void AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
 
 	/*Intermediaries*/
@@ -1589,10 +1982,14 @@
 	int  i;
 	int* doflist=NULL;
 
+	int    domaintype;
+	element->FindParam(&domaintype,DomainTypeEnum);
+
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
-	int numdof   = numnodes*2;
-
+	int numdof;
+	if(domaintype!=Domain2DverticalEnum)  numdof   = numnodes*2;
+	else			                          numdof   = numnodes*1;
 	/*Fetch dof list and allocate solution vectors*/
 	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 	IssmDouble* values  = xNew<IssmDouble>(numdof);
@@ -1603,16 +2000,18 @@
 	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
 
 	/*Transform solution in Cartesian Space*/
-	element->TransformSolutionCoord(&values[0],XYEnum);
+	if(domaintype!=Domain2DverticalEnum)	element->TransformSolutionCoord(&values[0],XYEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	for(i=0;i<numnodes;i++){
-		lambdax[i]=values[i*NDOF2+0];
-		lambday[i]=values[i*NDOF2+1];
-
+		if(domaintype!=Domain2DverticalEnum){
+			lambdax[i]=values[i*NDOF2+0];
+			lambday[i]=values[i*NDOF2+1];
+		}
+		else {lambdax[i]=values[i];lambday[i]=0;}
 		/*Check solution*/
 		if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
-		if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
+		if(domaintype!=Domain2DverticalEnum && xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
 	}
 
 	/*Add vx and vy as inputs to the tria element: */
@@ -1626,13 +2025,27 @@
 	xDelete<int>(doflist);
 }/*}}}*/
 void AdjointHorizAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
-	int          i,dim;
+	int          i,fe_FS;
 	int*         vdoflist=NULL;
 	int*         pdoflist=NULL;
 	IssmDouble   FSreconditioning;
 
-	element->FindParam(&dim,DomainDimensionEnum);
+	int      domaintype,dim;
+	element->FindParam(&domaintype,DomainTypeEnum);
+	switch(domaintype){
+		case Domain2DhorizontalEnum:
+			dim          = 3;
+			break;
+		case Domain2DverticalEnum:
+			dim          = 2;
+			break;
+		case Domain3DEnum:
+			dim          = 3;
+			break;
+		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+	}
 
+
 	/*Fetch number of nodes and dof for this finite element*/
 	int vnumnodes = element->NumberofNodesVelocity();
 	int pnumnodes = element->NumberofNodesPressure();
@@ -1664,9 +2077,11 @@
 
 	/*fill in all arrays: */
 	for(i=0;i<vnumnodes;i++){
-		lambdax[i] = values[i*NDOF3+0]; if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
-		lambday[i] = values[i*NDOF3+1]; if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
-		lambdaz[i] = values[i*NDOF3+2]; if(xIsNan<IssmDouble>(lambdaz[i])) _error_("NaN found in solution vector");
+		lambdax[i] = values[i*dim+0]; if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
+		lambday[i] = values[i*dim+1]; if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
+		if(dim==3){
+			lambdaz[i] = values[i*dim+2]; if(xIsNan<IssmDouble>(lambdaz[i])) _error_("NaN found in solution vector");
+		}
 	}
 	for(i=0;i<pnumnodes;i++){
 		lambdap[i] = values[vnumdof+i]; if(xIsNan<IssmDouble>(lambdap[i])) _error_("NaN found in solution vector");
@@ -1679,9 +2094,12 @@
 	/*Add vx and vy as inputs to the tria element: */
 	element->AddInput(AdjointxEnum,lambdax,element->VelocityInterpolation());
 	element->AddInput(AdjointyEnum,lambday,element->VelocityInterpolation());
-	element->AddInput(AdjointzEnum,lambdaz,element->VelocityInterpolation());
-	element->AddInput(AdjointpEnum,lambdap,element->PressureInterpolation());
+	if(domaintype!=Domain2DverticalEnum) element->AddInput(AdjointzEnum,lambdaz,element->VelocityInterpolation());
 
+	element->FindParam(&fe_FS,FlowequationFeFSEnum);
+	if(fe_FS!=LATaylorHoodEnum && fe_FS!=LACrouzeixRaviartEnum)	
+	 element->AddInput(AdjointpEnum,lambdap,element->PressureInterpolation());	
+
 	/*Free ressources:*/
 	xDelete<int>(vdoflist);
 	xDelete<int>(pdoflist);
