Index: /issm/trunk/test/Miscellaneous/Regularization/DomainOutline.exp
===================================================================
--- /issm/trunk/test/Miscellaneous/Regularization/DomainOutline.exp	(revision 6047)
+++ /issm/trunk/test/Miscellaneous/Regularization/DomainOutline.exp	(revision 6047)
@@ -0,0 +1,10 @@
+## Name:domainoutline
+## Icon:0
+# Points Count  Value
+5 1.
+# X pos Y pos
+0 0
+100000 0
+100000 100000
+0 100000
+0 0
Index: /issm/trunk/test/Miscellaneous/Regularization/Square.par
===================================================================
--- /issm/trunk/test/Miscellaneous/Regularization/Square.par	(revision 6047)
+++ /issm/trunk/test/Miscellaneous/Regularization/Square.par	(revision 6047)
@@ -0,0 +1,34 @@
+%Start defining model parameters here
+hmin=1000;
+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+10;
+md.surface=md.bed+md.thickness;
+
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag_coefficient=30*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.observed_temperature=(273-20)*ones(md.numberofgrids,1);
+
+md.rheology_B=paterson(md.observed_temperature);
+md.rheology_n=3*ones(md.numberofelements,1);
+
+%Deal with boundary conditions:
+md=SetIceSheetBC(md);
+pos=find(md.gridonboundary);
+md.spcvelocity(pos,1)=1;
+md.spcvelocity(pos,2)=1;
+md.spcvelocity(pos,3)=1;
+md.spcvelocity(pos,4)=0;
+md.spcvelocity(pos,5)=200;
+md.spcvelocity(pos,6)=0;
+
+md.eps_res=10^-5;
Index: /issm/trunk/test/Miscellaneous/Regularization/runme.m
===================================================================
--- /issm/trunk/test/Miscellaneous/Regularization/runme.m	(revision 6047)
+++ /issm/trunk/test/Miscellaneous/Regularization/runme.m	(revision 6047)
@@ -0,0 +1,116 @@
+steps=[0];
+
+if ismember(1,steps),
+	md=mesh(model,'DomainOutline.exp',2000);
+	md=geography(md,'','');
+	md=parameterize(md,'Square.par');
+	md=setelementstype(md,'macayeal','all');
+
+	%Run 1: given basal drag
+	centerx=50000;
+	centery=50000;
+	radius=50000;
+	pos=find(sqrt((md.x-centerx).^2+(md.y-centery).^2)<radius); 
+	md.drag_coefficient(pos)=md.drag_coefficient(pos)+50*cos(pi/(2*radius)*(sqrt((md.y(pos)-centery).^2+(md.x(pos)-centerx).^2)));
+
+	md.cluster=eval(oshostname);
+	md=solve(md,'analysis_type',DiagnosticSolutionEnum);
+	md.vx_obs=PatchToVec(md.results.DiagnosticSolution.Vx);
+	md.vy_obs=PatchToVec(md.results.DiagnosticSolution.Vy);
+	md.dummy=md.drag_coefficient;
+end
+
+if ismember(2,steps),
+	md.drag_coefficient(:)=30;
+	md.control_analysis=1;
+	md.nsteps=20;
+	md.cm_responses=SurfaceLogVelMisfitEnum*ones(md.nsteps,1);
+	md.optscal=30*ones(md.nsteps,1);
+	md.maxiter=15*ones(md.nsteps,1);
+	md.cm_jump=0.99*ones(md.nsteps,1);
+	md.cm_noisedmp=0;
+	md.weights=1*ones(md.numberofgrids,1);
+	md.eps_cm=NaN;
+	md=solve(md,'analysis_type',DiagnosticSolutionEnum);
+	md.results.DiagnosticSolution_NO=md.results.DiagnosticSolution;
+end
+
+if ismember(3,steps),
+	md.drag_coefficient(:)=30;
+	md.control_analysis=1;
+	md.nsteps=20;
+	md.cm_responses=SurfaceLogVelMisfitEnum*ones(md.nsteps,1);
+	md.optscal=30*ones(md.nsteps,1);
+	md.maxiter=15*ones(md.nsteps,1);
+	md.cm_jump=0.99*ones(md.nsteps,1);
+	md.cm_noisedmp=0.6*10^-9; %0.5 too small
+	md.weights=1*ones(md.numberofgrids,1);
+	md.eps_cm=NaN;
+	md=solve(md,'analysis_type',DiagnosticSolutionEnum);
+	md.results.DiagnosticSolution_YES1=md.results.DiagnosticSolution;
+end
+
+if ismember(4,steps),
+	md.drag_coefficient(:)=30;
+	md.control_analysis=1;
+	md.nsteps=20;
+	md.cm_responses=SurfaceLogVelMisfitEnum*ones(md.nsteps,1);
+	md.optscal=30*ones(md.nsteps,1);
+	md.maxiter=15*ones(md.nsteps,1);
+	md.cm_jump=0.99*ones(md.nsteps,1);
+	md.cm_noisedmp=5*10^-9;
+	md.weights=1*ones(md.numberofgrids,1);
+	md.eps_cm=NaN;
+	md=solve(md,'analysis_type',DiagnosticSolutionEnum);
+	md.results.DiagnosticSolution_YES2=md.results.DiagnosticSolution;
+end
+
+if ismember(5,steps),
+	md.drag_coefficient(:)=30;
+	md.control_analysis=1;
+	md.nsteps=20;
+	md.cm_responses=SurfaceLogVelMisfitEnum*ones(md.nsteps,1);
+	md.optscal=30*ones(md.nsteps,1);
+	md.maxiter=15*ones(md.nsteps,1);
+	md.cm_jump=0.99*ones(md.nsteps,1);
+	md.cm_noisedmp=8*10^-8;
+	md.weights=1*ones(md.numberofgrids,1);
+	md.eps_cm=NaN;
+	md=solve(md,'analysis_type',DiagnosticSolutionEnum);
+	md.results.DiagnosticSolution_YES3=md.results.DiagnosticSolution;
+end
+
+save model.mdl md
+
+printflag=true;
+plotmodel(md,'data',md.results.DiagnosticSolution_NO.DragCoefficient,'caxis#all',[10 90],'unit','km','axis','equal','xlim',[0 100],'ylim',[0 100]); set(gcf,'Color','w');
+if printflag, printmodel('SquareTiknonov0','png','margin','on','frame','off','marginsize',25,'resolution',1,'hardcopy','off'); end
+plotmodel(md,'data',md.results.DiagnosticSolution_YES1.DragCoefficient,'caxis#all',[10 90],'unit','km','axis','equal','xlim',[0 100],'ylim',[0 100]); set(gcf,'Color','w');
+if printflag, printmodel('SquareTiknonov1','png','margin','on','frame','off','marginsize',25,'resolution',1,'hardcopy','off'); end
+plotmodel(md,'data',md.results.DiagnosticSolution_YES2.DragCoefficient,'caxis#all',[10 90],'unit','km','axis','equal','xlim',[0 100],'ylim',[0 100]); set(gcf,'Color','w');
+if printflag, printmodel('SquareTiknonov2','png','margin','on','frame','off','marginsize',25,'resolution',1,'hardcopy','off'); end
+plotmodel(md,'data',md.results.DiagnosticSolution_YES3.DragCoefficient,'caxis#all',[10 90],'unit','km','axis','equal','xlim',[0 100],'ylim',[0 100]); set(gcf,'Color','w');
+if printflag, printmodel('SquareTiknonov3','png','margin','on','frame','off','marginsize',25,'resolution',1,'hardcopy','off'); end
+plotmodel(md,'data',md.dummy,'caxis#all',[10 90],'unit','km','axis','equal','xlim',[0 100],'ylim',[0 100]); set(gcf,'Color','w');
+if printflag, printmodel('SquareRealCoeff','png','margin','on','frame','off','marginsize',25,'resolution',1,'hardcopy','off'); end
+
+if false,
+	plotmodel(md,...
+		'data',PatchToVec(md.results.DiagnosticSolution_NO.Vy)-md.vy_obs,...
+		'data',PatchToVec(md.results.DiagnosticSolution_YES1.Vy)-md.vy_obs,...
+		'data',PatchToVec(md.results.DiagnosticSolution_YES2.Vy)-md.vy_obs,...
+		'data',PatchToVec(md.results.DiagnosticSolution_YES3.Vy)-md.vy_obs,...
+		'caxis#all',[-0.3 0.3]);
+end
+if false
+	plotmodel(md,...
+		'data',md.results.DiagnosticSolution_NO.DragCoefficient,...
+		'data',md.results.DiagnosticSolution_YES1.DragCoefficient,...
+		'data',md.results.DiagnosticSolution_YES2.DragCoefficient,...
+		'data',md.results.DiagnosticSolution_YES3.DragCoefficient,...
+		'data',md.dummy,...
+		'nlines',2,'ncols',3,'caxis#all',[10 90]);
+end
+
+%plotmodel(md,'data',md.results.DiagnosticSolution_NO.Vx,'data',md.results.DiagnosticSolution_YES1.Vx,'data',md.results.DiagnosticSolution_YES2.Vx,'data',md.vx_obs);
+%plotmodel(md,'data',PatchToVec(md.results.DiagnosticSolution_NO.Vy)-md.vy_obs,'data',PatchToVec(md.results.DiagnosticSolution_YES1.Vy)-md.vy_obs,'data',PatchToVec(md.results.DiagnosticSolution_YES2.Vy)-md.vy_obs);