Index: issm/trunk/src/c/Container/Inputs.cpp
===================================================================
--- issm/trunk/src/c/Container/Inputs.cpp	(revision 5639)
+++ issm/trunk/src/c/Container/Inputs.cpp	(revision 5640)
@@ -71,4 +71,32 @@
 }
 /*}}}*/
+/*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/
+void Inputs::GetParameterValue(double* pvalue,GaussTria* gauss, int enum_type){
+
+	vector<Object*>::iterator object;
+	Input* input=NULL;
+	bool   found=false;
+
+	/*Go through inputs and check whether any input with the same name is already in: */
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		input=(Input*)(*object); 
+		if (input->EnumType()==enum_type){
+			found=true;
+			break;
+		}
+	}
+
+	if (!found){
+		/*we could not find an input with the correct enum type. No defaults values were provided, 
+		 * error out: */
+		ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type));
+	}
+
+	/*Ok, we have an input if we made it here, request the input to return the values: */
+	input->GetParameterValue(pvalue,gauss);
+
+}
+/*}}}*/
 /*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue){{{1*/
 void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type,double defaultvalue){
Index: issm/trunk/src/c/Container/Inputs.h
===================================================================
--- issm/trunk/src/c/Container/Inputs.h	(revision 5639)
+++ issm/trunk/src/c/Container/Inputs.h	(revision 5640)
@@ -16,4 +16,5 @@
 class Input;
 class Node;
+class GaussTria;
 
 class Inputs: public DataSet{
@@ -43,4 +44,5 @@
 		void GetParameterValue(double* pvalue,int enum_type);
 		void GetParameterValue(double* pvalue,double* gauss,int enum_type);
+		void GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type);
 		void GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue);
 		void GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type);
Index: issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5639)
+++ issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5640)
@@ -2158,9 +2158,9 @@
 	switch(analysis_type){
 
-		case DiagnosticHorizAnalysisEnum: case DiagnosticVertAnalysisEnum: 
+		case DiagnosticHorizAnalysisEnum:
 
 			/*default vx,vy and vz: either observation or 0 */
 			if(!iomodel->vx){
-				if (false && iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts; // false TO BE DELETED (for NR CM)
+				if (iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts;
 				else                 for(i=0;i<6;i++)nodeinputs[i]=0;
 				this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs));
@@ -2169,5 +2169,5 @@
 			}
 			if(!iomodel->vy){
-				if (false && iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts; // false TO BE DELETED (for NR CM)
+				if (iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts;
 				else                 for(i=0;i<6;i++)nodeinputs[i]=0;
 				this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs));
@@ -4579,6 +4579,4 @@
 	double L[numdof];
 	double l1l2l3[3];
-	double alpha2_list[3];
-	double basalfriction_list[3]={0.0};
 	double basalfriction;
 	double epsilon[6];
Index: issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5639)
+++ issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5640)
@@ -688,20 +688,24 @@
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
 
-		/*Add Tikhonov regularization term to misfit*/
+		/*Add Tikhonov regularization term to misfit:
+		 *
+		 * WARNING: for now, the regularization is only taken into account by the gradient
+		 * the misfit reflects the response only!
+		 *
+		 * */
 		if (control_type==DragCoefficientEnum){
-			if (!shelf) return 0; //only shelf?... TO BE CHECKED
 			ISSMASSERT(drag_input);
 			drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
-			Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight;
+			//Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight;
 		}
 		else if (control_type==RheologyBbarEnum){
 			ISSMASSERT(Bbar_input);
 			Bbar_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss);
-			Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
+			//Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
 		}
 		else if (control_type==DhDtEnum){
 			ISSMASSERT(dhdt_input);
 			dhdt_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss);
-			Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
+			//Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
 		}
 		else{
@@ -1029,6 +1033,4 @@
 	/*nodal functions: */
 	double l1l2l3[3];
-	double    alpha2complement_list[numvertices];                          //TO BE DELETED
-	double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
 
 	/* strain rate: */
@@ -1068,9 +1070,4 @@
 	friction=new Friction("2d",inputs,matpar,analysis_type);
 
-	/*COMPUT alpha2_list (TO BE DELETED)*/
-	for(i=0;i<numvertices;i++){
-		friction->GetAlphaComplement(&alpha2complement_list[i],&gaussgrids[i][0],VxEnum,VyEnum);
-	}
-
 	/*Retrieve all inputs we will be needing: */
 	adjointx_input=inputs->GetInput(AdjointxEnum);
@@ -1087,7 +1084,6 @@
 
 		/*Build alpha_complement_list: */
-		//if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); // TO BE UNCOMMENTED
-		//else alpha_complement=0;
-		TriaRef::GetParameterValue(&alpha_complement,&alpha2complement_list[0],gauss); // TO BE DELETED
+		if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum);
+		else alpha_complement=0;
 	
 		/*Recover alpha_complement and k: */
@@ -2594,9 +2590,9 @@
 	switch(analysis_type){
 
-		case DiagnosticHorizAnalysisEnum: case DiagnosticVertAnalysisEnum:
+		case DiagnosticHorizAnalysisEnum:
 
 			/*default vx,vy and vz: either observation or 0 */
 			if(!iomodel->vx){
-				if (false && iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; //false TO BE DELETED (for NR CM)
+				if (iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts;
 				else                 for(i=0;i<3;i++)nodeinputs[i]=0;
 				this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs));
@@ -2605,5 +2601,5 @@
 			}
 			if(!iomodel->vy){
-				if (false && iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; //false TO BE DELETED (for NR CM)
+				if (iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts;
 				else                 for(i=0;i<3;i++)nodeinputs[i]=0;
 				this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs));
@@ -3256,6 +3252,4 @@
 	Friction *friction = NULL;
 	double    alpha2;
-	double    alpha2_list[numvertices];                                       //TO BE DELETED
-	double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
 
 	double MAXSLOPE=.06; // 6 %
@@ -3295,9 +3289,4 @@
 	friction=new Friction("2d",inputs,matpar,analysis_type);
 
-	/*COMPUT alpha2_list (TO BE DELETED)*/
-	for(i=0;i<numvertices;i++){
-		friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum);
-	}
-
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussTria(2);
@@ -3307,6 +3296,5 @@
 
 		/*Friction: */
-		// friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
-		TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED
+		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
 
 		// If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
@@ -3388,6 +3376,4 @@
 	Friction *friction = NULL;
 	double    alpha2;
-	double    alpha2_list[numvertices];                                       //TO BE DELETED
-	double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
 
 	double MAXSLOPE=.06; // 6 %
@@ -3426,9 +3412,4 @@
 	friction=new Friction("2d",inputs,matpar,analysis_type);
 
-	/*COMPUT alpha2_list (TO BE DELETED)*/
-	for(i=0;i<numvertices;i++){
-		friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum);
-	}
-
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussTria(2);
@@ -3438,6 +3419,5 @@
 
 		/*Friction: */
-		// friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
-		TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED
+		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
 
 		// If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
@@ -3517,6 +3497,4 @@
 	Friction *friction = NULL;
 	double    alpha2;
-	double    alpha2_list[numvertices];                                       //TO BE DELETED
-	double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
 
 	double MAXSLOPE=.06; // 6 %
@@ -3555,9 +3533,4 @@
 	friction=new Friction("2d",inputs,matpar,analysis_type);
 
-	/*COMPUT alpha2_list (TO BE DELETED)*/
-	for(i=0;i<numvertices;i++){
-		friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum);
-	}
-
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussTria(2);
@@ -3567,6 +3540,5 @@
 
 		/*Friction: */
-		// friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
-		TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED
+		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
 
 		// If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
@@ -5607,9 +5579,4 @@
 	double alpha2,vx,vy;
 	double geothermalflux_value;
-	double alpha2_list[numvertices];                                    //TO BE DELETED
-	double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
-	double vx_list[numvertices]; //TO BE DELETED
-	double vy_list[numvertices]; //TO BE DELETED
-	double basalfriction_list[numvertices]; //TO BE DELETED
 
 	/* gaussian points: */
@@ -5629,7 +5596,7 @@
 
 	/*inputs: */
-	Input* vx_input=NULL; //TO BE DELETED
-	Input* vy_input=NULL; //TO BE DELETED
-	Input* vz_input=NULL; //TO BE DELETED
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	Input* vz_input=NULL;
 	Input* geothermalflux_input=NULL;
 	
@@ -5656,14 +5623,4 @@
 	geothermalflux_input=inputs->GetInput(GeothermalFluxEnum);
 
-	/*COMPUT alpha2_list and basalfriction_list (TO BE DELETED)*/
-	for(i=0;i<numvertices;i++){
-		friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum); //TO BE DELETED
-	}
-	vx_input->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3); //TO BE DELETED
-	vy_input->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3); //TO BE DELETED
-	for(i=0;i<numvertices;i++){
-		basalfriction_list[i]=alpha2_list[i]*(pow(vx_list[i],(double)2.0)+pow(vy_list[i],(double)2.0)); //TO BE DELETED
-	}
-	
 	/* Ice/ocean heat exchange flux on ice shelf base */
 
@@ -5684,6 +5641,5 @@
 	
 		/*Friction: */
-		//friction->GetAlpha2(&alpha2,&gauss_coord[0],VxEnum,VyEnum,VzEnum);
-		TriaRef::GetParameterValue(&basalfriction,&basalfriction_list[0],gauss); // TO BE DELETED
+		friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
 		
 		/*Calculate scalar parameter*/
Index: issm/trunk/src/c/objects/Loads/Friction.cpp
===================================================================
--- issm/trunk/src/c/objects/Loads/Friction.cpp	(revision 5639)
+++ issm/trunk/src/c/objects/Loads/Friction.cpp	(revision 5640)
@@ -116,4 +116,65 @@
 }
 /*}}}*/
+/*FUNCTION Friction::GetAlpha2{{{1*/
+void Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
+
+	/*This routine calculates the basal friction coefficient 
+	  alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
+
+	/*diverse: */
+	double  r,s;
+	double  drag_p, drag_q;
+	double  gravity,rho_ice,rho_water;
+	double  Neff;
+	double  thickness,bed;
+	double  vx,vy,vz,vmag;
+	double  drag_coefficient;
+	double  alpha2;
+
+
+	/*Recover parameters: */
+	inputs->GetParameterValue(&drag_p,DragPEnum);
+	inputs->GetParameterValue(&drag_q,DragQEnum);
+	inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
+	inputs->GetParameterValue(&bed, gauss,BedEnum);
+	inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
+
+	/*Get material parameters: */
+	gravity=matpar->GetG();
+	rho_ice=matpar->GetRhoIce();
+	rho_water=matpar->GetRhoWater();
+
+	//compute r and q coefficients: */
+	r=drag_q/drag_p;
+	s=1./drag_p;
+
+	//From bed and thickness, compute effective pressure when drag is viscous:
+	Neff=gravity*(rho_ice*thickness+rho_water*bed);
+
+	/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
+	  the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
+	  for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
+	  replace it by Neff=0 (ie, equival it to an ice shelf)*/
+	if (Neff<0)Neff=0;
+
+	if(strcmp(element_type,"2d")==0){
+		inputs->GetParameterValue(&vx, gauss,vxenum);
+		inputs->GetParameterValue(&vy, gauss,vyenum);
+		vmag=sqrt(pow(vx,2)+pow(vy,2));
+	}
+	else if (strcmp(element_type,"3d")==0){
+		inputs->GetParameterValue(&vx, gauss,vxenum);
+		inputs->GetParameterValue(&vy, gauss,vyenum);
+		inputs->GetParameterValue(&vz, gauss,vzenum);
+		vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
+	}
+	else ISSMERROR("element_type %s not supported yet",element_type);
+
+	alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
+
+	/*Assign output pointers:*/
+	*palpha2=alpha2;
+}
+/*}}}*/
 /*FUNCTION Friction::GetAlphaComplement {{{1*/
 void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){
@@ -174,2 +235,60 @@
 }
 /*}}}*/
+/*FUNCTION Friction::GetAlphaComplement {{{1*/
+void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){
+
+	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
+	 * FrictionGetAlphaComplement is used in control methods on drag, and it computes: 
+	 * alpha_complement= Neff ^r * vel ^s*/
+
+	/*diverse: */
+	int     i;
+	double  Neff;
+	double  r,s;
+	double  vx;
+	double  vy;
+	double  vmag;
+	double  drag_p,drag_q;
+	double  drag_coefficient;
+	double  bed,thickness;
+	double  gravity,rho_ice,rho_water;
+	double  alpha_complement;
+
+	/*Recover parameters: */
+	inputs->GetParameterValue(&drag_p,DragPEnum);
+	inputs->GetParameterValue(&drag_q,DragQEnum);
+	inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
+	inputs->GetParameterValue(&bed, gauss,BedEnum);
+	inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
+
+	/*Get material parameters: */
+	gravity=matpar->GetG();
+	rho_ice=matpar->GetRhoIce();
+	rho_water=matpar->GetRhoWater();
+
+
+	//compute r and q coefficients: */
+	r=drag_q/drag_p;
+	s=1./drag_p;
+
+	//From bed and thickness, compute effective pressure when drag is viscous:
+	Neff=gravity*(rho_ice*thickness+rho_water*bed);
+
+	/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
+	  the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
+	  for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
+	  replace it by Neff=0 (ie, equival it to an ice shelf)*/
+	if (Neff<0)Neff=0;
+
+	//We need the velocity magnitude to evaluate the basal stress:
+	inputs->GetParameterValue(&vx, gauss,vxenum);
+	inputs->GetParameterValue(&vy, gauss,vyenum);
+	vmag=sqrt(pow(vx,2)+pow(vy,2));
+
+	alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
+
+	/*Assign output pointers:*/
+	*palpha_complement=alpha_complement;
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Loads/Friction.h
===================================================================
--- issm/trunk/src/c/objects/Loads/Friction.h	(revision 5639)
+++ issm/trunk/src/c/objects/Loads/Friction.h	(revision 5640)
@@ -28,5 +28,7 @@
 		void  Echo(void);
 		void  GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum);
+		void  GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
 		void  GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum);
+		void  GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum);
 
 };
