Index: ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp	(revision 18127)
+++ ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp	(revision 18128)
@@ -19,205 +19,226 @@
 #include "./types.h"
 #include "./isnan.h"
 
-void BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs){
+void BrentSearch(IssmDouble** pJ,OptPars optpars,IssmDouble* X0,IssmDouble (*f)(IssmDouble*,void*),IssmDouble (*g)(IssmDouble**,IssmDouble*,void*),void* usr){
 
 	/* This routine is optimizing a given function using Brent's method
 	 * (Golden or parabolic procedure)*/
 
 	/*Intermediary*/
+	int        iter;
 	IssmDouble si,gold,intervalgold,oldintervalgold;
-	IssmDouble parab_num,parab_den;
-	IssmDouble distance,cm_jump;
+	IssmDouble parab_num,parab_den,distance;
 	IssmDouble fxmax,fxmin,fxbest;
 	IssmDouble fx,fx1,fx2;
-	IssmDouble xmax,xmin,xbest;
-	IssmDouble x,x1,x2,xm;
+	IssmDouble x,x1,x2,xm,xbest;
 	IssmDouble tol1,tol2,seps;
 	IssmDouble tolerance = 1.e-4;
-	int        maxiter ,iter;
-	bool       loop= true,goldenflag;
 
 	/*Recover parameters:*/
-	xmin=optpars->xmin;
-	xmax=optpars->xmax;
-	maxiter=optpars->maxiter;
-	cm_jump=optpars->cm_jump;
+	int         nsteps  = optpars.nsteps;
+	int         nsize   = optpars.nsize;
+	IssmDouble  xmin    = optpars.xmin;
+	IssmDouble  xmax    = optpars.xmax;
+	int        *maxiter = optpars.maxiter;
+	IssmDouble *cm_jump = optpars.cm_jump;
 
-	/*initialize counter and get response at the boundaries*/
-	iter=0;
-	fxmin = (*f)(xmin,optargs);
-	if (xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
-	cout<<setprecision(5);
-	if(VerboseControl()) _printf0_("\n");
-	if(VerboseControl()) _printf0_("       Iteration         x           f(x)       Tolerance         Procedure\n");
-	if(VerboseControl()) _printf0_("\n");
-	if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmin<<"  "<<setw(12)<<fxmin<<"           N/A         boundary\n");
-	fxmax = (*f)(xmax,optargs);
-	if (xIsNan<IssmDouble>(fxmax)) _error_("Function evaluation returned NaN");
-	if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmax<<"  "<<setw(12)<<fxmax<<"           N/A         boundary\n");
+	/*Initialize gradient and controls*/
+	IssmDouble* G = NULL;
+	IssmDouble* J = xNew<IssmDouble>(nsteps);
+	IssmDouble* X = xNew<IssmDouble>(nsize);
 
-	/*test if jump option activated and xmin==0*/
-	if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && (fxmax/fxmin)<cm_jump){
-		*psearch_scalar=xmax;
-		*pJ=fxmax;
-		return;
-	}
+	/*start iterations*/
+	for(int n=0;n<nsteps;n++){
 
-	/*initialize optimization variables*/
-	seps=sqrt(DBL_EPSILON);    //precision of a IssmDouble
-	distance=0.0;              //new_x=old_x + distance
-	gold=0.5*(3.0-sqrt(5.0));  //gold = 1 - golden ratio
-	intervalgold=0.0;          //distance used by Golden procedure
+		/*Reset some variables*/
+		iter = 0;
+		xmin = 0.;
+		xmax = 1.;
+		bool loop = true;
+		cout<<setprecision(5);
 
-	/*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
-	x1=xmin+gold*(xmax-xmin);
-	x2=x1;
-	xbest=x1;
-	x=xbest;
+		/*Get current Gradient at xmin=0*/
+		if(VerboseControl()) _printf0_("\n" << "   step " << n+1 << "/" << nsteps << "\n");
+		fxmin = (*g)(&G,X0,usr); if(xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
+		if(VerboseControl()) _printf0_("\n");
+		if(VerboseControl()) _printf0_("       Iteration         x           f(x)       Tolerance         Procedure\n");
+		if(VerboseControl()) _printf0_("\n");
+		if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmin<<"  "<<setw(12)<<fxmin<<"           N/A         boundary\n");
+		
+		/*Get f(xmax)*/
+		for(int i=0;i<nsize;i++) X[i]=X0[i]+xmax*G[i];
+		fxmax = (*f)(X,usr); if (xIsNan<IssmDouble>(fxmax)) _error_("Function evaluation returned NaN");
+		if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmax<<"  "<<setw(12)<<fxmax<<"           N/A         boundary\n");
 
-	/*2: call the function to be evaluated*/
-	fxbest = (*f)(x,optargs);
-	if(xIsNan<IssmDouble>(fxbest)) _error_("Function evaluation returned NaN");
-	iter=iter+1;
+		/*test if jump option activated and xmin==0*/
+		if(!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && (fxmax/fxmin)<cm_jump[n]){
+			for(int i=0;i<nsize;i++) X0[i]=X0[i]+xmax*G[i];
+			xDelete<IssmDouble>(G);
+			J[n]=fxmax;
+			continue;
+		}
 
-	/*3: update the other variables*/
-	fx1=fxbest;
-	fx2=fxbest;
-	/*xm is always in the middle of a and b*/
-	xm=0.5*(xmin+xmax);                           
-	/*update tolerances*/
-	tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
-	tol2=2.0*tol1;
+		/*initialize optimization variables*/
+		seps=sqrt(DBL_EPSILON);    //precision of a IssmDouble
+		distance=0.0;              //new_x=old_x + distance
+		gold=0.5*(3.0-sqrt(5.0));  //gold = 1 - golden ratio
+		intervalgold=0.0;          //distance used by Golden procedure
 
-	/*4: print result*/
-	if(VerboseControl())
-	 _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<xbest<<"  "<<setw(12)<<fxbest<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<"         initial\n");
-	if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
-		if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump << "\n");
-		loop=false;
-	}
+		/*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
+		x1=xmin+gold*(xmax-xmin);
+		x2=x1;
+		xbest=x1;
+		x=xbest;
 
-	while(loop){
+		/*2: call the function to be evaluated*/
+		for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
+		fxbest = (*f)(X,usr); if(xIsNan<IssmDouble>(fxbest)) _error_("Function evaluation returned NaN");
+		iter++;
 
-		goldenflag=true;
+		/*3: update the other variables*/
+		fx1=fxbest;
+		fx2=fxbest;
+		/*xm is always in the middle of a and b*/
+		xm=0.5*(xmin+xmax);                           
+		/*update tolerances*/
+		tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
+		tol2=2.0*tol1;
 
-		// Is a parabolic fit possible ?
-		if (sqrt(pow(intervalgold,2))>tol1){
+		/*4: print result*/
+		if(VerboseControl())
+		 _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<xbest<<"  "<<setw(12)<<fxbest<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<"         initial\n");
+		if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
+			if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
+			loop=false;
+		}
 
-			// Yes, so fit parabola
-			goldenflag=false;
-			parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
-			parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
+		while(loop){
 
-			//reverse p if necessary
-			if(parab_den>0.0){ 
-				parab_num=-parab_num;
-			}
-			parab_den=sqrt(pow(parab_den,2));
-			oldintervalgold=intervalgold;
-			intervalgold=distance;
+			bool goldenflag=true;
 
-			// Is the parabola acceptable (we use seps here because in some configuration parab_num==parab_den*(xmax-xbest)
-			// and the result is not repeatable anymore
-			if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
-						(parab_num>parab_den*(xmin-xbest)+seps) && 
-						(parab_num<parab_den*(xmax-xbest)-seps)){
+			// Is a parabolic fit possible ?
+			if (sqrt(pow(intervalgold,2))>tol1){
 
-				// Yes, parabolic interpolation step
-				distance=parab_num/parab_den;
-				x=xbest+distance;
+				// Yes, so fit parabola
+				goldenflag=false;
+				parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
+				parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
 
-				// f must not be evaluated too close to min_x or max_x
-				if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
-					if ((xm-xbest)<0.0) si=-1;
-					else                si=1;
-					//compute new distance
-					distance=tol1*si;
+				//reverse p if necessary
+				if(parab_den>0.0){ 
+					parab_num=-parab_num;
 				}
+				parab_den=sqrt(pow(parab_den,2));
+				oldintervalgold=intervalgold;
+				intervalgold=distance;
+
+				// Is the parabola acceptable (we use seps here because in some configuration parab_num==parab_den*(xmax-xbest)
+				// and the result is not repeatable anymore
+				if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
+							(parab_num>parab_den*(xmin-xbest)+seps) && 
+							(parab_num<parab_den*(xmax-xbest)-seps)){
+
+					// Yes, parabolic interpolation step
+					distance=parab_num/parab_den;
+					x=xbest+distance;
+
+					// f must not be evaluated too close to min_x or max_x
+					if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
+						if ((xm-xbest)<0.0) si=-1;
+						else                si=1;
+						//compute new distance
+						distance=tol1*si;
+					}
+				}
+				else{
+					// Not acceptable, must do a golden section step
+					goldenflag=true;
+				}
 			}
-			else{
-				// Not acceptable, must do a golden section step
-				goldenflag=true;
-			}
-		}
 
-		//Golden procedure
-		if(goldenflag){
-			// compute the new distance d
-			if(xbest>=xm){
-				intervalgold=xmin-xbest;    
+			//Golden procedure
+			if(goldenflag){
+				// compute the new distance d
+				if(xbest>=xm){
+					intervalgold=xmin-xbest;    
+				}
+				else{ 
+					intervalgold=xmax-xbest;  
+				}
+				distance=gold*intervalgold;
 			}
-			else{ 
-				intervalgold=xmax-xbest;  
-			}
-			distance=gold*intervalgold;
-		}
 
-		// The function must not be evaluated too close to xbest
-		if(distance<0) si=-1;
-		else           si=1;
-		if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
-		else                           x=xbest+si*tol1;
+			// The function must not be evaluated too close to xbest
+			if(distance<0) si=-1;
+			else           si=1;
+			if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
+			else                           x=xbest+si*tol1;
 
-		//evaluate function on x
-		fx = (*f)(x,optargs);
-		if(xIsNan<IssmDouble>(fx)) _error_("Function evaluation returned NaN");
-		iter=iter+1;
+			//evaluate function on x
+			for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
+			fx = (*f)(X,usr); if(xIsNan<IssmDouble>(fx)) _error_("Function evaluation returned NaN");
+			iter++;
 
-		// Update a, b, xm, x1, x2, tol1, tol2
-		if (fx<=fxbest){
-			if (x>=xbest) xmin=xbest;
-			else          xmax=xbest;
-			x1=x2;    fx1=fx2;
-			x2=xbest; fx2=fxbest;
-			xbest=x;  fxbest=fx;
-		}
-		else{ // fx > fxbest
-			if (x<xbest) xmin=x;
-			else         xmax=x;
-			if ((fx<=fx2) || (x2==xbest)){
-				x1=x2; fx1=fx2;
-				x2=x;  fx2=fx;
+			// Update a, b, xm, x1, x2, tol1, tol2
+			if (fx<=fxbest){
+				if (x>=xbest) xmin=xbest;
+				else          xmax=xbest;
+				x1=x2;    fx1=fx2;
+				x2=xbest; fx2=fxbest;
+				xbest=x;  fxbest=fx;
 			}
-			else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
-				x1=x;  fx1=fx;
+			else{ // fx > fxbest
+				if (x<xbest) xmin=x;
+				else         xmax=x;
+				if ((fx<=fx2) || (x2==xbest)){
+					x1=x2; fx1=fx2;
+					x2=x;  fx2=fx;
+				}
+				else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
+					x1=x;  fx1=fx;
+				}
 			}
-		}
-		xm = 0.5*(xmin+xmax);
-		tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
-		tol2=2.0*tol1;
-		if(VerboseControl())
-		 _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<x<<"  "<<setw(12)<<fx<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<
-					 "         "<<(goldenflag?"golden":"parabolic")<<"\n");
+			xm = 0.5*(xmin+xmax);
+			tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
+			tol2=2.0*tol1;
+			if(VerboseControl())
+			 _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<x<<"  "<<setw(12)<<fx<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<
+						 "         "<<(goldenflag?"golden":"parabolic")<<"\n");
 
-		/*Stop the optimization?*/
-		if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
-			if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'tolx'=" << tolerance << "\n");
-			loop=false;
+			/*Stop the optimization?*/
+			if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
+				if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'tolx'=" << tolerance << "\n");
+				loop=false;
+			}
+			else if (iter>=maxiter[n]){
+				if(VerboseControl()) _printf0_("      exiting: Maximum number of iterations has been exceeded  ('maxiter'=" << maxiter[n] << ")\n");
+				loop=false;
+			}
+			else if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
+				if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
+				loop=false;
+			}
+			else{
+				//continue
+				loop=true;
+			}
+		}//end while
+
+		//Now, check that the value on the boundaries are not better than current fxbest
+		if (fxbest>fxmin){
+			xbest=optpars.xmin; fxbest=fxmin;
 		}
-		else if (iter>=maxiter){
-			if(VerboseControl()) _printf0_("      exiting: Maximum number of iterations has been exceeded  ('maxiter'=" << maxiter << ")\n");
-			loop=false;
+		if (fxbest>fxmax){
+			xbest=optpars.xmax; fxbest=fxmax;
 		}
-		else if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
-			if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump << "\n");
-			loop=false;
-		}
-		else{
-			//continue
-			loop=true;
-		}
-	}//end while
 
-	//Now, check that the value on the boundaries are not better than current fxbest
-	if (fxbest>fxmin){
-		xbest=optpars->xmin; fxbest=fxmin;
+		/*Assign output pointers: */
+		for(int i=0;i<nsize;i++) X0[i]=X0[i]+xbest*G[i];
+		xDelete<IssmDouble>(G);
+		J[n]=fxbest;
 	}
-	if (fxbest>fxmax){
-		xbest=optpars->xmax; fxbest=fxmax;
-	}
-
-	/*Assign output pointers: */
-	*psearch_scalar=xbest;
-	*pJ=fxbest;
+	
+	/*return*/
+	xDelete<IssmDouble>(X);
+	*pJ=J;
 }
Index: ../trunk-jpl/src/c/shared/Numerics/numerics.h
===================================================================
--- ../trunk-jpl/src/c/shared/Numerics/numerics.h	(revision 18127)
+++ ../trunk-jpl/src/c/shared/Numerics/numerics.h	(revision 18128)
@@ -29,7 +29,7 @@
 
 int         min(int a,int b);
 int         max(int a,int b);
-void        BrentSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs);
+void        BrentSearch(IssmDouble** pJ,OptPars optpars,IssmDouble* X0,IssmDouble (*f)(IssmDouble*,void*),IssmDouble (*g)(IssmDouble**,IssmDouble*,void*),void* usr);
 void        cross(IssmDouble *result,IssmDouble*vector1,IssmDouble*vector2);
 void        XZvectorsToCoordinateSystem(IssmDouble *T,IssmDouble*xzvectors);
 int         cubic(IssmDouble a, IssmDouble b, IssmDouble c, IssmDouble d, double X[3], int *num);
Index: ../trunk-jpl/src/c/shared/Numerics/OptPars.h
===================================================================
--- ../trunk-jpl/src/c/shared/Numerics/OptPars.h	(revision 18127)
+++ ../trunk-jpl/src/c/shared/Numerics/OptPars.h	(revision 18128)
@@ -9,10 +9,12 @@
 
 struct OptPars{
 
-	IssmDouble xmin;
-	IssmDouble xmax;
-	IssmDouble cm_jump;
-	int maxiter;
+	IssmDouble  xmin;
+	IssmDouble  xmax;
+	IssmDouble *cm_jump;
+	int* maxiter;
+	int  nsteps;
+	int  nsize;
 
 };
 
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 18127)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 18128)
@@ -15,10 +15,10 @@
 	int         num_control_type;
 	int         num_cm_responses;
 	int        *control_type     = NULL;
+	int        *maxiter          = NULL;
 	int        *cm_responses     = NULL;
 	IssmDouble *cm_jump          = NULL;
 	IssmDouble *optscal          = NULL;
-	IssmDouble *maxiter          = NULL;
 
 	/*retrieve some parameters: */
 	iomodel->Constant(&control_analysis,InversionIscontrolEnum);
@@ -55,7 +55,7 @@
 				iomodel->FetchData(&maxiter,NULL,NULL,InversionMaxiterPerStepEnum);
 				parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_control_type));
 				parameters->AddObject(new DoubleVecParam(InversionStepThresholdEnum,cm_jump,nsteps));
-				parameters->AddObject(new DoubleVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
+				parameters->AddObject(new IntVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
 				break;
 			case 1:/*TAO*/
 				parameters->AddObject(iomodel->CopyConstantObject(InversionFatolEnum));
@@ -82,8 +82,8 @@
 
 		xDelete<int>(control_type);
 		xDelete<int>(cm_responses);
+		xDelete<int>(maxiter);
 		iomodel->DeleteData(cm_jump,InversionStepThresholdEnum);
 		iomodel->DeleteData(optscal,InversionGradientScalingEnum);
-		iomodel->DeleteData(maxiter,InversionMaxiterPerStepEnum);
 	}
 }
Index: ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp	(revision 18127)
+++ ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp	(revision 18128)
@@ -6,18 +6,18 @@
 #include "../../shared/shared.h"
 #include "../../toolkits/toolkits.h"
 
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,IssmDouble* vector){
+void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector){
 
 	int  num_controls;
 	int *control_type = NULL;
 
 	/*Retrieve some parameters*/
-	parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
-	parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
+	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+	femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
 
 	for(int i=0;i<num_controls;i++){
-		for(int j=0;j<elements->Size();j++){
-			Element* element=(Element*)elements->GetObjectByOffset(j);
+		for(int j=0;j<femmodel->elements->Size();j++){
+			Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
 			element->SetControlInputsFromVector(vector,control_type[i],i);
 		}
 	}
@@ -25,14 +25,9 @@
 	xDelete<int>(control_type);
 }
 
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Vector<IssmDouble>* vector){
+void SetControlInputsFromVectorx(FemModel* femmodel,Vector<IssmDouble>* vector){
 
-	IssmDouble* serial_vector=NULL;
-
-	serial_vector=vector->ToMPISerial();
-
-	SetControlInputsFromVectorx(elements,nodes, vertices, loads, materials, parameters,serial_vector);
-
-	/*Free ressources:*/
+	IssmDouble* serial_vector=vector->ToMPISerial();
+	SetControlInputsFromVectorx(femmodel,serial_vector);
 	xDelete<IssmDouble>(serial_vector);
 }
Index: ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h
===================================================================
--- ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h	(revision 18127)
+++ ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h	(revision 18128)
@@ -7,7 +7,7 @@
 #include "../../classes/classes.h"
 
 /* local prototypes: */
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,Vector<IssmDouble>* vector);
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,IssmDouble* vector);
+void SetControlInputsFromVectorx(FemModel* femmodel,Vector<IssmDouble>* vector);
+void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector);
 
 #endif 
Index: ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp	(revision 18127)
+++ ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp	(revision 18128)
@@ -49,6 +49,8 @@
 		gradient->AXPY(gradient_list[i],1.0);
 		delete gradient_list[i];
 	}
+	//gradient->Echo();
+	//_error_("S");
 
 	/*Check that gradient is clean*/
 	norm_inf=gradient->Norm(NORM_INF);
Index: ../trunk-jpl/src/c/cores/controlvalidation_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/controlvalidation_core.cpp	(revision 18127)
+++ ../trunk-jpl/src/c/cores/controlvalidation_core.cpp	(revision 18128)
@@ -72,7 +72,7 @@
 		for(int i=0;i<n;i++) X[i] = X0[i] + alpha;
 
 		/*Calculate j(k+alpha delta k) */
-		SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+		SetControlInputsFromVectorx(femmodel,X);
 		solutioncore(femmodel);
 		femmodel->CostFunctionx(&j,NULL,NULL);
 
Index: ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 18127)
+++ ../trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 18128)
@@ -105,7 +105,7 @@
 	}
 
 	/*Get solution*/
-	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+	SetControlInputsFromVectorx(femmodel,X);
 	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
 	femmodel->OutputControlsx(&femmodel->results);
 	femmodel->results->AddObject(new GenericExternalResult<double>(JEnum,f,1,0));
@@ -130,10 +130,9 @@
 	int         solution_type;
 	FemModel   *femmodel  = (FemModel*)dzs;
 
-	/*Recover responses*/
-	int         num_responses;
-	int        *responses = NULL;
-	femmodel->parameters->FindParam(&responses,&num_responses,InversionCostFunctionsEnum);
+	/*Recover number of cost functions responses*/
+	int num_responses;
+	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
 
 	/*Constrain input vector*/
 	IssmDouble  *XL = NULL;
@@ -146,7 +145,7 @@
 	}
 
 	/*Update control input*/
-	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+	SetControlInputsFromVectorx(femmodel,X);
 
 	/*Recover some parameters*/
 	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
@@ -169,7 +168,6 @@
 
 	if(indic==0){
 		/*dry run, no gradient required*/
-		xDelete<int>(responses);
 		xDelete<IssmDouble>(XU);
 		xDelete<IssmDouble>(XL);
 		return;
@@ -192,7 +190,6 @@
 	}
 
 	/*Clean-up and return*/
-	xDelete<int>(responses);
 	xDelete<IssmDouble>(XU);
 	xDelete<IssmDouble>(XL);
 }
Index: ../trunk-jpl/src/c/cores/control_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/control_core.cpp	(revision 18127)
+++ ../trunk-jpl/src/c/cores/control_core.cpp	(revision 18128)
@@ -10,15 +10,19 @@
 #include "../solutionsequences/solutionsequences.h"
 
 /*Local prototypes*/
-bool controlconvergence(IssmDouble J, IssmDouble tol_cm);
-IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs);
+IssmDouble FormFunction(IssmDouble* X,void* usr);
+IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usr);
+typedef struct {
+	FemModel* femmodel;
+	int       nsize;
+} AppCtx;
 
 void control_core(FemModel* femmodel){
 
 	int     i;
 
 	/*parameters: */
-	int        num_controls;
+	int        num_controls,nsize;
 	int        nsteps;
 	IssmDouble tol_cm;
 	int        solution_type;
@@ -26,11 +30,10 @@
 	bool       dakota_analysis = false;
 
 	int        *control_type   = NULL;
-	IssmDouble *maxiter        = NULL;
+	int*        maxiter        = NULL;
 	IssmDouble *cm_jump        = NULL;
 
 	/*intermediary: */
-	IssmDouble search_scalar = 1.;
 	OptPars    optpars;
 
 	/*Solution and Adjoint core pointer*/
@@ -60,37 +63,41 @@
 	if(VerboseControl()) _printf0_("   preparing initial solution\n");
 	if(isFS) solutioncore(femmodel);
 
-	/*Initialize cost function: */
-	J=xNew<IssmDouble>(nsteps);
+	/*Get initial guess*/
+	Vector<IssmDouble> *Xpetsc = NULL;
+	GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
+	IssmDouble* X0 = Xpetsc->ToMPISerial();
+	Xpetsc->GetSize(&nsize);
+	delete Xpetsc;
 
 	/*Initialize some of the BrentSearch arguments: */
-	optpars.xmin=0; optpars.xmax=1;
+	optpars.xmin    = 0; 
+	optpars.xmax    = 1;
+	optpars.nsteps  = nsteps;
+	optpars.nsize   = nsize;
+	optpars.maxiter = maxiter;
+	optpars.cm_jump = cm_jump;
 
-	/*Start looping: */
-	for(int n=0;n<nsteps;n++){
+	/*Initialize function argument*/
+	AppCtx usr;
+	usr.femmodel = femmodel;
+	usr.nsize    = nsize;
 
-		/*Display info*/
-		if(VerboseControl()) _printf0_("\n" << "   control method step " << n+1 << "/" << nsteps << "\n");
+	/*Call Brent optimization*/
+	BrentSearch(&J,optpars,X0,&FormFunction,&FormFunctionGradient,(void*)&usr);
 
-
-		/*In steady state inversion, compute new temperature field now*/
-		if(solution_type==SteadystateSolutionEnum) solutioncore(femmodel);
-
-		if(VerboseControl()) _printf0_("   compute adjoint state:\n");
-		adjointcore(femmodel);
-		gradient_core(femmodel,n,search_scalar==0.);
-
-		if(VerboseControl()) _printf0_("   optimizing along gradient direction\n");
-		optpars.maxiter=reCast<int,IssmDouble>(maxiter[n]); optpars.cm_jump=cm_jump[n];
-		BrentSearch(&search_scalar,J+n,&optpars,&objectivefunction,(void*)femmodel);
-
-		if(VerboseControl()) _printf0_("   updating parameter using optimized search scalar\n"); //true means update save controls
-		InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true);
-
-		if(controlconvergence(J[n],tol_cm)) break;
+	if(VerboseControl()) _printf0_("   preparing final solution\n");
+	IssmDouble  *XL = NULL;
+	IssmDouble  *XU = NULL;
+	GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
+	GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+	for(long i=0;i<nsize;i++){
+		if(X0[i]>XU[i]) X0[i]=XU[i];
+		if(X0[i]<XL[i]) X0[i]=XL[i];
 	}
-
-	if(VerboseControl()) _printf0_("   preparing final solution\n");
+	xDelete<IssmDouble>(XU);
+	xDelete<IssmDouble>(XL);
+	SetControlInputsFromVectorx(femmodel,X0);
 	femmodel->parameters->SetParam(true,SaveResultsEnum);
 	solutioncore(femmodel);
 
@@ -110,84 +117,187 @@
 
 	/*Free ressources: */
 	xDelete<int>(control_type);
-	xDelete<IssmDouble>(maxiter);
+	xDelete<int>(maxiter);
 	xDelete<IssmDouble>(cm_jump);
 	xDelete<IssmDouble>(J);
+	xDelete<IssmDouble>(X0);
 }
-bool controlconvergence(IssmDouble J, IssmDouble tol_cm){
 
-	bool converged=false;
+IssmDouble FormFunction(IssmDouble* X,void* usrvoid){
 
-	/*Has convergence been reached?*/
-	if (!xIsNan<IssmDouble>(tol_cm) && J<tol_cm){
-		converged=true;
-		if(VerboseConvergence()) _printf0_("      Convergence criterion reached: J = " << J << " < " << tol_cm);
-	}
-
-	return converged;
-}
-
-IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs){
-
 	/*output: */
 	IssmDouble J;
 
 	/*parameters: */
-	int        solution_type,analysis_type;
-	bool       isFS       = false;
+	int        solution_type,analysis_type,num_responses;
+	bool       isFS           = false;
 	bool       conserve_loads = true;
-	FemModel  *femmodel       = (FemModel*)optargs;
+	AppCtx*    usr = (AppCtx*)usrvoid;
+	FemModel  *femmodel  = usr->femmodel;
+	int        nsize     = usr->nsize;
 
 	/*Recover parameters: */
 	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
 	femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
 
-	/*set analysis type to compute velocity: */
+	/*Constrain input vector*/
+	IssmDouble  *XL = NULL;
+	IssmDouble  *XU = NULL;
+	GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
+	GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+	for(long i=0;i<nsize;i++){
+		if(X[i]>XU[i]) X[i]=XU[i];
+		if(X[i]<XL[i]) X[i]=XL[i];
+	}
+
+	/*Update control input*/
+	SetControlInputsFromVectorx(femmodel,X);
+
+	/*solve forward: */
 	switch(solution_type){
 		case SteadystateSolutionEnum:
+			femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+			stressbalance_core(femmodel);	//We need a 3D velocity!! (vz is required for the next thermal run)
+			break;
 		case StressbalanceSolutionEnum:
 			femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+			solutionsequence_nonlinear(femmodel,conserve_loads); 
 			break;
 		case BalancethicknessSolutionEnum:
 			femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+			solutionsequence_linear(femmodel); 
 			break;
 		case BalancethicknessSoftSolutionEnum:
-			femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+			/*NOTHING*/
 			break;
 		case Balancethickness2SolutionEnum:
 			femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
+			solutionsequence_linear(femmodel); 
 			break;
 		default:
 			_error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
 	}
 
-	/*update parameter according to scalar: */ //false means: do not save control
-	InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,false);
+	/*Compute misfit for this velocity field.*/
+	IssmDouble* Jlist = NULL;
+	femmodel->CostFunctionx(&J,&Jlist,NULL);
+	//_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
 
-	/*Run stressbalance with updated inputs: */
-	if (solution_type==SteadystateSolutionEnum){
-		stressbalance_core(femmodel);	//We need a 3D velocity!! (vz is required for the next thermal run)
+	/*Retrieve objective functions independently*/
+	//for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
+	//_printf0_("\n");
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(XU);
+	xDelete<IssmDouble>(XL);
+	xDelete<IssmDouble>(Jlist);
+	return J;
+}
+IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usrvoid){
+
+	/*output: */
+	IssmDouble J;
+
+	/*parameters: */
+	void (*adjointcore)(FemModel*)=NULL;
+	int         solution_type,analysis_type,num_responses,num_controls,numvertices;
+	bool        isFS           = false;
+	bool        conserve_loads = true;
+	IssmDouble *scalar_list    = NULL;
+	IssmDouble *Jlist          = NULL;
+	IssmDouble *G              = NULL;
+	IssmDouble *norm_list      = NULL;
+	AppCtx     *usr            = (AppCtx*)usrvoid;
+	FemModel   *femmodel       = usr->femmodel;
+	int         nsize          = usr->nsize;
+
+	/*Recover parameters: */
+	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+	femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+	femmodel->parameters->FindParam(&scalar_list,NULL,NULL,InversionGradientScalingEnum);
+	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);   _assert_(num_controls);
+	numvertices=femmodel->vertices->NumberOfVertices();
+
+	/*Constrain input vector*/
+	IssmDouble  *XL = NULL;
+	IssmDouble  *XU = NULL;
+	GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
+	GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+	for(long i=0;i<nsize;i++){
+		if(X[i]>XU[i]) X[i]=XU[i];
+		if(X[i]<XL[i]) X[i]=XL[i];
 	}
-	else if (solution_type==StressbalanceSolutionEnum){
-		solutionsequence_nonlinear(femmodel,conserve_loads); 
+
+	/*Update control input*/
+	SetControlInputsFromVectorx(femmodel,X);
+
+	/*Compute Adjoint*/
+	AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
+	adjointcore(femmodel);
+
+	/*Compute gradient*/
+	Gradjx(&G,&norm_list,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+
+	/*Compute scaling factor*/
+	IssmDouble scalar = scalar_list[0]/norm_list[0];
+	for(int i=1;i<num_controls;i++) scalar=min(scalar,scalar_list[i]/norm_list[i]);
+
+	/*Constrain Gradient*/
+	for(int i=0;i<num_controls;i++){
+		for(int j=0;j<numvertices;j++){
+			G[i*numvertices+j] = scalar*G[i*numvertices+j];
+		}
 	}
-	else if (solution_type==BalancethicknessSolutionEnum){
-		solutionsequence_linear(femmodel); 
+
+	/*Needed for output results (FIXME: should be placed 6 lines below)*/
+	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
+
+	for(long i=0;i<nsize;i++){
+		if(X[i]>=XU[i]) G[i]=0.;
+		if(X[i]<=XL[i]) G[i]=0.;
 	}
-	else if (solution_type==Balancethickness2SolutionEnum){
-		solutionsequence_linear(femmodel); 
+
+	/*solve forward: (FIXME: not needed actually...)*/
+	switch(solution_type){
+		case SteadystateSolutionEnum:
+			femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+			stressbalance_core(femmodel);	//We need a 3D velocity!! (vz is required for the next thermal run)
+			break;
+		case StressbalanceSolutionEnum:
+			femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+			solutionsequence_nonlinear(femmodel,conserve_loads); 
+			break;
+		case BalancethicknessSolutionEnum:
+			femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+			solutionsequence_linear(femmodel); 
+			break;
+		case BalancethicknessSoftSolutionEnum:
+			/*NOTHING*/
+			break;
+		case Balancethickness2SolutionEnum:
+			femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
+			solutionsequence_linear(femmodel); 
+			break;
+		default:
+			_error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
 	}
-	else if (solution_type==BalancethicknessSoftSolutionEnum){
-		/*Don't do anything*/
-	}
-	else{
-		_error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
-	}
 
 	/*Compute misfit for this velocity field.*/
-	femmodel->CostFunctionx(&J,NULL,NULL);
+	femmodel->CostFunctionx(&J,&Jlist,NULL);
+	//_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
+	//for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
+	//_printf0_("\n");
 
-	/*Free ressources:*/
+	/*Clean-up and return*/
+	xDelete<IssmDouble>(XU);
+	xDelete<IssmDouble>(XL);
+	xDelete<IssmDouble>(norm_list);
+	xDelete<IssmDouble>(scalar_list);
+	xDelete<IssmDouble>(Jlist);
+	*pG = G;
 	return J;
 }
Index: ../trunk-jpl/src/c/cores/controltao_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/controltao_core.cpp	(revision 18127)
+++ ../trunk-jpl/src/c/cores/controltao_core.cpp	(revision 18128)
@@ -92,7 +92,7 @@
 	TaoGetSolutionVector(tao,&X->pvector->vector);
 	G=new Vector<IssmDouble>(0); VecFree(&G->pvector->vector);
 	TaoGetGradientVector(tao,&G->pvector->vector);
-	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+	SetControlInputsFromVectorx(femmodel,X);
 	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
 	femmodel->OutputControlsx(&femmodel->results);
 	femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,user.J,maxiter+3,1,1,0));
@@ -112,11 +112,11 @@
 	TaoDestroy(&tao);
 	TaoFinalize();
 }
-int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *userCtx){
+int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *uservoid){
 
 	/*Retreive arguments*/
 	int                  solution_type;
-	AppCtx              *user            = (AppCtx *)userCtx;
+	AppCtx              *user            = (AppCtx *)uservoid;
 	FemModel            *femmodel        = user->femmodel;
 	Vector<IssmDouble>  *gradient        = NULL;
 	Vector<IssmDouble>  *X               = NULL;
@@ -126,7 +126,7 @@
 
 	/*Set new variable*/
 	//VecView(X,PETSC_VIEWER_STDOUT_WORLD);
-	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+	SetControlInputsFromVectorx(femmodel,X);
 	delete X;
 
 	/*Recover some parameters*/
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18127)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18128)
@@ -983,7 +983,6 @@
 		this->inputs->AddInput(new IntInput(ApproximationEnum,reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[index])));
 	}
 
-
 	/*Control Inputs*/
 	if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
 		for(i=0;i<num_control_type;i++){
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18127)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18128)
@@ -2850,6 +2850,7 @@
 	GradientIndexing(&vertexpidlist[0],control_index);
 
 	/*Get input (either in element or material)*/
+	if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
 	Input* input=inputs->GetInput(control_enum);
 	if(!input) _error_("Input " << EnumToStringx(control_enum) << " not found in element");
 
@@ -2864,10 +2865,21 @@
 void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
 
 	IssmDouble  values[NUMVERTICES];
-	int     vertexpidlist[NUMVERTICES];
-	Input  *input     = NULL;
-	Input  *new_input = NULL;
+	int         vertexpidlist[NUMVERTICES],control_init;
+	Input      *input     = NULL;
+	Input      *new_input = NULL;
 
+	/*Specific case for depth averaged quantities*/
+	control_init=control_enum;
+	if(control_enum==MaterialsRheologyBbarEnum){
+		control_enum=MaterialsRheologyBEnum;
+		if(!IsOnBase()) return;
+	}
+	if(control_enum==DamageDbarEnum){
+		control_enum=DamageDEnum;
+		if(!IsOnBase()) return;
+	}
+
 	/*Get out if this is not an element input*/
 	if(!IsInput(control_enum)) return;
 
@@ -2875,7 +2887,7 @@
 	GradientIndexing(&vertexpidlist[0],control_index);
 
 	/*Get values on vertices*/
-	for (int i=0;i<NUMVERTICES;i++){
+	for(int i=0;i<NUMVERTICES;i++){
 		values[i]=vector[vertexpidlist[i]];
 	}
 	new_input = new PentaInput(control_enum,values,P1Enum);
@@ -2886,6 +2898,13 @@
 	}
 
 	((ControlInput*)input)->SetInput(new_input);
+
+	if(control_init==MaterialsRheologyBbarEnum){
+		this->InputExtrude(control_enum);
+	}
+	if(control_init==DamageDbarEnum){
+		this->InputExtrude(control_enum);
+	}
 }
 /*}}}*/
 
Index: ../trunk-jpl/src/m/classes/inversion.m
===================================================================
--- ../trunk-jpl/src/m/classes/inversion.m	(revision 18127)
+++ ../trunk-jpl/src/m/classes/inversion.m	(revision 18128)
@@ -24,55 +24,55 @@
 		thickness_obs               = NaN
 	end
 	methods
-         function createxml(obj,fid) % {{{
-            fprintf(fid, '<!-- inversion -->\n');            
-                    
-            % inversion parameters
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="inversion parameters">','<section name="inversion" />');                    
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="iscontrol" type="',class(obj.iscontrol),'" default="',convert2str(obj.iscontrol),'">','     <section name="inversion" />','     <help> is inversion activated? </help>','  </parameter>');
-            
-            % incompleteadjoing drop-down (0 or 1)
-            fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="incomplete_adjoint" type="alternative" optional="false">','     <section name="inversion" />','     <help> 1: linear viscosity, 0: non-linear viscosity </help>');
-            fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n%s\n','       <option value="1" type="string" default="false"> </option>','</parameter>');
-            
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="control_parameters" type="',class(obj.control_parameters),'" default="',convert2str(obj.control_parameters),'">','     <section name="inversion" />','     <help> ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''} </help>','  </parameter>');
-                
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="nsteps" type="',class(obj.nsteps),'" default="',convert2str(obj.nsteps),'">','     <section name="inversion" />','     <help> number of optimization searches </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions" type="',class(obj.cost_functions),'" default="',convert2str(obj.cost_functions),'">','     <section name="inversion" />','     <help> indicate the type of response for each optimization step  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions_coefficients" type="',class(obj.cost_functions_coefficients),'" default="',convert2str(obj.cost_functions_coefficients),'">','     <section name="inversion" />','     <help> cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter </help>','  </parameter>');
-                
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_function_threshold" type="',class(obj.cost_function_threshold),'" default="',convert2str(obj.cost_function_threshold),'">','     <section name="inversion" />','     <help> misfit convergence criterion. Default is 1%, NaN if not applied </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter_per_step" type="',class(obj.maxiter_per_step),'" default="',convert2str(obj.maxiter_per_step),'">','     <section name="inversion" />','     <help> maximum iterations during each optimization step  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="gradient_scaling" type="',class(obj.gradient_scaling),'" default="',convert2str(obj.gradient_scaling),'">','     <section name="inversion" />','     <help> scaling factor on gradient direction during optimization, for each optimization step </help>','  </parameter>');
-               
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="step_threshold" type="',class(obj.step_threshold),'" default="',convert2str(obj.step_threshold),'">','     <section name="inversion" />','     <help> decrease threshold for misfit, default is 30% </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="min_parameters" type="',class(obj.min_parameters),'" default="',convert2str(obj.min_parameters),'">','     <section name="inversion" />','     <help> absolute minimum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="max_parameters" type="',class(obj.max_parameters),'" default="',convert2str(obj.max_parameters),'">','     <section name="inversion" />','     <help> absolute maximum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
-               
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vx_obs" type="',class(obj.vx_obs),'" default="',convert2str(obj.vx_obs),'">','     <section name="inversion" />','     <help> observed velocity x component [m/yr] </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vy_obs" type="',class(obj.vy_obs),'" default="',convert2str(obj.vy_obs),'">','     <section name="inversion" />','     <help> observed velocity y component [m/yr]  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vel_obs" type="',class(obj.vel_obs),'" default="',convert2str(obj.vel_obs),'">','     <section name="inversion" />','     <help> observed velocity magnitude [m/yr] </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="thickness_obs" type="',class(obj.thickness_obs),'" default="',convert2str(obj.thickness_obs),'">','     <section name="inversion" />','     <help> observed thickness [m]) </help>','  </parameter>');
-               
-            fprintf(fid,'%s\n%s\n','</frame>');    
-            
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Available cost functions">','<section name="inversion" />');                    
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAbsVelMisfit" type="','string','" default="','101','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceRelVelMisfit" type="','string','" default="','102','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVelMisfit" type="','string','" default="','103','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-                
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVxVyMisfit" type="','string','" default="','104','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAverageVelMisfit" type="','string','" default="','105','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsMisfit" type="','string','" default="','106','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-               
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="DragCoefficientAbsGradient" type="','string','" default="','107','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="RheologyBbarAbsGradient" type="','string','" default="','108','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsGradient" type="','string','" default="','109','">','     <section name="inversion" />','     <help> </help>','  </parameter>');
-               
-            fprintf(fid,'%s\n%s\n','</frame>');    
-        
-        end % }}}       
+		function createxml(obj,fid) % {{{
+			fprintf(fid, '<!-- inversion -->\n');            
+
+			% inversion parameters
+			fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="inversion parameters">','<section name="inversion" />');                    
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="iscontrol" type="',class(obj.iscontrol),'" default="',convert2str(obj.iscontrol),'">','     <section name="inversion" />','     <help> is inversion activated? </help>','  </parameter>');
+
+			% incompleteadjoing drop-down (0 or 1)
+			fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="incomplete_adjoint" type="alternative" optional="false">','     <section name="inversion" />','     <help> 1: linear viscosity, 0: non-linear viscosity </help>');
+			fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
+			fprintf(fid,'%s\n%s\n','       <option value="1" type="string" default="false"> </option>','</parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="control_parameters" type="',class(obj.control_parameters),'" default="',convert2str(obj.control_parameters),'">','     <section name="inversion" />','     <help> ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''} </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="nsteps" type="',class(obj.nsteps),'" default="',convert2str(obj.nsteps),'">','     <section name="inversion" />','     <help> number of optimization searches </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions" type="',class(obj.cost_functions),'" default="',convert2str(obj.cost_functions),'">','     <section name="inversion" />','     <help> indicate the type of response for each optimization step  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions_coefficients" type="',class(obj.cost_functions_coefficients),'" default="',convert2str(obj.cost_functions_coefficients),'">','     <section name="inversion" />','     <help> cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_function_threshold" type="',class(obj.cost_function_threshold),'" default="',convert2str(obj.cost_function_threshold),'">','     <section name="inversion" />','     <help> misfit convergence criterion. Default is 1%, NaN if not applied </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter_per_step" type="',class(obj.maxiter_per_step),'" default="',convert2str(obj.maxiter_per_step),'">','     <section name="inversion" />','     <help> maximum iterations during each optimization step  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="gradient_scaling" type="',class(obj.gradient_scaling),'" default="',convert2str(obj.gradient_scaling),'">','     <section name="inversion" />','     <help> scaling factor on gradient direction during optimization, for each optimization step </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="step_threshold" type="',class(obj.step_threshold),'" default="',convert2str(obj.step_threshold),'">','     <section name="inversion" />','     <help> decrease threshold for misfit, default is 30% </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="min_parameters" type="',class(obj.min_parameters),'" default="',convert2str(obj.min_parameters),'">','     <section name="inversion" />','     <help> absolute minimum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="max_parameters" type="',class(obj.max_parameters),'" default="',convert2str(obj.max_parameters),'">','     <section name="inversion" />','     <help> absolute maximum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vx_obs" type="',class(obj.vx_obs),'" default="',convert2str(obj.vx_obs),'">','     <section name="inversion" />','     <help> observed velocity x component [m/yr] </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vy_obs" type="',class(obj.vy_obs),'" default="',convert2str(obj.vy_obs),'">','     <section name="inversion" />','     <help> observed velocity y component [m/yr]  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vel_obs" type="',class(obj.vel_obs),'" default="',convert2str(obj.vel_obs),'">','     <section name="inversion" />','     <help> observed velocity magnitude [m/yr] </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="thickness_obs" type="',class(obj.thickness_obs),'" default="',convert2str(obj.thickness_obs),'">','     <section name="inversion" />','     <help> observed thickness [m]) </help>','  </parameter>');
+
+			fprintf(fid,'%s\n%s\n','</frame>');    
+
+			fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Available cost functions">','<section name="inversion" />');                    
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAbsVelMisfit" type="','string','" default="','101','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceRelVelMisfit" type="','string','" default="','102','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVelMisfit" type="','string','" default="','103','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVxVyMisfit" type="','string','" default="','104','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAverageVelMisfit" type="','string','" default="','105','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsMisfit" type="','string','" default="','106','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="DragCoefficientAbsGradient" type="','string','" default="','107','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="RheologyBbarAbsGradient" type="','string','" default="','108','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsGradient" type="','string','" default="','109','">','     <section name="inversion" />','     <help> </help>','  </parameter>');
+
+			fprintf(fid,'%s\n%s\n','</frame>');    
+
+		end % }}}       
 		function obj = inversion(varargin) % {{{
 			switch nargin
 				case 0
@@ -195,7 +195,7 @@
 			WriteData(fid,'object',obj,'fieldname','incomplete_adjoint','format','Boolean');
 			if ~obj.iscontrol, return; end
 			WriteData(fid,'object',obj,'fieldname','nsteps','format','Integer');
-			WriteData(fid,'object',obj,'fieldname','maxiter_per_step','format','DoubleMat','mattype',3);
+			WriteData(fid,'object',obj,'fieldname','maxiter_per_step','format','IntMat','mattype',3);
 			WriteData(fid,'object',obj,'fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
 			WriteData(fid,'object',obj,'fieldname','gradient_scaling','format','DoubleMat','mattype',3);
 			WriteData(fid,'object',obj,'fieldname','cost_function_threshold','format','Double');