Index: /issm/trunk/scripts/runme.m
===================================================================
--- /issm/trunk/scripts/runme.m	(revision 4985)
+++ /issm/trunk/scripts/runme.m	(revision 4986)
@@ -1,6 +1,6 @@
-%AGU 09: second attempt (after Moca09) at  control methods on entire Antarctica.
+%SaeRISE Antarctica
 
 %Run steps
-steps=[3];
+steps=[1 2 3];
 
 %steps available: 
@@ -13,5 +13,5 @@
 modeldatapath='./../../Data/';
 %clustername=oshostname();
-clustername='cosmos';
+clustername='wilkes';
 type='2d';
 
@@ -24,5 +24,6 @@
 	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;
+	cluster.name=oshostname; cluster.np=8;   cluster.queue='long';  cluster.time=12*60;
+	%cluster.name=oshostname; cluster.np=16;   cluster.queue='long';  cluster.time=12*60;
 end%}}}
 
@@ -30,52 +31,19 @@
 t1=clock;
 
-%Mesh generation{{{1
+%Mesh generation {{{1
 num=1;
 if ismember(num,steps),
 	message(rs,num)
+	
+	%Mesh parameters
+	resolution=15000;
 
-	%Mesh parameters
-	domain =['DomainOutline.exp'];
-	hmin=3000;
-	hmax=15*10^3;
-	gradation=2;
-	err=[4 40];
+	%Initial mesh
+	md=model;
+	md=mesh(md,'DomainOutline.exp',resolution);
 
-	%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
+%Parameterization {{{1
 num=num+1;
 if ismember(num,steps),
@@ -83,22 +51,20 @@
 	md=loadmodel(rs,'Mesh generation');
 
-	%{
-	landcoverfile=[modeldatapath '/landcover.mat'];
-	load(landcoverfile);
-	land_type=InterpFromGridToMesh(x1,y1,landcover,md.x,md.y,0);
+	%Find ice sheets and ice shelves
+	landcoveringpath=[modeldatapath '/thkmask.mat'];
+	land_type=InterpFromFile(md.x,md.y,landcoveringpath,0);
 	gridonicesheet=zeros(md.numberofgrids,1);
-	gridonicesheet(find(land_type>2))=1;
+	gridonicesheet(find(land_type>0))=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');
+	elementoniceshelf=1-elementonicesheet;
+	md=geography(md,elementoniceshelf,'');
+
+	md=parameterize(md,'Antarctica.par');
 	md=setelementstype(md,'macayeal','all');
-	md=setcluster(md,cluster);
 	
 	savemodel(rs,num,md);
-end%}}}
-%Control Method drag 1{{{1
+end %}}}
+%Control method drag 1 {{{
 num=num+1;
 if ismember(num,steps),
@@ -106,24 +72,20 @@
 	md=loadmodel(rs,'Parameterization');
 
-	md.eps_res=0.01;
-	md.eps_rel=0.1;
+	disp('Compute 2d drag coefficient');
+
 	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.drag(:)=30;
+	md.control_analysis=1;
+	md.nsteps=10;
+	md.maxiter=10*ones(md.nsteps,1);
+	fit=[0;0;2;0];
+	md.fit=repmat(fit,md.nsteps,1);
+	md.fit(1:30)=2;
+	md.fit(31:60)=0;
+	md.fit(md.nsteps+1:end)=[];
+	md.optscal=30*ones(md.nsteps,1);
 	md.cm_jump=.99*ones(md.nsteps,1);
 	md.control_type='drag';
@@ -131,5 +93,4 @@
 	md.cm_min=1;
 	md.cm_max=200;
-	md.cm_noisedmp=10^-15; %13 is too much
 
 	%spc all 0 velocities.
@@ -137,14 +98,11 @@
 	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%}}}
+end %}}}
 
 %timing