Index: /issm/trunk/scripts/runme.m
===================================================================
--- /issm/trunk/scripts/runme.m	(revision 3206)
+++ /issm/trunk/scripts/runme.m	(revision 3206)
@@ -0,0 +1,151 @@
+%AGU 09: second attempt (after Moca09) at  control methods on entire Antarctica.
+
+%Run steps
+steps=[3];
+
+%steps available: 
+rs=runsteps('Models/','SeaRISE.Greenland.step#');
+rs=addstep(rs,'Mesh generation');
+rs=addstep(rs,'Parameterization');
+rs=addstep(rs,'Control Method drag 1');
+
+%Hard coded parameters
+modeldatapath='./../../Data/';
+%clustername=oshostname();
+clustername='cosmos';
+type='2d';
+
+%Cluster parameters{{{1
+if strcmpi(clustername,'pfe');
+	cluster.name='pfe';      cluster.np=8*70; cluster.queue='long';  cluster.time=12*60;
+elseif strcmpi(clustername,'cosmos')
+	cluster.name='cosmos';   cluster.np=100;  cluster.queue='longq'; cluster.time=12*60;
+elseif strcmpi(clustername,'larsen'),
+	cluster.name='larsen';   cluster.np=8;    cluster.queue='long';  cluster.time=12*60;
+else
+	cluster.name=oshostname; cluster.np=16;   cluster.queue='long';  cluster.time=12*60;
+end%}}}
+
+%some timing
+t1=clock;
+
+%Mesh generation{{{1
+num=1;
+if ismember(num,steps),
+	message(rs,num)
+
+	%Mesh parameters
+	domain =['DomainOutline.exp'];
+	hmin=3000;
+	hmax=15*10^3;
+	gradation=2;
+	err=[4 40];
+
+	%Initil mesh creation
+	disp('   Interpolating fields'); 
+	md=model;
+	md=bamg(md,'domain',domain,'hmax',hmin,'splitcorner',0,'KeepVertices',0);
+
+	%Mesh adaptation
+	disp('   Interpolating fields'); 
+	vxpath   =[modeldatapath '/surfvelx.mat'];
+	vypath   =[modeldatapath '/surfvely.mat'];
+	md.vx_obs=InterpFromFile(md.x,md.y,vxpath,0);
+	md.vy_obs=InterpFromFile(md.x,md.y,vypath,0);
+	md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
+	thicknesspath=[modeldatapath '/thk.mat'];
+	md.thickness=InterpFromFile(md.x,md.y,thicknesspath,0);
+
+	%Complete with balanced velocities
+	surfacepath  =[modeldatapath '/usrf.mat'];
+	balancedvelpath=[modeldatapath '/balvelmag.mat'];
+	md.surface=InterpFromFile(md.x,md.y,surfacepath,0);
+	[md.surface_slopex md.surface_slopey]=slope(md);
+	md.surface_slopex=averaging(md,md.surface_slopex,0);
+	md.surface_slopey=averaging(md,md.surface_slopey,0);
+	balancedvelocities=InterpFromFile(md.x,md.y,balancedvelpath,0);
+	velx=balancedvelocities.*md.surface_slopex./sqrt(md.surface_slopex.^2+md.surface_slopey.^2);
+	vely=balancedvelocities.*md.surface_slopey./sqrt(md.surface_slopex.^2+md.surface_slopey.^2);
+	pos=find(md.vel_obs==0);
+	md.vx_obs(pos)=velx(pos);
+	md.vy_obs(pos)=vely(pos);
+	md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
+
+	disp('   Mesh adaptation'); 
+	md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',[md.vel_obs md.thickness],'err',[4 40],'splitcorner',0,'KeepVertices',0);
+	
+	savemodel(rs,num,md);
+end %}}}
+%Parameterization{{{1
+num=num+1;
+if ismember(num,steps),
+	message(rs,num)
+	md=loadmodel(rs,'Mesh generation');
+
+	%{
+	landcoverfile=[modeldatapath '/landcover.mat'];
+	load(landcoverfile);
+	land_type=InterpFromGridToMesh(x1,y1,landcover,md.x,md.y,0);
+	gridonicesheet=zeros(md.numberofgrids,1);
+	gridonicesheet(find(land_type>2))=1;
+	elementonicesheet=zeros(md.numberofelements,1);
+	elementonicesheet(find(sum(gridonicesheet(md.elements(:,:)),2)==3))=1;
+	md=geography(md,double(~elementonicesheet),'');
+	%}
+	md=geography(md,'','');
+	md=parameterize(md,'Greenland.par');
+	md=setelementstype(md,'macayeal','all');
+	md=setcluster(md,cluster);
+	
+	savemodel(rs,num,md);
+end%}}}
+%Control Method drag 1{{{1
+num=num+1;
+if ismember(num,steps),
+	message(rs,num)
+	md=loadmodel(rs,'Parameterization');
+
+	md.eps_res=0.01;
+	md.eps_rel=0.1;
+	md.eps_abs=NaN;
+	md.eps_cm=NaN;
+	md.verbose=0;
+
+	%50  (2-0) cmmax=120: misfit(md)  = 8.6417
+	%100 (2-3) cmmax=120: misfit(md)  = 7.4280
+	%100 (2-3) cmmax=200: misfit(md)  = 7.4146
+	%100 (4-3) cmmax=200: misfit(md)  = 8.6725
+	%100 (2-0) cmmax=200: misfit(md)  = 8.4770
+	%100 (2-3) cmmax=200: misfit(md)  = 7.3538
+	md=parametercontroldrag(md,'nsteps',100);
+	%fit=[0;0;2;2];
+	%md.fit=repmat(fit,md.nsteps,1);
+	md.fit=[];
+	md.fit(1:floor(md.nsteps/2),1)=2;
+	md.fit(floor(md.nsteps/2):md.nsteps,1)=3;
+	md.optscal(1:floor(md.nsteps/2))=500;
+	md.optscal(floor(md.nsteps/2):md.nsteps)=100;
+	md.cm_jump=.99*ones(md.nsteps,1);
+	md.control_type='drag';
+	md.waitonlock=inf;
+	md.cm_min=1;
+	md.cm_max=200;
+	md.cm_noisedmp=10^-15; %13 is too much
+
+	%spc all 0 velocities.
+	pos=find(md.vel_obs==0);
+	npos=length(pos);
+	md.spcvelocity(pos,:)=[ones(npos,1) ones(npos,1) ones(npos,1) zeros(npos,1) zeros(npos,1) zeros(npos,1)];
+	in=ContourToMesh(md.elements,md.x,md.y,expread('trackscontours.exp',1),'node',1);
+	md.weights(find(in))=0.1;
+
+	md=setcluster(md,cluster);
+	md.control_analysis=1;
+	md=solve(md,'analysis_type','diagnostic');
+	md.drag=md.results.diagnostic.parameter;
+
+	savemodel(rs,num,md);
+end%}}}
+
+%timing
+t2=clock; disp(sprintf('\n   %s\n',['run elapsed time: ' num2str(etime(t2,t1)) ' seconds']));