Index: /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 2210)
+++ /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 2211)
@@ -80,4 +80,10 @@
 		param= new Param(count,"tolx",DOUBLE);
 		param->SetDouble(iomodel->tolx);
+		parameters->AddObject(param);
+
+		/*eps_cm: */
+		count++;
+		param= new Param(count,"eps_cm",DOUBLE);
+		param->SetDouble(iomodel->eps_cm);
 		parameters->AddObject(param);
 
Index: /issm/trunk/src/c/ModelProcessorx/IoModel.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/IoModel.cpp	(revision 2210)
+++ /issm/trunk/src/c/ModelProcessorx/IoModel.cpp	(revision 2211)
@@ -113,4 +113,5 @@
 	iomodel->artificial_diffusivity=0;
 	iomodel->nsteps=0;
+	iomodel->eps_cm=0;
 	iomodel->tolx=0;
 	iomodel->maxiter=NULL;
@@ -327,4 +328,5 @@
 	IoModelFetchData((void**)&iomodel->artificial_diffusivity,NULL,NULL,iomodel_handle,"artificial_diffusivity","Integer",NULL);
 	IoModelFetchData((void**)&iomodel->nsteps,NULL,NULL,iomodel_handle,"nsteps","Integer",NULL);
+	IoModelFetchData((void**)&iomodel->eps_cm,NULL,NULL,iomodel_handle,"eps_cm","Scalar",NULL);
 	IoModelFetchData((void**)&iomodel->tolx,NULL,NULL,iomodel_handle,"tolx","Scalar",NULL);
 	IoModelFetchData((void**)&iomodel->mincontrolconstraint,NULL,NULL,iomodel_handle,"mincontrolconstraint","Scalar",NULL);
Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 2210)
+++ /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 2211)
@@ -111,4 +111,5 @@
 	int     artificial_diffusivity;
 	int     nsteps;
+	double  eps_cm;
 	double  tolx;
 	double* maxiter;
Index: /issm/trunk/src/c/parallel/control_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/control_core.cpp	(revision 2210)
+++ /issm/trunk/src/c/parallel/control_core.cpp	(revision 2211)
@@ -30,4 +30,5 @@
 	double* maxiter=NULL;
 	double* cmjump=NULL;
+	double  eps_cm;
 	double  tolx;
 	double* param_g=NULL;
@@ -48,4 +49,5 @@
 	int sub_analysis_type;
 	int debug=0;
+	int convergence=0;
 	int numberofnodes;
 
@@ -61,4 +63,5 @@
 	model->FindParam(&maxiter,"maxiter");
 	model->FindParam(&cmjump,"cmjump");
+	model->FindParam(&eps_cm,"eps_cm");
 	model->FindParam(&tolx,"tolx");
 	model->FindParam(&mincontrolconstraint,"mincontrolconstraint");
@@ -124,4 +127,27 @@
 		xfree((void**)&grad_g_double);
 
+		/*Has convergence been reached?*/
+		if (!isnan(eps_cm)){
+			i=n-2;
+			//go through the previous misfits(starting from n-2)
+			while(i>=0){
+				if (fit[i]==fit[n]){
+					//convergence test only if we have the same misfits
+					if ((J[i]-J[n])/J[n] <= eps_cm){
+						//convergence if convergence criteria fullfilled
+						convergence=1;
+						_printf_("%s%g%s%g\n","      Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],"<",eps_cm);
+					}
+					else{
+						_printf_("%s%g%s%g\n","      Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],">",eps_cm);
+					}
+					break;
+				}
+				i=i-1;
+			}
+		}
+		//stop if convergence has been reached
+		if(convergence) break;
+
 		//some temporary saving
 		/*if (((n+1)%5)==0){
Index: /issm/trunk/src/m/classes/@model/model.m
===================================================================
--- /issm/trunk/src/m/classes/@model/model.m	(revision 2210)
+++ /issm/trunk/src/m/classes/@model/model.m	(revision 2211)
@@ -178,4 +178,5 @@
 	md.nsteps=0;
 	md.maxiter=[];
+	md.eps_cm=0;
 	md.tolx=0;
 	md.optscal=[];
Index: /issm/trunk/src/m/classes/@model/setdefaultparameters.m
===================================================================
--- /issm/trunk/src/m/classes/@model/setdefaultparameters.m	(revision 2210)
+++ /issm/trunk/src/m/classes/@model/setdefaultparameters.m	(revision 2211)
@@ -190,5 +190,10 @@
 
 %cmjump
-md.cmjump=.7*ones(md.nsteps,1); %30 % decrement.
+md.cmjump=.7*ones(md.nsteps,1); %30 per cent decrement.
+
+%eps_cm is a criteria to stop the control methods.
+%if J[n]-J[n-1]/J[n] < criteria, the control run stops
+%NaN if not applied
+md.eps_cm=0.01; %1 percent
 
 %minimum velocity to avoid the misfit to be singular
Index: /issm/trunk/src/m/classes/public/display/displaycontrol.m
===================================================================
--- /issm/trunk/src/m/classes/public/display/displaycontrol.m	(revision 2210)
+++ /issm/trunk/src/m/classes/public/display/displaycontrol.m	(revision 2211)
@@ -13,4 +13,5 @@
 	fielddisplay(md,'control_type','parameter where inverse control is carried out; ex: ''drag'', or ''B'')');
 	fielddisplay(md,'nsteps','number of optimization searches');
+	fielddisplay(md,'eps_cm','misfit convergence criterion. Default is 1%, NaN if not applied');
 	fielddisplay(md,'optscal','scaling factor on gradient direction during optimization, for each optimization step');
 	fielddisplay(md,'fit','''absolute: 0'', ''relative: 1'', or ''logarithmic: 2''. default is ''absolute: 0'', for each optimization steps');
Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m
===================================================================
--- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 2210)
+++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 2211)
@@ -351,6 +351,6 @@
 
 	%LENGTH CONTROL FIELDS
-	if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps)
-		disp('maxiter, optscal and fit must have the length specified by nsteps')
+	if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps | length(md.cmjump)~=md.nsteps)
+		disp('maxiter, optscal, fit and cmjump must have the length specified by nsteps')
 		bool=0;return;
 	end
Index: /issm/trunk/src/m/classes/public/marshall.m
===================================================================
--- /issm/trunk/src/m/classes/public/marshall.m	(revision 2210)
+++ /issm/trunk/src/m/classes/public/marshall.m	(revision 2211)
@@ -109,4 +109,5 @@
 WriteData(fid,md.artificial_diffusivity,'Integer','artificial_diffusivity');
 WriteData(fid,md.nsteps,'Integer','nsteps');
+WriteData(fid,md.eps_cm,'Scalar','eps_cm');
 WriteData(fid,md.tolx,'Scalar','tolx');
 WriteData(fid,md.maxiter,'Mat','maxiter');
Index: /issm/trunk/src/m/solutions/cielo/control_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/control_core.m	(revision 2210)
+++ /issm/trunk/src/m/solutions/cielo/control_core.m	(revision 2211)
@@ -16,4 +16,8 @@
 dim=model.parameters.dim;
 isstokes=model.parameters.isstokes;
+
+%recover other parameters
+eps_cm=model.parameters.eps_cm;
+fit=model.parameters.fit;
 
 %initialize control parameters
@@ -71,10 +75,29 @@
 	param_g=ControlConstrain(param_g,model.parameters);
 
-	%visualize direction.
-	if model.parameters.plot,
-		plot_newdistribution;
+	disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
+
+	%Has convergence been reached?
+	convergence=0;
+	if ~isnan(eps_cm),
+		i=n-2;
+		%go through the previous misfits(starting from n-2)
+		while (i>=1),
+			if (fit(i)==fit(n)),
+				%convergence test only if we have the same misfits
+				if ((c(i).J-c(n).J)/c(n).J <= eps_cm),
+					%convergence if convergence criteria fullfilled
+					convergence=1;
+					displaystring(debug,'\n%s%g%s%g\n','      Convergence criterion: dJ/J = ',(c(i).J-c(n).J)/c(n).J,'<',eps_cm);
+				else
+					displaystring(debug,'\n%s%g%s%g\n','      Convergence criterion: dJ/J = ',(c(i).J-c(n).J)/c(n).J,'>',eps_cm);
+				end
+				break;
+			end
+			i=i-1;                                                                                                                                         
+		end
 	end
+	%stop if convergence has been reached                                                                                                               
+	if (convergence), break; end
 
-	disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
 end