Index: /issm/trunk/test/Validation/ControlMethods/Square.par
===================================================================
--- /issm/trunk/test/Validation/ControlMethods/Square.par	(revision 4849)
+++ /issm/trunk/test/Validation/ControlMethods/Square.par	(revision 4850)
@@ -34,5 +34,5 @@
 %Deal with boundary conditions:
 md=SetIceShelfBC(md,'Front.exp');
-md.dirichletvalues_diag(:,2)=500;
+md.spcvelocity(:,5)=500;
 
 %Parallel options
Index: /issm/trunk/test/Validation/ControlMethods/Squaredrag.par
===================================================================
--- /issm/trunk/test/Validation/ControlMethods/Squaredrag.par	(revision 4850)
+++ /issm/trunk/test/Validation/ControlMethods/Squaredrag.par	(revision 4850)
@@ -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_coefficient=200*ones(md.numberofgrids,1); %q=1.
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag_coefficient(md.elements(pos,:))=0;
+md.drag_p=ones(md.numberofelements,1);
+md.drag_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.rheology_B=paterson(md.observed_temperature);
+md.rheology_n=3*ones(md.numberofelements,1);
+
+%Deal with boundary conditions:
+md=SetMarineIceSheetBC(md,'Front.exp');
+md.spcvelocity(:,5)=500;
+
+%Parallel options
+md.np=3;
+md.time=50;
+md.waitonlock=1;
Index: /issm/trunk/test/Validation/ControlMethods/runmeB.m
===================================================================
--- /issm/trunk/test/Validation/ControlMethods/runmeB.m	(revision 4849)
+++ /issm/trunk/test/Validation/ControlMethods/runmeB.m	(revision 4850)
@@ -24,5 +24,5 @@
 		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);
+			md.rheology_B=paterson(md.observed_temperature);
 			count=count+1;
 		end
@@ -31,10 +31,13 @@
 	%diagnostic with this B
 	md.cluster='wilkes'; 
-	md=solve(md,'analysis_type','diagnostic');
+	md=solve(md,'analysis_type',DiagnosticSolutionEnum);
 
 	%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
+	vx_obs=zeros(md.numberofgrids,1);
+	vx_obs(md.results.DiagnosticSolution.Vx.index)=md.results.DiagnosticSolution.Vx.value;
+	vy_obs=zeros(md.numberofgrids,1);
+	vy_obs(md.results.DiagnosticSolution.Vy.index)=md.results.DiagnosticSolution.Vy.value;
+	B=md.rheology_B;
+	save ObservedVelocities vx_obs vy_obs
 	save ExactB B
 	clear md
@@ -54,15 +57,17 @@
 	%plug computed velocities
 	load ObservedVelocities
-	md.vx_obs=vx_obs; md.vy_obs=vy_obs; md.vel_obs=vel_obs;
+	md.vx_obs=vx_obs; md.vy_obs=vy_obs;
 
 	%control method
 	md.nsteps=5;
-	md.control_type='B';
-	md.mincontrolconstraint=10^7;
-	md.maxcontrolconstraint=1.2*10^9;
+	md.control_type='rheology_B';
+	md.cm_min=10^7;
+	md.cm_max=1.2*10^9;
+	md.cm_jump=0.99*ones(md.nsteps,1);
 	md.plot=1;
-	md.debug=0;
-	md.optscal=10^8*ones(md.nsteps,1); md.optscal(1)=2*10^8;
+	md.verbose=0;
+	md.optscal=2*10^8*ones(md.nsteps,1); md.optscal(1)=5*10^8;
 	md.fit=0*ones(md.nsteps,1);
+	md.weights=1*ones(md.numberofgrids,1);
 	md.maxiter=10*ones(md.nsteps,1);
 	md.eps_rel=0.0001;
@@ -71,7 +76,5 @@
 
 	md.np=8;
-	%md.cluster='wilkes';
-	%md=solve(md,'analysis_type','control');
-	md=solve(md,'analysis_type','control');
-	%md=solve(md,'analysis_type','control','package','macayeal');
+	md.control_analysis=1;
+	md=solve(md,'analysis_type',DiagnosticSolutionEnum);
 end
Index: /issm/trunk/test/Validation/ControlMethods/runmedrag.m
===================================================================
--- /issm/trunk/test/Validation/ControlMethods/runmedrag.m	(revision 4849)
+++ /issm/trunk/test/Validation/ControlMethods/runmedrag.m	(revision 4850)
@@ -1,3 +1,3 @@
-step=[2];
+step=[1 2];
 division=2;
 density=70000;
@@ -23,5 +23,5 @@
 	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);
+			md.drag_coefficient(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
@@ -29,11 +29,13 @@
 
 	%diagnostic with this drag
-	md=solve(md,'analysis_type','diagnostic');
-	%md.cluster='wilkes'; md=solve(md,'analysis_type','diagnostic_horiz');
+	md=solve(md,'analysis_type',DiagnosticSolutionEnum);
 
 	%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
+	vx_obs=zeros(md.numberofgrids,1);
+	vx_obs(md.results.DiagnosticSolution.Vx.index)=md.results.DiagnosticSolution.Vx.value;
+	vy_obs=zeros(md.numberofgrids,1);
+	vy_obs(md.results.DiagnosticSolution.Vy.index)=md.results.DiagnosticSolution.Vy.value;
+	drag=md.drag_coefficient;
+	save ObservedVelocities vx_obs vy_obs
 	save Exactdrag drag
 	clear md
@@ -53,13 +55,15 @@
 	%plug computed velocities
 	load ObservedVelocities
-	md.vx_obs=vx_obs; md.vy_obs=vy_obs; md.vel_obs=vel_obs;
+	md.vx_obs=vx_obs; md.vy_obs=vy_obs;
 
 	%control method
 	md.nsteps=10;
-	md.control_type='drag';
-	md.mincontrolconstraint=0;
-	md.maxcontrolconstraint=100;
+	md.control_type='drag_coefficient';
+	md.cm_min=0;
+	md.cm_max=100;
+	md.cm_jump=0.99*ones(md.nsteps,1);
 	md.plot=1;
-	md.debug=0;
+	md.verbose=0;
+	md.weights=1*ones(md.numberofgrids,1);
 	md.optscal=180*ones(md.nsteps,1);
 	md.fit=2*ones(md.nsteps,1);
@@ -69,6 +73,5 @@
 
 %md.cluster='wilkes';
-%md=solve(md,'analysis_type','control');
-md=solve(md,'analysis_type','control');
-%md=solve(md,'analysis_type','control','package','macayeal');
+md.control_analysis=1;
+md=solve(md,'analysis_type',DiagnosticSolutionEnum);
 end
