Index: sm/trunk/test/Validation/ControlMethods/runme.m
===================================================================
--- /issm/trunk/test/Validation/ControlMethods/runme.m	(revision 409)
+++ 	(revision )
@@ -1,76 +1,0 @@
-step=[2];
-division=2;
-density=70000;
-startup
-
-%step1: initial velocity -> observed velocity
-if ismember(1,step),
-	disp(sprintf('\nSTEP 1 : Build observed velocity\n'));
-
-	%create model
-	md=model;
-	md=mesh(md,'DomainOutline.exp',density);
-	md=geography(md,'all','');
-	md=parameterize(md,'Square.par');
-	md=setelementstype(md,'macayeal','all');
-
-	%tweak B
-	maxx=max(md.x); minx=min(md.x); Lx=maxx-minx;
-	maxy=max(md.y); miny=min(md.y); Ly=maxy-miny;
-	temperature=linspace(213,263,division^2);
-	temperature=temperature(randperm(division^2));
-	count=1;
-	for i=1:division
-		for j=1:division
-			md.observed_temperature(find(md.x>=(minx+(i-1)*Lx/division) & md.x<=(minx+i*Lx/division) & md.y>=(miny+(j-1)*Ly/division) & md.y<=(miny+j*Ly/division)))=temperature(count);
-			md.B=paterson(md.observed_temperature);
-			count=count+1;
-		end
-	end
-
-	%diagnostic with this B
-	md=solve(md,'diagnostic','cielo');
-	%md.cluster='wilkes'; md=solve(md,'diagnostic_horiz','cielo');
-
-	%save observed velocities and exact B
-	vx_obs=md.vx; vy_obs=md.vy; vel_obs=md.vel;
-	B=md.B;
-	save ObservedVelocities vx_obs vy_obs vel_obs
-	save ExactB B
-	clear md
-end
-
-%step2: control method
-if ismember(2,step),
-	disp(sprintf('\nSTEP 2 : control method\n'));
-
-	%create model
-	md=model;
-	md=mesh(md,'DomainOutline.exp',density);
-	md=geography(md,'all','');
-	md=parameterize(md,'Square.par');
-	md=setelementstype(md,'macayeal','all');
-
-	%plug computed velocities
-	load ObservedVelocities
-	md.vx_obs=vx_obs; md.vy_obs=vy_obs; md.vel_obs=vel_obs;
-
-	%control method
-	md.nsteps=1;
-	md.control_type='B';
-	md.mincontrolconstraint=0;
-	md.maxcontrolconstraint=1.2*10^9;
-	md.plot=0;
-	md.debug=1;
-	md.optscal=10^8*ones(md.nsteps,1); md.optscal(1)=2*10^8;
-	md.fit=0*ones(md.nsteps,1);
-	md.maxiter=1*ones(md.nsteps,1);
-	md.eps_rel=0.01;
-	md.eps_abs=NaN;
-
-	md.np=1;
-md.cluster='wilkes';
-%	md=solve(md,'control','ice');
-md=solve(md,'control','cielo');
-%md=solve(md,'control','macayeal');
-end
Index: /issm/trunk/test/Validation/ControlMethods/runmeB.m
===================================================================
--- /issm/trunk/test/Validation/ControlMethods/runmeB.m	(revision 410)
+++ /issm/trunk/test/Validation/ControlMethods/runmeB.m	(revision 410)
@@ -0,0 +1,77 @@
+step=[2];
+division=2;
+density=100000;
+startup
+
+%step1: initial velocity -> observed velocity
+if ismember(1,step),
+	disp(sprintf('\nSTEP 1 : Build observed velocity\n'));
+
+	%create model
+	md=model;
+	md=mesh(md,'DomainOutline.exp',density);
+	md=geography(md,'all','');
+	md=parameterize(md,'Square.par');
+	md=setelementstype(md,'macayeal','all');
+
+	%tweak B
+	maxx=max(md.x); minx=min(md.x); Lx=maxx-minx;
+	maxy=max(md.y); miny=min(md.y); Ly=maxy-miny;
+	temperature=linspace(213,263,division^2);
+	temperature=temperature(randperm(division^2));
+	count=1;
+	for i=1:division
+		for j=1:division
+			md.observed_temperature(find(md.x>=(minx+(i-1)*Lx/division) & md.x<=(minx+i*Lx/division) & md.y>=(miny+(j-1)*Ly/division) & md.y<=(miny+j*Ly/division)))=temperature(count);
+			md.B=paterson(md.observed_temperature);
+			count=count+1;
+		end
+	end
+
+	%diagnostic with this B
+	md.cluster='wilkes'; 
+	md=solve(md,'diagnostic','cielo');
+
+	%save observed velocities and exact B
+	vx_obs=md.vx; vy_obs=md.vy; vel_obs=md.vel;
+	B=md.B;
+	save ObservedVelocities vx_obs vy_obs vel_obs
+	save ExactB B
+	clear md
+end
+
+%step2: control method
+if ismember(2,step),
+	disp(sprintf('\nSTEP 2 : control method\n'));
+
+	%create model
+	md=model;
+	md=mesh(md,'DomainOutline.exp',density);
+	md=geography(md,'all','');
+	md=parameterize(md,'Square.par');
+	md=setelementstype(md,'macayeal','all');
+
+	%plug computed velocities
+	load ObservedVelocities
+	md.vx_obs=vx_obs; md.vy_obs=vy_obs; md.vel_obs=vel_obs;
+
+	%control method
+	md.nsteps=5;
+	md.control_type='B';
+	md.mincontrolconstraint=10^7;
+	md.maxcontrolconstraint=1.2*10^9;
+	md.plot=1;
+	md.debug=0;
+	md.optscal=10^8*ones(md.nsteps,1); md.optscal(1)=2*10^8;
+	md.fit=0*ones(md.nsteps,1);
+	md.maxiter=10*ones(md.nsteps,1);
+	md.eps_rel=0.0001;
+	md.eps_abs=NaN;
+	md.viscosity_overshoot=0;
+
+	md.np=8;
+	%md.cluster='wilkes';
+	%md=solve(md,'control','ice');
+	md=solve(md,'control','cielo');
+	%md=solve(md,'control','macayeal');
+end
Index: /issm/trunk/test/Validation/ControlMethods/runmedrag.m
===================================================================
--- /issm/trunk/test/Validation/ControlMethods/runmedrag.m	(revision 410)
+++ /issm/trunk/test/Validation/ControlMethods/runmedrag.m	(revision 410)
@@ -0,0 +1,74 @@
+step=[2];
+division=2;
+density=70000;
+startup
+
+%step1: initial velocity -> observed velocity
+if ismember(1,step),
+	disp(sprintf('\nSTEP 1 : Build observed velocity\n'));
+
+	%create model
+	md=model;
+	md=mesh(md,'DomainOutline.exp',density);
+	md=geography(md,'','');
+	md=parameterize(md,'Squaredrag.par');
+	md=setelementstype(md,'macayeal','all');
+
+	%tweak drag
+	maxx=max(md.x); minx=min(md.x); Lx=maxx-minx;
+	maxy=max(md.y); miny=min(md.y); Ly=maxy-miny;
+	drag=linspace(10,50,division^2);
+	drag=drag(randperm(division^2));
+	count=1;
+	for i=1:division
+		for j=1:division
+			md.drag(find(md.x>=(minx+(i-1)*Lx/division) & md.x<=(minx+i*Lx/division) & md.y>=(miny+(j-1)*Ly/division) & md.y<=(miny+j*Ly/division)))=drag(count);
+			count=count+1;
+		end
+	end
+
+	%diagnostic with this drag
+	md=solve(md,'diagnostic','cielo');
+	%md.cluster='wilkes'; md=solve(md,'diagnostic_horiz','cielo');
+
+	%save observed velocities and exact B
+	vx_obs=md.vx; vy_obs=md.vy; vel_obs=md.vel;
+	drag=md.drag;
+	save ObservedVelocities vx_obs vy_obs vel_obs
+	save Exactdrag drag
+	clear md
+end
+
+%step2: control method
+if ismember(2,step),
+	disp(sprintf('\nSTEP 2 : control method\n'));
+
+	%create model
+	md=model;
+	md=mesh(md,'DomainOutline.exp',density);
+	md=geography(md,'','');
+	md=parameterize(md,'Squaredrag.par');
+	md=setelementstype(md,'macayeal','all');
+
+	%plug computed velocities
+	load ObservedVelocities
+	md.vx_obs=vx_obs; md.vy_obs=vy_obs; md.vel_obs=vel_obs;
+
+	%control method
+	md.nsteps=10;
+	md.control_type='drag';
+	md.mincontrolconstraint=0;
+	md.maxcontrolconstraint=100;
+	md.plot=1;
+	md.debug=0;
+	md.optscal=180*ones(md.nsteps,1);
+	md.fit=2*ones(md.nsteps,1);
+	md.maxiter=10*ones(md.nsteps,1);
+	md.eps_rel=0.001;
+	md.eps_abs=NaN;
+
+%md.cluster='wilkes';
+%md=solve(md,'control','ice');
+md=solve(md,'control','cielo');
+%md=solve(md,'control','macayeal');
+end