Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 16846)
+++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 16847)
@@ -30,5 +30,597 @@
 }/*}}}*/
 ElementVector* AdjointHorizAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+
+	int approximation;
+	element->GetInputValue(&approximation,ApproximationEnum);
+	switch(approximation){
+		case SSAApproximationEnum: 
+			return CreatePVectorSSA(element);
+		case HOApproximationEnum: 
+			return CreatePVectorHO(element);
+		case FSApproximationEnum: 
+			return CreatePVectorFS(element);
+		case NoneApproximationEnum:
+			return NULL;
+		default:
+			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
+	}
+}/*}}}*/
+ElementVector* AdjointHorizAnalysis::CreatePVectorFS(Element* element){/*{{{*/
+
+	/*Nothing to be done if not on surface*/
+	if(!element->IsOnSurface()) return NULL;
+
+	/*Intermediaries */
+	int        num_responses,i,meshtype,dim;
+	IssmDouble thickness,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(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2; break;
+		case Mesh3DEnum:         dim = 3; break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->GetNumberOfNodesVelocity();
+	int pnumnodes = element->GetNumberOfNodesPressure();
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+	/*Initialize Element vector and vectors*/
+	ElementVector* pe     = element->NewElementVector(FSApproximationEnum);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesTop(&xyz_list_top);
+	element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+	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);
+	IssmDouble epsvel  = 2.220446049250313e-16;
+	IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/
+
+	/*Get Surface if required by one response*/
+	for(int resp=0;resp<num_responses;resp++){
+		if(responses[resp]==SurfaceAverageVelMisfitEnum){
+			element->GetInputValue(&S,SurfaceAreaEnum); break;
+		}
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGaussTop(4);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
+		element->NodalFunctionsVelocity(vbasis,gauss);
+
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vxobs_input->GetInputValue(&vxobs,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]);
+
+			switch(responses[resp]){
+				case SurfaceAbsVelMisfitEnum:
+					/*
+					 *      1  [           2              2 ]
+					 * J = --- | (u - u   )  +  (v - v   )  |
+					 *      2  [       obs            obs   ]
+					 *
+					 *        dJ
+					 * DU = - -- = (u   - u )
+					 *        du     obs
+					 */
+					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]; 
+					}
+					break;
+				case SurfaceRelVelMisfitEnum:
+					/*
+					 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
+					 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
+					 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
+					 *              obs                        obs                      
+					 *
+					 *        dJ     \bar{v}^2
+					 * DU = - -- = ------------- (u   - u )
+					 *        du   (u   + eps)^2    obs
+					 *               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]; 
+					}
+					break;
+				case SurfaceLogVelMisfitEnum:
+					/*
+					 *                 [        vel + eps     ] 2
+					 * J = 4 \bar{v}^2 | log ( -----------  ) |  
+					 *                 [       vel   + eps    ]
+					 *                            obs
+					 *
+					 *        dJ                 2 * log(...)
+					 * DU = - -- = - 4 \bar{v}^2 -------------  u
+					 *        du                 vel^2 + eps
+					 *            
+					 */
+					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]; 
+					}
+					break;
+				case SurfaceAverageVelMisfitEnum:
+					/*
+					 *      1                    2              2
+					 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
+					 *      S                obs            obs
+					 *
+					 *        dJ      1       1 
+					 * DU = - -- = - --- ----------- * 2 (u - u   )
+					 *        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]; 
+					}
+					break;
+				case SurfaceLogVxVyMisfitEnum:
+					/*
+					 *      1            [        |u| + eps     2          |v| + eps     2  ]
+					 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
+					 *      2            [       |u    |+ eps              |v    |+ eps     ]
+					 *                              obs                       obs
+					 *        dJ                              1      u                             1
+					 * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
+					 *        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]; 
+					}
+					break;
+				case DragCoefficientAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case ThicknessAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case ThicknessAlongGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case ThicknessAcrossGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case RheologyBbarAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				default:
+					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
+			}
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,cs_list);
+
+	/*Clean up and return*/
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(xyz_list_top);
+	xDelete<IssmDouble>(vbasis);
+	delete gauss;
+	return pe;
+	_error_("S");
+}/*}}}*/
+ElementVector* AdjointHorizAnalysis::CreatePVectorHO(Element* element){/*{{{*/
+
+	/*Nothing to be done if not on surface*/
+	if(!element->IsOnSurface()) return NULL;
+
+	/*Intermediaries */
+	int        num_responses,i;
+	IssmDouble thickness,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;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and vectors*/
+	ElementVector* pe    = element->NewElementVector(HOApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesTop(&xyz_list_top);
+	element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+	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);
+	IssmDouble epsvel  = 2.220446049250313e-16;
+	IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/
+
+	/*Get Surface if required by one response*/
+	for(int resp=0;resp<num_responses;resp++){
+		if(responses[resp]==SurfaceAverageVelMisfitEnum){
+			element->GetInputValue(&S,SurfaceAreaEnum); break;
+		}
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGaussTop(4);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
+		element->NodalFunctions(basis, gauss);
+
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vxobs_input->GetInputValue(&vxobs,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]);
+
+			switch(responses[resp]){
+				case SurfaceAbsVelMisfitEnum:
+					/*
+					 *      1  [           2              2 ]
+					 * J = --- | (u - u   )  +  (v - v   )  |
+					 *      2  [       obs            obs   ]
+					 *
+					 *        dJ
+					 * DU = - -- = (u   - u )
+					 *        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]; 
+					}
+					break;
+				case SurfaceRelVelMisfitEnum:
+					/*
+					 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
+					 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
+					 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
+					 *              obs                        obs                      
+					 *
+					 *        dJ     \bar{v}^2
+					 * DU = - -- = ------------- (u   - u )
+					 *        du   (u   + eps)^2    obs
+					 *               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]; 
+					}
+					break;
+				case SurfaceLogVelMisfitEnum:
+					/*
+					 *                 [        vel + eps     ] 2
+					 * J = 4 \bar{v}^2 | log ( -----------  ) |  
+					 *                 [       vel   + eps    ]
+					 *                            obs
+					 *
+					 *        dJ                 2 * log(...)
+					 * DU = - -- = - 4 \bar{v}^2 -------------  u
+					 *        du                 vel^2 + eps
+					 *            
+					 */
+					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]; 
+					}
+					break;
+				case SurfaceAverageVelMisfitEnum:
+					/*
+					 *      1                    2              2
+					 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
+					 *      S                obs            obs
+					 *
+					 *        dJ      1       1 
+					 * DU = - -- = - --- ----------- * 2 (u - u   )
+					 *        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]; 
+					}
+					break;
+				case SurfaceLogVxVyMisfitEnum:
+					/*
+					 *      1            [        |u| + eps     2          |v| + eps     2  ]
+					 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
+					 *      2            [       |u    |+ eps              |v    |+ eps     ]
+					 *                              obs                       obs
+					 *        dJ                              1      u                             1
+					 * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
+					 *        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]; 
+					}
+					break;
+				case DragCoefficientAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case ThicknessAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case ThicknessAlongGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case ThicknessAcrossGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case RheologyBbarAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				default:
+					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
+			}
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,XYEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list_top);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+
+}/*}}}*/
+ElementVector* AdjointHorizAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int      meshtype;
+	Element* basalelement;
+
+	/*Get basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Intermediaries */
+	int         num_responses,i;
+	IssmDouble  thickness,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  = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = basalelement->GetNumberOfNodes();
+
+	/*Initialize Element vector and vectors*/
+	ElementVector* pe    = basalelement->NewElementVector(SSAApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	basalelement->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+	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);
+	IssmDouble epsvel  = 2.220446049250313e-16;
+	IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/
+
+	/*Get Surface if required by one response*/
+	for(int resp=0;resp<num_responses;resp++){
+		if(responses[resp]==SurfaceAverageVelMisfitEnum){
+			basalelement->GetInputValue(&S,SurfaceAreaEnum); break;
+		}
+	}
+
+	/* 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);
+
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctions(basis, gauss);
+
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vxobs_input->GetInputValue(&vxobs,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]);
+
+			switch(responses[resp]){
+				case SurfaceAbsVelMisfitEnum:
+					/*
+					 *      1  [           2              2 ]
+					 * J = --- | (u - u   )  +  (v - v   )  |
+					 *      2  [       obs            obs   ]
+					 *
+					 *        dJ
+					 * DU = - -- = (u   - u )
+					 *        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]; 
+					}
+					break;
+				case SurfaceRelVelMisfitEnum:
+					/*
+					 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
+					 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
+					 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
+					 *              obs                        obs                      
+					 *
+					 *        dJ     \bar{v}^2
+					 * DU = - -- = ------------- (u   - u )
+					 *        du   (u   + eps)^2    obs
+					 *               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]; 
+					}
+					break;
+				case SurfaceLogVelMisfitEnum:
+					/*
+					 *                 [        vel + eps     ] 2
+					 * J = 4 \bar{v}^2 | log ( -----------  ) |  
+					 *                 [       vel   + eps    ]
+					 *                            obs
+					 *
+					 *        dJ                 2 * log(...)
+					 * DU = - -- = - 4 \bar{v}^2 -------------  u
+					 *        du                 vel^2 + eps
+					 *            
+					 */
+					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]; 
+					}
+					break;
+				case SurfaceAverageVelMisfitEnum:
+					/*
+					 *      1                    2              2
+					 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
+					 *      S                obs            obs
+					 *
+					 *        dJ      1       1 
+					 * DU = - -- = - --- ----------- * 2 (u - u   )
+					 *        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]; 
+					}
+					break;
+				case SurfaceLogVxVyMisfitEnum:
+					/*
+					 *      1            [        |u| + eps     2          |v| + eps     2  ]
+					 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
+					 *      2            [       |u    |+ eps              |v    |+ eps     ]
+					 *                              obs                       obs
+					 *        dJ                              1      u                             1
+					 * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
+					 *        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]; 
+					}
+					break;
+				case DragCoefficientAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case ThicknessAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case ThicknessAlongGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case ThicknessAcrossGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				case RheologyBbarAbsGradientEnum:
+					/*Nothing in P vector*/
+					break;
+				default:
+					_error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
+			}
+		}
+	}
+
+	/*Transform coordinate system*/
+	basalelement->TransformLoadVectorCoord(pe,XYEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	delete gauss;
+	return pe;
 }/*}}}*/
 void AdjointHorizAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h	(revision 16846)
+++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h	(revision 16847)
@@ -23,4 +23,7 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
+		ElementVector* CreatePVectorSSA(Element* element);
+		ElementVector* CreatePVectorHO(Element* element);
+		ElementVector* CreatePVectorFS(Element* element);
 		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16846)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16847)
@@ -55,4 +55,5 @@
 		virtual void   FindParam(int* pvalue,int paramenum)=0;
 		virtual void   FindParam(IssmDouble* pvalue,int paramenum)=0;
+		virtual void   FindParam(int** pvalues,int* psize,int paramenum)=0;
 		virtual int    FiniteElement(void)=0;
 		virtual IssmDouble GetMaterialParameter(int enum_in)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16846)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16847)
@@ -869,4 +869,9 @@
 void Penta::FindParam(IssmDouble* pvalue,int paramenum){
 	this->parameters->FindParam(pvalue,paramenum);
+}
+/*}}}*/
+/*FUNCTION Penta::FindParam(int** pvalues,int* psize,int paramenum){{{*/
+void Penta::FindParam(int** pvalues,int* psize,int paramenum){
+	this->parameters->FindParam(pvalues,psize,paramenum);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16846)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16847)
@@ -78,4 +78,5 @@
 		void   FindParam(int* pvalue,int paramenum);
 		void   FindParam(IssmDouble* pvalue,int paramenum);
+		void   FindParam(int** pvalues,int* psize,int paramenum);
 		int    FiniteElement(void);
 		void   SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp	(revision 16846)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp	(revision 16847)
@@ -98,4 +98,9 @@
 }
 /*}}}*/
+/*FUNCTION Seg::FindParam(int** pvalues,int* psize,int paramenum){{{*/
+void Seg::FindParam(int** pvalues,int* psize,int paramenum){
+	this->parameters->FindParam(pvalues,psize,paramenum);
+}
+/*}}}*/
 /*FUNCTION Seg::FiniteElement{{{*/
 int Seg::FiniteElement(void){
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16846)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16847)
@@ -90,4 +90,5 @@
 		void        FindParam(int* pvalue,int paramenum);
 		void        FindParam(IssmDouble* pvalue,int paramenum);
+		void        FindParam(int** pvalues,int* psize,int paramenum);
 		int         FiniteElement(void);
 		Element*    GetBasalElement(void){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16846)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16847)
@@ -848,4 +848,9 @@
 void Tria::FindParam(IssmDouble* pvalue,int paramenum){
 	this->parameters->FindParam(pvalue,paramenum);
+}
+/*}}}*/
+/*FUNCTION Tria::FindParam(int** pvalues,int* psize,int paramenum){{{*/
+void Tria::FindParam(int** pvalues,int* psize,int paramenum){
+	this->parameters->FindParam(pvalues,psize,paramenum);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16846)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16847)
@@ -83,4 +83,5 @@
 		void        FindParam(int* pvalue,int paramenum);
 		void        FindParam(IssmDouble* pvalue,int paramenum);
+		void        FindParam(int** pvalues,int* psize,int paramenum);
 		int         FiniteElement(void);
 		Element*    GetBasalElement(void){_error_("not implemented yet");};