Index: /issm/trunk/src/m/solutions/control_core.m
===================================================================
--- /issm/trunk/src/m/solutions/control_core.m	(revision 4536)
+++ /issm/trunk/src/m/solutions/control_core.m	(revision 4537)
@@ -28,4 +28,7 @@
 	options.MaxIter=femmodel.parameters.MaxIter;
 
+	%Initialize misfits with a vector of zeros
+	J=zeros(nsteps,1);
+
 	for n=1:nsteps,
 
@@ -48,13 +51,17 @@
 
 		displaystring(verbose,'\n%s',['      optimizing along gradient direction...']);
-		[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,femmodel,param_g,c(n).grad_g,n,femmodel.parameters);
+		[search_scalar J(n)]=ControlOptimization('objectivefunctionC',0,1,options,femmodel,n,femmodel.parameters);
 
-		displaystring(verbose,'\n%s',['      updating parameter using optimized search scalar...']);
-		param_g=param_g+search_scalar*femmodel.parameters.Optscal(n)*grad_g;
+		displaystring('\n%s',['      updating parameter using optimized search scalar:']);
+		scalar=search_scalar*optscal(n);
+		femmodel.elements=InputAXPY(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,scalar,ControlParameterEnum);
 
-		displaystring(verbose,'\n%s',['      constraining the new distribution...']);
-		param_g=ControlConstrain(param_g,femmodel.parameters);
+		displaystring('\n%s',['      constraning the new distribution...']);
+		[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlConstrain(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,cm_min,cm_max);
 
-		disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
+		displaystring('\n%s',['      save new parameter...']);
+		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,ControlParameterEnum);
+
+		disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(J(n))]);
 
 		%Has convergence been reached?
Index: /issm/trunk/src/m/solutions/controlconvergence.m
===================================================================
--- /issm/trunk/src/m/solutions/controlconvergence.m	(revision 4536)
+++ /issm/trunk/src/m/solutions/controlconvergence.m	(revision 4537)
@@ -1,3 +1,3 @@
-function converged=controlconvergence(J,fit,eps_cm,n)
+function convergence=controlconvergence(J,fit,eps_cm,n)
 %CONTROLCONVERGENCE - determine the convergence of control_core solution
 %
Index: /issm/trunk/src/m/solutions/objectivefunctionC.m
===================================================================
--- /issm/trunk/src/m/solutions/objectivefunctionC.m	(revision 4536)
+++ /issm/trunk/src/m/solutions/objectivefunctionC.m	(revision 4537)
@@ -1,25 +1,41 @@
-function J =objectivefunctionC(search_scalar,models,p_g,grad_g,n,analysis_type,sub_analysis_type);
+function J =objectivefunctionC(search_scalar,femmodel,n);
+%OBJECTIVEFUNCTIONC - objective function that return a parameter for a certain function
 
-%recover active model.
-m=models.active;
+conserve_loads=true;
+%recover some parameters
+optscal=femmodel.parameters.OptScal(n);
+fit=femmodel.parameters.Fit(n);
+control_type=femmodel.parameters.ControlType;
+control_steady=femmodel.parameters.ControlSteady;
+analysis_type=femmodel.parameters.AnalysisType;
+isstokes=femmodel.parameters.IsStokes;
+cm_min=femmodel.parameters.CmMin;
+cm_max=femmodel.parameters.CmMax;
 
-%recover some parameters
-optscal=m.parameters.Optscal(n);
-fit=m.parameters.Fit(n);
-control_type=m.parameters.ControlType;
-analysis_type=m.parameters.AnalysisType;
+%set current configuration
+if isstokes,
+	femmodel=SetCurrentConfiguration(femmodel,DiagnsoticStokesAnalysisEnum);
+else
+	femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
+end
 
-%Update along gradient using scalar supplied by fmincon optimization routine
-parameter=p_g+search_scalar*optscal*grad_g;
+%Use ControlParameterEnum input to  reinitialize our input parameter:
+femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ControlParameterEnum,control_type);
 
-%Plug parameter into inputs
-inputs=add(inputs,m.parameters.ControlType,parameter,'doublevec',1,m.parameters.NumberOfNodes);
+%Use search scalar to shoot parameter in the gradient direction:
+scalar=search_scalar*optscal;
+femmodel.elements=InputAXPY(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,scalar,ControlParameterEnum);
 
-%Run diagnostic with updated parameters.
-u_g=diagnostic_core_nonlinear(m,analysis_type,sub_analysis_type);
+%Constrain:
+[femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlConstrain(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,cm_min,cm_max);
 
-%add velocity to inputs.
-inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.NumberOfDofsPerNode,m.parameters.NumberOfNodes);
+%Run diagnostic with updated inputs:
+if(control_steady==0),
+	femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads);  %true means we conserve loads at each diagnostic run
+else
+	femmodel=diagnostic_core(femmodel);  %We need a 3D velocity!! (vz is required for the next thermal run)
+end
 
-%Compute misfit for this velocity field. 
-J=CostFunction(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type);
+%Compute misfit for this velocity field
+[femmodel.elements,femmodel.loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices, femmodel.loads,femmodel.materials,femmodel.parameters,fit,FitEnum);
+J=CostFunction(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials, femmodel.parameters);