Index: /issm/workshop/17_SeaRISE/runme.m
===================================================================
--- /issm/workshop/17_SeaRISE/runme.m	(revision 11048)
+++ /issm/workshop/17_SeaRISE/runme.m	(revision 11048)
@@ -0,0 +1,174 @@
+steps=[1];
+
+%Hard coded parameters
+cluster=generic('name',oshostname,'np',2);
+%What directory is you nc file in?
+modeldatapath=['./'];
+
+
+%Run Steps
+if sum(find(steps==1))
+
+	%Mesh parameters
+	domain =['./Exp_Par/DomainOutline.exp'];
+	hmax=400000;
+	hmin=80000;
+	gradation=1.7;
+	err=[15 200];
+
+	%Initial mesh creation
+	disp('   Interpolating fields'); 
+	md=bamg(model,'domain',domain,'hmax',hmin,'splitcorner',1,'KeepVertices',0);
+
+	ncdata=[modeldatapath '/Greenland_5km_dev1.2.nc'];
+
+	%Get x1 and y1 from ncdata
+	%Use ncdisp('file') to see an ncdump
+
+
+
+
+	%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
+
+
+
+
+	disp('   Mesh adaptation'); 
+	md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',...
+		[md.inversion.vel_obs md.geometry.thickness],'err',err,...
+		'splitcorner',1,'KeepVertices',0);
+	[md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,+1,39,71);
+
+	name='./Models/Greenland.Mesh_generation';
+	save(name,'md');
+end
+if sum(find(steps==2)) 
+
+	%Set up your Parameter file.  Start by copying your Jks.par
+
+	md=loadmodel('./Models/Greenland.Mesh_generation');
+
+	md=setmask(md,'','');
+	md=parameterize(md,'./Exp_Par/Greenland.par');
+
+	name='./Models/Greenland.Parameterization'; 
+	save(name,'md');
+end
+if sum(find(steps==3)) 
+
+	md=loadmodel('./Models/Greenland.Parameterization');
+	md=setflowequation(md,'macayeal','all');
+
+	md.diagnostic.restol=0.01;
+	md.diagnostic.reltol=0.1;
+
+	md=parametercontroldrag(md,'nsteps',50);
+	md.inversion.step_threshold=.99*ones(md.inversion.nsteps,1);
+	md.inversion.gradient_scaling(:)=100;
+
+
+	%Set up inversion control parameters
+
+
+
+
+
+
+
+	md=setcluster(md,cluster);
+	md.inversion.iscontrol=1;
+	md.verbose=verbose('solution',true,'control',true);
+	md=solve(md,DiagnosticSolutionEnum);
+
+	%Save results in md.friction.coefficient, md.initialization.vx and md.initialization.vy
+	%before save
+
+
+
+
+	name='./Models/Greenland.Control_drag'; 
+	save(name,'md');
+end
+if sum(find(steps==4))
+
+	md=loadmodel('./Models/Greenland.Control_drag');
+
+	md=setflowequation(md,'macayeal','all');
+	md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1);
+
+	ncdata=[modeldatapath '/Greenland_5km_dev1.2.nc'];
+	x1=double(ncread(ncdata,'x1'));
+	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)
+
+
+
+	md.thermal.spctemperature=[md.initialization.temperature;1]; 
+
+	%Step of 0.2, for 20 years, output every year
+
+
+
+
+	md.inversion.iscontrol=0;
+	md.transient.requested_outputs=IceVolumeEnum();
+
+	md=setcluster(md,cluster);
+	md.verbose=verbose('solution',true,'module',true,'convergence',true);
+
+	md=solve(md,TransientSolutionEnum);
+
+	name='./Models/Greenland.Transient'; 
+	save(name,'md');
+end% 
+if sum(find(steps==5))
+
+	md=loadmodel('./Models/Greenland.Transient');
+
+	%Planview plots
+	plotmodel(md,'data',md.results.TransientSolution(end).Vel,'caxis',[1e-1 5000],...
+		'log', 10, 'title','Velocity (m/y)', ...
+		'data', md.results.TransientSolution(1).SurfaceforcingsMassBalance, ...
+		'title', 'Surface mass balance (m/y)', ...
+		'data',md.results.TransientSolution(end).Thickness,...
+		'title','Thickness (m)', ...
+		'data',md.results.TransientSolution(end).Surface, ...
+		'title', 'Surface (m)');
+
+	%Line Plots
+
+	figure
+
+	%Plot surface mass balance
+	surfmb=[]; for i=1:100; surfmb=[surfmb ...
+		md.results.TransientSolution(i).SurfaceforcingsMassBalance]; end
+	subplot(3,1,1); plot([0.2:0.2:20],mean(surfmb)); title('Mean Surface mass balance');
+
+	%Plot velocity
+
+
+
+	%Plot Volume
+
+
+
+end%