Index: /issm/workshop/17_SeaRISE/runme.m
===================================================================
--- /issm/workshop/17_SeaRISE/runme.m	(revision 11052)
+++ /issm/workshop/17_SeaRISE/runme.m	(revision 11053)
@@ -1,11 +1,10 @@
-steps=[1];
+steps=[1:5];
 
 %Hard coded parameters
 cluster=generic('name',oshostname,'np',2);
-%What directory is you nc file in?
-modeldatapath=['./'];
-
+modeldatapath=['/Users/issm/Desktop/issm/workshop/14_JKS/Greenland5km_v1.2'];
 
 %Run Steps
+%Mesh{{{1
 if sum(find(steps==1))
 
@@ -22,27 +21,18 @@
 
 	ncdata=[modeldatapath '/Greenland_5km_dev1.2.nc'];
+	%Use ncdisp('file') to see an ncdump
+	x1=double(ncread(ncdata,'x1'));
+	y1=double(ncread(ncdata,'y1'));
 
-	%Get x1 and y1 from ncdata
-	%Use ncdisp('file') to see an ncdump
+	velx=double(ncread(ncdata,'surfvelx'));
+	vely=double(ncread(ncdata,'surfvely'));
 
+	md.inversion.vx_obs=InterpFromGridToMesh(x1,y1,velx',md.mesh.x,md.mesh.y,0);
+	md.inversion.vy_obs=InterpFromGridToMesh(x1,y1,vely',md.mesh.x,md.mesh.y,0);
+	md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
 
-
-
-	%Get velx and vely from ncdata
-
-
-
-
-
-	%Interpolate grid to mesh velx and vely into md.inversion vx_obs, vy_obs, and
-	%calculate vel_obs
-
-
-
-
-	%Read thickness and Interpolate into md.geometry.thickness
-
-
-
+	thickness=double(ncread(ncdata,'thk'));
+	md.geometry.thickness=InterpFromGridToMesh(x1,y1,thickness',...
+		md.mesh.x,md.mesh.y,0);
 
 	disp('   Mesh adaptation'); 
@@ -54,8 +44,7 @@
 	name='./Models/Greenland.Mesh_generation';
 	save(name,'md');
-end
+end% }}}
+%Parameterization{{{1
 if sum(find(steps==2)) 
-
-	%Set up your Parameter file.  Start by copying your Jks.par
 
 	md=loadmodel('./Models/Greenland.Mesh_generation');
@@ -66,5 +55,6 @@
 	name='./Models/Greenland.Parameterization'; 
 	save(name,'md');
-end
+end%}}}
+%Control_drag{{{1
 if sum(find(steps==3)) 
 
@@ -79,12 +69,13 @@
 	md.inversion.gradient_scaling(:)=100;
 
-
-	%Set up inversion control parameters
-
-
-
-
-
-
+	md.inversion.cost_functions=[];
+	md.inversion.cost_functions(1:floor(md.inversion.nsteps),1)=103;
+	md.inversion.cost_functions=[103*ones(md.inversion.nsteps,1)  ...
+		501*ones(md.inversion.nsteps,1)];
+	md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,2);
+	md.inversion.cost_functions_coefficients(:,2)=8*10^-15;
+	md.inversion.control_parameters={'FrictionCoefficient'};
+	md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
+	md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);
 
 	md=setcluster(md,cluster);
@@ -93,13 +84,12 @@
 	md=solve(md,DiagnosticSolutionEnum);
 
-	%Save results in md.friction.coefficient, md.initialization.vx and md.initialization.vy
-	%before save
-
-
-
+	md.friction.coefficient=md.results.DiagnosticSolution.FrictionCoefficient;
+	md.initialization.vx=md.results.DiagnosticSolution.Vx;
+	md.initialization.vy=md.results.DiagnosticSolution.Vy;
 
 	name='./Models/Greenland.Control_drag'; 
 	save(name,'md');
-end
+end%}}}
+%Transient{{{1
 if sum(find(steps==4))
 
@@ -113,20 +103,15 @@
 	y1=double(ncread(ncdata,'y1'));
 
-
-	%Get the surface mass balance from your nc file (Hint, your par file
-	%can help).  Interpolate and then Make sure you convert to m ice/year.
-
-
-	
-	%Set up forcing (forcing time in final row)
-
-
+	smb=double(ncread(ncdata,'smb'));
+	smb=InterpFromGridToMesh(x1,y1,smb',md.mesh.x,md.mesh.y,0)*...
+		1000/md.materials.rho_ice;
+	smb=[smb smb smb-1.0];
+	md.surfaceforcings.mass_balance=[smb; [1 10 20]];
 
 	md.thermal.spctemperature=[md.initialization.temperature;1]; 
 
-	%Step of 0.2, for 20 years, output every year
-
-
-
+	md.timestepping.time_step=0.2;
+	md.timestepping.final_time=20;
+	md.settings.output_frequency=1;
 
 	md.inversion.iscontrol=0;
@@ -140,5 +125,6 @@
 	name='./Models/Greenland.Transient'; 
 	save(name,'md');
-end% 
+end% }}}
+%Plot{{{1
 if sum(find(steps==5))
 
@@ -165,10 +151,10 @@
 
 	%Plot velocity
-
-
+	vel=[]; for i=1:100; vel=[vel md.results.TransientSolution(i).Vel]; end
+	subplot(3,1,2); plot([0.2:0.2:20],mean(vel)); title('Mean Velocity');
 
 	%Plot Volume
-
-
-
-end% 
+	volume=[]; for i=1:100; volume=[volume md.results.TransientSolution(i).IceVolume]; end
+	subplot(3,1,3); plot([0.2:0.2:20],volume); title('Ice Volume');
+	xlabel('years')
+end% }}}