Index: /issm/trunk/test/Validation/ControlMethods/DomainOutline.exp
===================================================================
--- /issm/trunk/test/Validation/ControlMethods/DomainOutline.exp	(revision 161)
+++ /issm/trunk/test/Validation/ControlMethods/DomainOutline.exp	(revision 161)
@@ -0,0 +1,10 @@
+## Name:DomainOutline
+## Icon:0
+# Points Count  Value
+5 1.000000
+# X pos Y pos
+0 0
+1000000 0
+1000000 1000000
+0 1000000
+0 0
Index: /issm/trunk/test/Validation/ControlMethods/Front.exp
===================================================================
--- /issm/trunk/test/Validation/ControlMethods/Front.exp	(revision 161)
+++ /issm/trunk/test/Validation/ControlMethods/Front.exp	(revision 161)
@@ -0,0 +1,10 @@
+## Name:
+## Icon:0
+# Points Count	Value
+5	1.
+# X pos	Y pos
+-1000 990000
+1001000  990000
+1001000  1001000
+-1000  1001000
+-1000 990000
Index: /issm/trunk/test/Validation/ControlMethods/Square.par
===================================================================
--- /issm/trunk/test/Validation/ControlMethods/Square.par	(revision 161)
+++ /issm/trunk/test/Validation/ControlMethods/Square.par	(revision 161)
@@ -0,0 +1,41 @@
+%Start defining model parameters here
+
+%dynamics
+md.dt=1*md.yts; %1 year
+md.ndt=md.dt*10; 
+md.artificial_diffusivity=1;
+
+disp('      creating thickness');
+hmin=300;
+hmax=1000;
+ymin=min(md.y);
+ymax=max(md.y);
+md.thickness=hmax+(hmin-hmax)*(md.y-ymin)/(ymax-ymin);
+md.bed=-md.rho_ice/md.rho_water*md.thickness;
+md.surface=md.bed+md.thickness;
+
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag=200*ones(md.numberofgrids,1); %q=1.
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag(md.elements(pos,:))=0;
+md.p=ones(md.numberofelements,1);
+md.q=ones(md.numberofelements,1);
+md.viscosity_overshoot=0.0;
+
+disp('      creating temperature');
+md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
+
+disp('      creating flow law paramter');
+md.B=paterson(md.observed_temperature);
+md.n=3*ones(md.numberofelements,1);
+
+%Deal with boundary conditions:
+md=SetIceShelfBC(md,'Front.exp');
+md.dirichletvalues_diag(:,2)=500;
+
+%Parallel options
+md.np=8;
+md.time=50;
+md.waitonlock=1;
Index: /issm/trunk/test/Validation/ControlMethods/runme.m
===================================================================
--- /issm/trunk/test/Validation/ControlMethods/runme.m	(revision 161)
+++ /issm/trunk/test/Validation/ControlMethods/runme.m	(revision 161)
@@ -0,0 +1,67 @@
+step=[2];
+division=2;
+
+%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',100000);
+	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');
+
+	%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',100000);
+	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=10;
+	md.control_type='B';
+	md.plot=1;
+	md.debug=1;
+	md.optscal=10^7*ones(md.nsteps,1);
+	md.fit=0*ones(md.nsteps,1);
+	md.maxiter=30*ones(md.nsteps,1);
+
+	%md.cluster='wilkes';
+	md.eps_abs=NaN;
+	md=solve(md,'control','cielo');
+end