%Organizers section:  {{{
%org=organizer('repository','./Models','prefix','LarsenC.final.','steps',steps); %using GSHHG meshing.
org=organizer('repository','./Models','prefix','LarsenC.final.highres.','steps',steps); %using GSHHG meshing.
addpath Functions Code

clustername=oshostname(); 
%clustername='pfe';
printresults=1;
loadresults=0;
ncores=8;
interactive=1;
%}}}
%Cluster parameters{{{
pfe=pfe;
pfe.login='elarour';
%pfe.numnodes=20;
%pfe.cpuspernode=20;
pfe.numnodes=20;
pfe.cpuspernode=10;
pfe.time=120*60;
pfe.interactive=interactive; 
pfe.port=1025;
pfe.processor='ivy';
pfe.codepath='/nobackupnfs2/elarour/trunk-larour-NatGeoScience2016/build-fw/bin';
pfe.executionpath='/nobackupp2/elarour/SlrIceCoupled2017/Execution';

cluster_generic=generic('name',oshostname,'np',12,'verbose',0);
cluster_generic.codepath=[issmdir '/build-fw/bin'];

if strcmpi(clustername,'pfe'),
	cluster=pfe;
elseif strcmpi(clustername,'local_pfe'),
	cluster=local_pfe_fw;
elseif strcmpi(clustername,oshostname()),
	cluster=cluster_generic;
else
	error('cluster not supported yet');
end% }}}
%Data paths %{{{
fontsize=18;
linewidth=2;
modeldatapath='/Users/larour/ModelData';
shppath='/Users/larour/issm-jpl/proj-group/qgis/LarsenC';

gracetimeseries='/Users/larour/issm-jpl/proj-group/GRACE/Data/GRACE_SS41z_n719_p25deg_remOC_GRN_GIAg.mat';

%antactica
firnpath  ='/Users/larour/ModelData/SeaRISE/Antarctica5km_shelves_v1.0/firn.mat';

%greenland: 
MARname=[modeldatapath '/MAR3.5/'];
searise_gpath           =[modeldatapath '/SeaRISE/Greenland5km_v1.2'];
seariseg_vxpath         =[searise_gpath '/surfvelx.mat'];
seariseg_vypath         =[searise_gpath '/surfvely.mat'];
seariseg_temperature    =[searise_gpath '/surftemp.mat'];
seariseg_heatflux       =[searise_gpath '/bheatflx.mat'];
rignotg_velpath=[modeldatapath '/RignotGreenland2012Vel/RignotGreenland2012Vel.mat'];
bamberg_velpath=[modeldatapath '/BamberGreenlandVel2km/velocity.mat'];
boxg_smb=[modeldatapath '/BoxGreenlandSMB2011/smb150.mat'];
boxg_smbmean=[modeldatapath '/BoxGreenlandSMB2014/meansmb170_int.mat'];
gshhsshapefile=[modeldatapath '/Gshhg/Shp/GSHHS_shp/c/GSHHS_c_L1.shp'];

%global
etopo=[modeldatapath '/Etopo2v2/ETOPO2v2c_f4.nc'];
openstreetgreenland=[modeldatapath '/OpenStreet/coastlines-split-4326/coastlines_z5_3413.shp'];
openstreetantarctica=[modeldatapath '/OpenStreet/coastlines-split-3031/coastline_z5-antarctica-3031.shp'];
rgigreenland=[modeldatapath '/Rgi/05_rgi60_GreenlandPeriphery/05_rgi60_GreenlandPeriphery-3413.shp'];
rgiantarctica=[modeldatapath '/Rgi/19_rgi60_AntarcticSubantarctic/19_rgi60_AntarcticSubantarctic-3031.shp'];
holesshpgreenland='/Users/larour/issm-jpl/proj-group/qgis/Greenland/RGI_1km.shp';
holesshpgreenland5km='/Users/larour/issm-jpl/proj-group/qgis/Greenland/RGI_5km_DownSelect5km.shp';
holesshpantarctica='/Users/larour/issm-jpl/proj-group/qgis/Antarctica/RGI_1km_OutsideAntNoHolesSimp1km.shp';

%velocities and surface: 
%velpath=[shppath '/ANNUAL/2010-07-01_2011-06-31_LC/vel_final_2010-07-01_2011-06-31_LC.nc'];
%velpath=[shppath '/ANNUAL/2016-07-01_2017-06-31_LC/vel_final_2016-07-01_2017-06-31_LC.nc'];
%velpath=[shppath '/ANNUAL/2013-07-01_2014-06-31_LC/vel_final_2013-07-01_2014-06-31_LC.nc'];
velpath=[shppath '/ANNUAL/2014-07-01_2015-06-31_LC/vel_final_2014-07-01_2015-06-31_LC.nc'];
surfacepath='/Users/larour/issm-jpl/proj-larour/LarsenC-2020/Data/TDX_ANTPEN450m_msl_v18july2017.tif';
%rift='./Exp/Rift2014.shp';
riftpassive=[shppath '/Rift2014_Passive.shp'];
rift=[shppath '/Rift2014_Active.shp'];
%rift=[shppath '/Rift.shp'];
meltratepath='./Data/LarsenC_CS2_2013-ERS1_1994_thickness_difference.tif';

%domain='./Exp/IceShelf.shp';
%domain='./Exp/IceShelf2010.shp';
domain='./Exp/IceShelf2013.shp';
%domain=[shppath '/DomainOutlineMerged.shp'];

%}}}
if perform(org,'Mesh'), % {{{

	%mesh parameters:  
	hmax=20000; hmin=1500; err=[1.5 .05]; gradation=1.5; loop=3; 
	%hmax=50000; hmin=2500; err=[5.5 1.5]; gradation=1.5; loop=2; 

	
	md=bamg(model,'domain',domain,'rifts',rift,'hmin',hmin,'KeepVertices',0,'MaxCornerAngle',0.0000000001);
	
	%loop on refinement: 
	for i=1:loop,

		%interpolate  vel and thickness{{{
		x=ncread(velpath,'x'); y=ncread(velpath,'y'); vx=ncread(velpath,'VX'); vy=ncread(velpath,'VY');
		velgrid=sqrt(vx.^2+vy.^2);
		vel=InterpFromGridToMesh(x,flipud(y),flipud(velgrid'),md.mesh.x,md.mesh.y,NaN);
		 
		pos=find(isnan(vel) | vel>1700 );  pos2=find(vel<=1700 & ~isnan(vel));
		vel(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),vel(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');
		
		info=gdalinfo(surfacepath);
		x1=info.xmin:info.dx:info.xmax; x1=x1';
		y1=info.ymin:info.dy:info.ymax; y1=y1';
		demread=imread(surfacepath);
		demread=flipud(demread);
		surface   = InterpFromGridToMesh(x1,y1,demread,md.mesh.x,md.mesh.y,NaN);

		%interpolate NaN values away: 
		surface=griddata(md.mesh.x(~isnan(surface)),md.mesh.y(~isnan(surface)),surface(~isnan(surface)),md.mesh.x,md.mesh.y,'nearest');
		%}}}
		%refine tips {{{
		tipres=50;
		tips=rifttipsonmesh(md,rift);
		hVertices=NaN*ones(md.mesh.numberofvertices,1);
		hVertices(tips)=tipres/10;
		%}}}
		
		md=bamg(md,'hVertices',hVertices,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',[vel surface],'err',err,'splitcorner',1,'MaxCornerAngle',1e-15,'KeepVertices',1);
	end

	md=meshprocessrifts(md,rift);


	%fix holes: {{{
	fix=1;
	if fix,
		for r=1:length(md.rifts.riftstruct),

			pos=md.rifts.riftstruct(r).tips; 
			for i=3:4,
				tip=pos(i);  rep=pos(i-2);
				pos2=find(md.mesh.elements==tip); md.mesh.elements(pos2)=rep;
			end
			md.rifts.riftstruct(r).tips= md.rifts.riftstruct(r).tips(1:2);
		end
		n=length(md.rifts.riftstruct);

		md.mesh.x=md.mesh.x(1:end-2*n);
		md.mesh.y=md.mesh.y(1:end-2*n);
		md.mesh.numberofvertices= md.mesh.numberofvertices-2*n;
		md.mesh.vertexonboundary=md.mesh.vertexonboundary(1:end-2*n);
	end % }}}

	%interpolate  vel and thickness{{{
	x=ncread(velpath,'x'); y=ncread(velpath,'y'); vx=ncread(velpath,'VX'); vy=ncread(velpath,'VY');
	velgrid=sqrt(vx.^2+vy.^2);
	vel=InterpFromGridToMesh(x,flipud(y),flipud(velgrid'),md.mesh.x,md.mesh.y,NaN);

	pos=find(isnan(vel) | vel>1700 );  pos2=find(vel<=1700 & ~isnan(vel));
	vel(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),vel(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');

	info=gdalinfo(surfacepath);
	x1=info.xmin:info.dx:info.xmax; x1=x1';
	y1=info.ymin:info.dy:info.ymax; y1=y1';
	demread=imread(surfacepath);
	demread=flipud(demread);
	surface   = InterpFromGridToMesh(x1,y1,demread,md.mesh.x,md.mesh.y,NaN);

	%interpolate NaN values away: 
	surface=griddata(md.mesh.x(~isnan(surface)),md.mesh.y(~isnan(surface)),surface(~isnan(surface)),md.mesh.x,md.mesh.y,'nearest');
	%}}}

	md.geometry.surface=surface;
	md.inversion.vel_obs=vel;
	savemodel(org,md);

end % }}}
if perform(org,'Velocities'), % {{{

	md=loadmodel(org,'Mesh');

	vel2000=[shppath '/2000_test_vel/Antarctica_ice_velocity_2000_450m_v1.nc'];
	vel2010=[shppath '/ANNUAL/2010-07-01_2011-06-31_LC/vel_final_2010-07-01_2011-06-31_LC.nc'];
	vel2013=[shppath '/ANNUAL/2013-07-01_2014-06-31_LC/vel_final_2013-07-01_2014-06-31_LC.nc'];
	vel2014=[shppath '/ANNUAL/2014-07-01_2015-06-31_LC/vel_final_2014-07-01_2015-06-31_LC.nc'];
	vel2017=[shppath '/ANNUAL/2016-07-01_2017-06-31_LC/vel_final_2016-07-01_2017-06-31_LC.nc'];

	vel=vel2000;
	x=ncread(vel,'x'); y=ncread(vel,'y'); vx=ncread(vel,'VX'); vy=ncread(vel,'VY');
	velgrid=sqrt(vx.^2+vy.^2);
	vel=InterpFromGridToMesh(x,flipud(y),flipud(velgrid'),md.mesh.x,md.mesh.y,NaN);

	pos=find(isnan(vel) | vel>1700 );  pos2=find(vel<=1700 & ~isnan(vel));
	vel(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),vel(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');
	vel2000=vel;
	
	vel=vel2010;
	x=ncread(vel,'x'); y=ncread(vel,'y'); vx=ncread(vel,'VX'); vy=ncread(vel,'VY');
	velgrid=sqrt(vx.^2+vy.^2);
	vel=InterpFromGridToMesh(x,flipud(y),flipud(velgrid'),md.mesh.x,md.mesh.y,NaN);

	pos=find(isnan(vel) | vel>1700 );  pos2=find(vel<=1700 & ~isnan(vel));
	vel(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),vel(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');
	vel2010=vel;
	
	
	vel=vel2013;
	x=ncread(vel,'x'); y=ncread(vel,'y'); vx=ncread(vel,'VX'); vy=ncread(vel,'VY');
	velgrid=sqrt(vx.^2+vy.^2);
	vel=InterpFromGridToMesh(x,flipud(y),flipud(velgrid'),md.mesh.x,md.mesh.y,NaN);

	pos=find(isnan(vel) | vel>1700 );  pos2=find(vel<=1700 & ~isnan(vel));
	vel(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),vel(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');
	vel2013=vel;
	
	
	vel=vel2014;
	x=ncread(vel,'x'); y=ncread(vel,'y'); vx=ncread(vel,'VX'); vy=ncread(vel,'VY');
	velgrid=sqrt(vx.^2+vy.^2);
	vel=InterpFromGridToMesh(x,flipud(y),flipud(velgrid'),md.mesh.x,md.mesh.y,NaN);

	pos=find(isnan(vel) | vel>1700 );  pos2=find(vel<=1700 & ~isnan(vel));
	vel(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),vel(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');
	vel2014=vel;



	vel=vel2017;
	x=ncread(vel,'x'); y=ncread(vel,'y'); vx=ncread(vel,'VX'); vy=ncread(vel,'VY');
	velgrid=sqrt(vx.^2+vy.^2);
	vel=InterpFromGridToMesh(x,flipud(y),flipud(velgrid'),md.mesh.x,md.mesh.y,NaN);

	pos=find(isnan(vel) | vel>1700 );  pos2=find(vel<=1700 & ~isnan(vel));
	vel(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),vel(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');
	vel2017=vel;

	plotmodel(md,'data',vel2000,'title','2000',...
	'data',vel2010,'title','2010',...
	'data',vel2013,'title','2013',...
	'data',vel2014,'title','2014',...
	'data',vel2017,'title','2017',...
	'caxis#all',[0 1000],'nlines',3,'ncols',2);


	savemodel(org,md);

end % }}}
if perform(org,'Parameterization'),% {{{

	md=loadmodel(org,'Mesh');
	continent='Antarctica';

	%masks and grounding line:  %{{{
	md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1);
	flags=ContourToNodes(md.mesh.x,md.mesh.y,[shppath '/IceFront.shp'],0);
	pos=find(flags & md.mesh.vertexonboundary); 
	md.mask.ice_levelset(pos)=0;
	
	%figure out initial grounding line position: 
	land_type=interpBedmap2(md.mesh.x,md.mesh.y,'icemask_grounded_and_shelves');
	
	pos=find(isnan(land_type)); 
	pos2=find(~isnan(land_type)); 
	for i=1:length(pos),
		ind=find_point(md.mesh.x(pos2),md.mesh.y(pos2),md.mesh.x(pos(i)),md.mesh.y(pos(i)));
		land_type(pos(i))=land_type(pos2(ind));
	end

	gridoniceshelf=zeros(md.mesh.numberofvertices,1);
	gridoniceshelf(land_type>0.9)=land_type(land_type>0.9);
	gridoniceshelf(gridoniceshelf>0)=-1;
	gridoniceshelf(gridoniceshelf~=-1)=1;

	md.mask.groundedice_levelset=gridoniceshelf;
	%}}}
	%latlong:  % {{{
	[md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,'EPSG:3031','EPSG:4326'); 
	md.mesh.epsg=3031;
	%}}}
	%geometry: {{{
	di=md.materials.rho_ice/md.materials.rho_water;

	disp('      reading TSX dem');
	dem='./Data/TDX_ANTPEN450m_msl_v18july2017.tif';
	info=gdalinfo(dem);
	x1=info.xmin:info.dx:info.xmax; x1=x1';
	y1=info.ymin:info.dy:info.ymax; y1=y1';

	%Interpolate:
	demread=imread(dem);
	demread=flipud(demread);
	%pos=find(demread<-80000); demread(pos)=NaN;
	surface   = InterpFromGridToMesh(x1,y1,demread,md.mesh.x,md.mesh.y,NaN);

	%interpolate NaN values away: 
	surface=griddata(md.mesh.x(~isnan(surface)),md.mesh.y(~isnan(surface)),surface(~isnan(surface)),md.mesh.x,md.mesh.y,'nearest');
	
	%thickness:
	thickness=1/(1-di)*surface;

	disp('      reading firn layer');
	load(firnpath)
	firn_layer=InterpFromGridToMesh(x1,y1,firn,md.mesh.x,md.mesh.y,0);
	thickness=thickness-firn_layer;
	pos=find(thickness<10); thickness(pos)=10;

	%surface and base: 
	surface=(1-di)*thickness; 
	base=-di*thickness;
	bed=base-100;

	md.geometry.thickness=thickness;
	md.geometry.surface=surface;
	md.geometry.base=base;
	md.geometry.bed=bed;
	% }}}
	%velocities:  %{{{
	%interpolate  velpath {{{
	x=ncread(velpath,'x'); y=ncread(velpath,'y'); vxgrid=ncread(velpath,'VX'); vygrid=ncread(velpath,'VY');
	velgrid=sqrt(vxgrid.^2+vygrid.^2);
	vel=InterpFromGridToMesh(x,flipud(y),flipud(velgrid'),md.mesh.x,md.mesh.y,NaN);
	vx=InterpFromGridToMesh(x,flipud(y),flipud(vxgrid'),md.mesh.x,md.mesh.y,NaN);
	vy=InterpFromGridToMesh(x,flipud(y),flipud(vygrid'),md.mesh.x,md.mesh.y,NaN);

	pos=find(isnan(vel) | vel>1700 );  pos2=find(vel<=1700 & ~isnan(vel));
	vel(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),vel(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');
	vx(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),vx(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');
	vy(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),vy(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest');
	%}}}

	md.inversion.vx_obs=vx;
	md.inversion.vy_obs=vy;
	md.inversion.vel_obs=vel;
	md.initialization.vx=md.inversion.vx_obs;
	md.initialization.vy=md.inversion.vy_obs;
	md.initialization.vz=zeros(md.mesh.numberofvertices,1);
	md.initialization.vel=md.inversion.vel_obs;
	 
	% }}} 
	%friction:  {{{
	[sx,sy,s]=slope(md); sslope=averaging(md,s,0);
	
	pos=find(md.inversion.vel_obs==0);
	vel=max(md.inversion.vel_obs,0.1); vel(pos)=1;
	
	md.friction.coefficient=sqrt(md.materials.rho_ice*md.geometry.thickness.*(sslope)./((md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.base).*vel/md.constants.yts));
	md.friction.coefficient=min(md.friction.coefficient,200);

	min_drag_coeff=35;
	min_drag_coeff_outlets=5;
	threshhold_drag_coeff=70;
	background_drag_coeff=200;

	maxvel=30;
	pos=find(md.friction.coefficient<min_drag_coeff);
	md.friction.coefficient(pos)=min_drag_coeff;

	%Take care of iceshelves: no drag md.drag
	pos=find(md.mask.groundedice_levelset<0);
	md.friction.coefficient(pos)=0;
	md.friction.p=ones(md.mesh.numberofelements,1);
	md.friction.q=ones(md.mesh.numberofelements,1);
	% }}}
	%Temperatures and surface mass balance: {{{ 
	smbpath        =['/Users/larour/ModelData/RACMO2Antarctica/annsmbfill33.mat'];
	temperaturepath=['/Users/larour/ModelData/RACMO2Antarctica/annt33.mat'];

	load(temperaturepath);
	index=delaunay(x1,y1);
	md.initialization.temperature=InterpFromMeshToMesh2d(index,x1(:),y1(:),annt33(:),md.mesh.x,md.mesh.y,'default',273.15);
	md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);

	load(smbpath);
	smb=InterpFromMeshToMesh2d(index,x1(:),y1(:),annsmbfill33(:),md.mesh.x,md.mesh.y,'default',NaN);
	md.smb.mass_balance = [smb; 1];
	%}}}
	%Mechanical boundary conditions: {{{
	%intialize spc boundary conditions:
	md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
	md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
	md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
	md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
	
	%constrain boundaries of basin to observed velocities: 
	pos=find(md.mesh.vertexonboundary & md.mask.ice_levelset<0); 
	md.stressbalance.spcvx(pos)=md.inversion.vx_obs(pos);
	md.stressbalance.spcvy(pos)=md.inversion.vy_obs(pos);
	md.stressbalance.spcvz(pos)=0;

	pos=find(md.mask.groundedice_levelset>0); 
	md.stressbalance.spcvx(pos)=md.inversion.vx_obs(pos);
	md.stressbalance.spcvy(pos)=md.inversion.vy_obs(pos);
	md.stressbalance.spcvz(pos)=0;

	%prognostic
	md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
	% }}}
	%Thermal boundary conditions: {{{
	searisepath  =[modeldatapath '/SeaRISE/Antarctica5km_shelves_v1.0'];
	heatfluxpath   =[searisepath '/bheatflx_fox.mat'];

	%melt rates: {{{
	md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
	md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
	load('./Data/monthly_mean_melt_adj09.mat'); fw2=fw_flux;
	load('./Data/monthly_mean_melt_adj04.mat'); fw1=fw_flux;
	fw_flux(:,:,end+1:end+size(fw2,3))=fw2;
	load('./Data/grid_info_cp09.mat')
	[xi,yi]= ll2xy(lat,lon,-1,0,71); index=delaunay(xi,yi);
	mrate_mesh=InterpFromMeshToMesh2d(index,xi(:),yi(:),-1*reshape(squeeze(mean(fw_flux,3)),prod(size(xi)),1),md.mesh.x,md.mesh.y,'default',0);
	pos=find(md.mask.groundedice_levelset<0 & mrate_mesh~=0);                                                                                   md.basalforcings.floatingice_melting_rate(pos)=mrate_mesh(pos);
	% }}}

	md.thermal.spctemperature=[md.initialization.temperature;1]; %impose observed temperature on surface
	load(heatfluxpath)
	md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,bheatflx_fox,md.mesh.x,md.mesh.y,0);
	%}}}
	%Slr: {{{
	md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
	md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
	md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
	%}}}
	% material properties: {{{
	disp('      creating flow law paramter');
	md.materials.rheology_B=paterson(md.initialization.temperature-4);
	md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
	%}}}
	%diverse: {{{
	md.miscellaneous.name='LarsenC';
	% }}}
	%rifts: %{{{
	for i=1:length(md.rifts.riftstruct),
		md.rifts.riftstruct(i).fill='Melange';
		md.rifts.riftstruct(i).fraction=0;
		md.rifts.riftstruct(i).fractionincrement=.1;
		md.stressbalance.rift_penalty_lock=2;
		md.stressbalance.rift_penalty_threshold=0;
	end
	%}}}
	%flowequation: % {{{
	md=setflowequation(md,'SSA','all');
	pos=find(md.flowequation.vertex_equation==0); 
	md.flowequation.vertex_equation(pos)=2;
	%}}}
	
	savemodel(org,md) ;
end 
% }}}
if perform(org,'MeltRateBamber'),% {{{

	md=loadmodel(org,'Parameterization');
	md.cluster=cluster;
	
	info=gdalinfo(meltratepath);
	x1=info.xmin:info.dx:info.xmax; x1=x1';
	y1=info.ymin:info.dy:info.ymax; y1=y1';
	meltread=imread(meltratepath);
	meltread=flipud(meltread);
	meltrate   = InterpFromGridToMesh(x1,y1,meltread,md.mesh.x,md.mesh.y,NaN);

	meltrate(find(isinf(meltrate)))=NaN;
	meltrate=griddata(md.mesh.x(~isnan(meltrate)),md.mesh.y(~isnan(meltrate)),meltrate(~isnan(meltrate)),md.mesh.x,md.mesh.y,'nearest');

	meltrate=meltrate/(2013-1994+1);

	savedata(org,meltrate);
end% }}}
if perform(org,'StressBalance'),% {{{

	md=loadmodel(org,'Parameterization');
	md.cluster=cluster;
	
	md.materials.rheology_B(:)=1.7*1e8;
	md.settings.solver_residue_threshold=5*1e-6;

	md=solve(md,'sb');

	md.initialization.vx=md.results.StressbalanceSolution.Vx;
	md.initialization.vy=md.results.StressbalanceSolution.Vy;
	md.initialization.vel=md.results.StressbalanceSolution.Vel;

	savemodel(org,md);
end% }}}
if perform(org,'ControlPreThermal'),% {{{

	md=loadmodel(org,'Stressbalance');
	md.verbose=verbose(0);
	md.cluster=cluster;

	md=setBcontrol(md,27,6e-3,.5e-25); 

	md.materials.rheology_B(:)=1.5e8;

	md.miscellaneous.name='ControlPreThermalB';
	md=solve(md,'Stressbalance');
	md.materials.rheology_B=md.results.StressbalanceSolution.MaterialsRheologyBbar;
	
	%Another stress balance without control to get initial velocity right: 
	md.inversion.iscontrol=0;
	md.miscellaneous.name='ControlPreThermalB';
	md=solve(md,'Stressbalance');
	md.initialization.vx=md.results.StressbalanceSolution.Vx;
	md.initialization.vy=md.results.StressbalanceSolution.Vy;
	md.initialization.vel=md.results.StressbalanceSolution.Vel;

	savemodel(org,md);
end% }}}
if perform(org,'ThinningCurve'),% {{{

	md=loadmodel(org,'ControlPreThermal');
	md.cluster=cluster;
	di=md.materials.rho_ice/md.materials.rho_water;
	md.settings.solver_residue_threshold=1e-4;
	
	N=100;

	%percentage reduction in ice rigidity: 
	t_pmp=md.materials.meltingpoint-md.materials.beta*md.initialization.pressure;
	Bi=zeros(md.mesh.numberofvertices,1);
	Bf=zeros(md.mesh.numberofvertices,1);

	thinning=[0:.1:1,1.1:.2:3,3.2:.5:10,11:1:30];
	nmodels=length(thinning);
	models=cell(nmodels,1);
	penetration=zeros(nmodels,1);
	for j=1:nmodels,
		ex=thinning(j);
		if mod(j,5)==0,
			disp(sprintf('%i/%i',j,nmodels));
		end
		deltaH=ex*ones(md.mesh.numberofvertices,1);

		for i=1:md.mesh.numberofvertices,
			h=md.geometry.base(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bi(i)=mean(paterson(t));
			
			h=md.geometry.base(i)+deltaH(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bf(i)=mean(paterson(t));
		end
		scale=1+(Bf-Bi)./Bi;
		scale=griddata(md.mesh.x(~isnan(scale)),md.mesh.y(~isnan(scale)),scale(~isnan(scale)),md.mesh.x,md.mesh.y,'nearest');

		md.geometry.thickness=md.geometry.thickness-deltaH;
		pos=find(md.geometry.thickness<1); md.geometry.thickness(pos)=1;
		md.geometry.surface=(1-di)*md.geometry.thickness;
		md.geometry.base=-di*md.geometry.thickness;

		md.materials.rheology_B=md.materials.rheology_B.*scale;
		
		md=solve(md,'Stressbalance');
		md.initialization.vx=md.results.StressbalanceSolution.Vx;
		md.initialization.vy=md.results.StressbalanceSolution.Vy;
		md.initialization.vel=md.results.StressbalanceSolution.Vel;

		penetration(j)=meanriftdeltavel(md,1);
		models{j}=md;
	end
	
	penetration_std=penetration; 
	savedata(org,penetration_std,models,thinning);
end% }}}
if perform(org,'ThinningCurveMelangeOnly'),% {{{

	md=loadmodel(org,'ControlPreThermal');
	md.cluster=cluster;
	di=md.materials.rho_ice/md.materials.rho_water;
	md.settings.solver_residue_threshold=1e-4;
	
	N=100;

	%percentage reduction in ice rigidity: 
	t_pmp=md.materials.meltingpoint-md.materials.beta*md.initialization.pressure;
	Bi=zeros(md.mesh.numberofvertices,1);
	Bf=zeros(md.mesh.numberofvertices,1);

	
	thinning=[0:.1:1,1.1:.2:3,3.2:.5:10,11:1:30];
	nmodels=length(thinning);
	models=cell(nmodels,1);
	penetration=zeros(nmodels,1);
	for j=1:nmodels,
		ex=thinning(j);
		if mod(j,5)==0,
			disp(sprintf('%i/%i',j,nmodels));
		end

		deltaH=ex*ones(md.mesh.numberofvertices,1);

		%melange only:
		pos=find(md.geometry.thickness>250);  deltaH(pos)=0;

		for i=1:md.mesh.numberofvertices,
			h=md.geometry.base(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bi(i)=mean(paterson(t));
			
			h=md.geometry.base(i)+deltaH(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bf(i)=mean(paterson(t));
		end
		scale=1+(Bf-Bi)./Bi;
		scale=griddata(md.mesh.x(~isnan(scale)),md.mesh.y(~isnan(scale)),scale(~isnan(scale)),md.mesh.x,md.mesh.y,'nearest');

		md.geometry.thickness=md.geometry.thickness-deltaH;
		pos=find(md.geometry.thickness<1); md.geometry.thickness(pos)=1;
		md.geometry.surface=(1-di)*md.geometry.thickness;
		md.geometry.base=-di*md.geometry.thickness;

		md.materials.rheology_B=md.materials.rheology_B.*scale;
		
		md=solve(md,'Stressbalance');
		md.initialization.vx=md.results.StressbalanceSolution.Vx;
		md.initialization.vy=md.results.StressbalanceSolution.Vy;
		md.initialization.vel=md.results.StressbalanceSolution.Vel;

		penetration(j)=meanriftdeltavel(md,1);
		models{j}=md;
	end
	
	penetration_mel=penetration;
	savedata(org,penetration_mel,models,thinning);
end% }}}
if perform(org,'ThinningCurveIceShelfOnly'),% {{{

	md=loadmodel(org,'ControlPreThermal');
	md.cluster=cluster;
	di=md.materials.rho_ice/md.materials.rho_water;
	md.settings.solver_residue_threshold=1e-4;
	
	N=100;

	%percentage reduction in ice rigidity: 
	t_pmp=md.materials.meltingpoint-md.materials.beta*md.initialization.pressure;
	Bi=zeros(md.mesh.numberofvertices,1);
	Bf=zeros(md.mesh.numberofvertices,1);
	
	
	thinning=[0:.1:1,1.1:.2:3,3.2:.5:10,11:1:30];
	nmodels=length(thinning);
	models=cell(nmodels,1);
	penetration=zeros(nmodels,1);
	for j=1:nmodels,
		ex=thinning(j);
		if mod(j,5)==0,
			disp(sprintf('%i/%i',j,nmodels));
		end

		deltaH=ex*ones(md.mesh.numberofvertices,1);

		%ice shelf only:
		pos=find(md.geometry.thickness<250);  deltaH(pos)=0;

		for i=1:md.mesh.numberofvertices,
			h=md.geometry.base(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bi(i)=mean(paterson(t));
			
			h=md.geometry.base(i)+deltaH(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bf(i)=mean(paterson(t));
		end
		scale=1+(Bf-Bi)./Bi;
		scale=griddata(md.mesh.x(~isnan(scale)),md.mesh.y(~isnan(scale)),scale(~isnan(scale)),md.mesh.x,md.mesh.y,'nearest');

		md.geometry.thickness=md.geometry.thickness-deltaH;
		pos=find(md.geometry.thickness<1); md.geometry.thickness(pos)=1;
		md.geometry.surface=(1-di)*md.geometry.thickness;
		md.geometry.base=-di*md.geometry.thickness;

		md.materials.rheology_B=md.materials.rheology_B.*scale;
		
		md=solve(md,'Stressbalance');
		md.initialization.vx=md.results.StressbalanceSolution.Vx;
		md.initialization.vy=md.results.StressbalanceSolution.Vy;
		md.initialization.vel=md.results.StressbalanceSolution.Vel;

		penetration(j)=meanriftdeltavel(md,1);
		models{j}=md;
	end
	
	penetration_ice=penetration;
	savedata(org,penetration_ice,models,thinning);
end% }}}
if perform(org,'ThinningCurveBamber'),% {{{

	md=loadmodel(org,'ControlPreThermal');
	md.cluster=cluster;
	di=md.materials.rho_ice/md.materials.rho_water;
	loaddata(org,'MeltRateBamber');
	md.settings.solver_residue_threshold=1e-4;
	
	N=100;

	%percentage reduction in ice rigidity: 
	t_pmp=md.materials.meltingpoint-md.materials.beta*md.initialization.pressure;
	Bi=zeros(md.mesh.numberofvertices,1);
	Bf=zeros(md.mesh.numberofvertices,1);

	thinning=[0:.1:1,1.1:.2:3,3.2:.5:10,11:1:30];
	nmodels=31;
	models=cell(nmodels,1);
	penetration=zeros(nmodels,1);
	for j=1:nmodels,

		if mod(j,5)==0,
			disp(sprintf('%i/%i',j,nmodels));
		end
	
		deltaH=meltrate*(j-1);

		for i=1:md.mesh.numberofvertices,
			h=md.geometry.base(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bi(i)=mean(paterson(t));
			
			h=md.geometry.base(i)+deltaH(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bf(i)=mean(paterson(t));
		end
		scale=1+(Bf-Bi)./Bi;
		scale=griddata(md.mesh.x(~isnan(scale)),md.mesh.y(~isnan(scale)),scale(~isnan(scale)),md.mesh.x,md.mesh.y,'nearest');

		md.geometry.thickness=md.geometry.thickness-deltaH;
		pos=find(md.geometry.thickness<1); md.geometry.thickness(pos)=1;
		md.geometry.surface=(1-di)*md.geometry.thickness;
		md.geometry.base=-di*md.geometry.thickness;

		md.materials.rheology_B=md.materials.rheology_B.*scale;
		
		md=solve(md,'Stressbalance');
		md.initialization.vx=md.results.StressbalanceSolution.Vx;
		md.initialization.vy=md.results.StressbalanceSolution.Vy;
		md.initialization.vel=md.results.StressbalanceSolution.Vel;

		penetration(j)=meanriftdeltavel(md,1);
		models{j}=md;
	end
	
	penetration_std=penetration; 
	savedata(org,penetration_std,models,thinning);
end% }}}
if perform(org,'ThinningCurveMelangeOnlyBamber'),% {{{

	md=loadmodel(org,'ControlPreThermal');
	loaddata(org,'MeltRateBamber');
	md.cluster=cluster;
	di=md.materials.rho_ice/md.materials.rho_water;
	md.settings.solver_residue_threshold=1e-4;
	
	N=100;

	%percentage reduction in ice rigidity: 
	t_pmp=md.materials.meltingpoint-md.materials.beta*md.initialization.pressure;
	Bi=zeros(md.mesh.numberofvertices,1);
	Bf=zeros(md.mesh.numberofvertices,1);

	
	thinning=[0:.1:1,1.1:.2:3,3.2:.5:10,11:1:30];
	nmodels=31;
	models=cell(nmodels,1);
	penetration=zeros(nmodels,1);
	for j=1:nmodels,
	
		if mod(j,5)==0,
			disp(sprintf('%i/%i',j,nmodels));
		end
	
		deltaH=meltrate*(j-1);
		
		%melange only:
		pos=find(md.geometry.thickness>250);  deltaH(pos)=0;

		for i=1:md.mesh.numberofvertices,
			h=md.geometry.base(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bi(i)=mean(paterson(t));
			
			h=md.geometry.base(i)+deltaH(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bf(i)=mean(paterson(t));
		end
		scale=1+(Bf-Bi)./Bi;
		scale=griddata(md.mesh.x(~isnan(scale)),md.mesh.y(~isnan(scale)),scale(~isnan(scale)),md.mesh.x,md.mesh.y,'nearest');

		md.geometry.thickness=md.geometry.thickness-deltaH;
		pos=find(md.geometry.thickness<1); md.geometry.thickness(pos)=1;
		md.geometry.surface=(1-di)*md.geometry.thickness;
		md.geometry.base=-di*md.geometry.thickness;

		md.materials.rheology_B=md.materials.rheology_B.*scale;
		
		md=solve(md,'Stressbalance');
		md.initialization.vx=md.results.StressbalanceSolution.Vx;
		md.initialization.vy=md.results.StressbalanceSolution.Vy;
		md.initialization.vel=md.results.StressbalanceSolution.Vel;

		penetration(j)=meanriftdeltavel(md,1);
		models{j}=md;
	end
	
	penetration_mel=penetration;
	savedata(org,penetration_mel,models,thinning);
end% }}}
if perform(org,'ThinningCurveIceShelfOnlyBamber'),% {{{

	md=loadmodel(org,'ControlPreThermal');
	md.cluster=cluster;
	di=md.materials.rho_ice/md.materials.rho_water;
	md.settings.solver_residue_threshold=1e-4;
	
	N=100;

	%percentage reduction in ice rigidity: 
	t_pmp=md.materials.meltingpoint-md.materials.beta*md.initialization.pressure;
	Bi=zeros(md.mesh.numberofvertices,1);
	Bf=zeros(md.mesh.numberofvertices,1);
	
	
	thinning=[0:.1:1,1.1:.2:3,3.2:.5:10,11:1:30];
	nmodels=31;
	models=cell(nmodels,1);
	penetration=zeros(nmodels,1);
	for j=1:nmodels,
	
		if mod(j,5)==0,
			disp(sprintf('%i/%i',j,nmodels));
		end

		deltaH=meltrate*(j-1);

		%ice shelf only:
		pos=find(md.geometry.thickness<250);  deltaH(pos)=0;

		for i=1:md.mesh.numberofvertices,
			h=md.geometry.base(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bi(i)=mean(paterson(t));
			
			h=md.geometry.base(i)+deltaH(i):md.geometry.thickness(i)/N:md.geometry.surface(i); 
			t=h; t=t_pmp(i)-(h-md.geometry.base(i))/md.geometry.thickness(i)*(t_pmp(i)-md.initialization.temperature(i));
			Bf(i)=mean(paterson(t));
		end
		scale=1+(Bf-Bi)./Bi;
		scale=griddata(md.mesh.x(~isnan(scale)),md.mesh.y(~isnan(scale)),scale(~isnan(scale)),md.mesh.x,md.mesh.y,'nearest');

		md.geometry.thickness=md.geometry.thickness-deltaH;
		pos=find(md.geometry.thickness<1); md.geometry.thickness(pos)=1;
		md.geometry.surface=(1-di)*md.geometry.thickness;
		md.geometry.base=-di*md.geometry.thickness;

		md.materials.rheology_B=md.materials.rheology_B.*scale;
		
		md=solve(md,'Stressbalance');
		md.initialization.vx=md.results.StressbalanceSolution.Vx;
		md.initialization.vy=md.results.StressbalanceSolution.Vy;
		md.initialization.vel=md.results.StressbalanceSolution.Vel;

		penetration(j)=meanriftdeltavel(md,1);
		models{j}=md;
	end
	
	penetration_ice=penetration;
	savedata(org,penetration_ice,models,thinning);
end% }}}

%Article figures:
if perform(org,'Figure2a'),% {{{ Data

	close all;
	if loadresults,
		md=loadmodel(org,'ControlPreThermal');
		load Data/colors.mat
		%md=googlemaps(md);
		md=radarpower(md);
	end
	close all;
	
	plotmodel(md,'data',md.results.StressbalanceSolution.Vel-md.inversion.vel_obs,'data',1e-8*md.materials.rheology_B,'caxis',[-50 50],'fontsize#all',fontsize+2,'figposition','halfright','nlines',2,'ncols',1,'axis#all','square off image','googlemapsoff#all','on','overlay#all','on','xlim#all',[-2.35 -1.99]*1e6,'ylim#all',[.97 1.3]*1e6,...
	'gridded#all',1,...
	'shpdisp#all',{riftpassive,rift},'shpstyle#all',{'k-','g-'},'linewidth#all',2,'alpha#all',2,'contrast#all',2,...
	'colorbarpos',[.67 .785 .015 .1],'colorbartitle','V-V_{obs}(m/yr)','FontColor','w',...
	'colorbarpos',[.67 .465 .015 .1],'colorbartitle','B(1e8 Pa.s^{1/3})','FontColor','w',...
	'colorbarfontsize#all',15,...
	'colormap#all',lbmap(62,'Redblue'));

%.39
	dy=.15;

	subplot(2,1,1);
	text(-2.33*1e6,1.27*1e6,15,'a','Color','k','FontSize',35)
	text(-2.063*1e6,1.15*1e6,15,'1','Color','g','FontSize',20)
	text(-2.11*1e6,1.17*1e6,15,'1a','Color','k','FontSize',20)
	text(-2.07*1e6,1.13*1e6,15,'2','Color','g','FontSize',20)
	text(-2.075*1e6,1.11*1e6,15,'2a','Color','k','FontSize',20)
	text(-2.14*1e6,1.128*1e6,15,'2b','Color','k','FontSize',20)
	text(-2.092*1e6,1.095*1e6,15,'2c','Color','k','FontSize',20)
	text(-2.132*1e6,1.1*1e6,15,'3','Color','g','FontSize',20)
	text(-2.25*1e6,1.255*1e6,15,'4','Color','g','FontSize',20)

	subplot(2,1,2); a=get(gca,'Position'); a(2)=a(2)+dy; set(gca,'YDir','normal','Position',a);
	text(-2.33*1e6,1.27*1e6,15,'b','Color','k','FontSize',35)

	text(-2.063*1e6,1.15*1e6,15,'1','Color','g','FontSize',20)
	text(-2.11*1e6,1.17*1e6,15,'1a','Color','k','FontSize',20)
	text(-2.07*1e6,1.13*1e6,15,'2','Color','g','FontSize',20)
	text(-2.075*1e6,1.11*1e6,15,'2a','Color','k','FontSize',20)
	text(-2.14*1e6,1.128*1e6,15,'2b','Color','k','FontSize',20)
	text(-2.092*1e6,1.095*1e6,15,'2c','Color','k','FontSize',20)
	text(-2.132*1e6,1.1*1e6,15,'3','Color','g','FontSize',20)
	text(-2.25*1e6,1.255*1e6,15,'4','Color','g','FontSize',20)

	arrow(-2.19*1e6,1.16*1e6,-2.1102*1e6,1.1639*1e6,'ratio',.1,'widthratio',.01);
	arrow(-2.19*1e6,1.16*1e6,-2.1207*1e6,1.1513*1e6,'ratio',.1,'widthratio',.02);
	arrow(-2.19*1e6,1.16*1e6,-2.1322*1e6,1.1376*1e6,'ratio',.1,'widthratio',.03);
	text(-2.23*1e6,1.16*1e6,15,{'Suture','Zone'},'Color','k','FontSize',16)
	
	if printresults,
		export_fig(['Article/figures/' org.steps(org.currentstep).string '.pdf']);
	end


end% }}}
if perform(org,'Figure2b'),% {{{ Setup a)

	close all;

	md=loadmodel(org,'Parameterization');
	xl=[-2.12 -2.03]*1e6;
	yl=[1.12 1.16]*1e6;

	A=[-2.132    1.2135]*1e6;
	B=[-1.999    1.2135]*1e6;

	profile=[tempname '.exp'];
	a.x=[xl(1), xl(2), xl(2), xl(1) ,xl(1)];
	a.y=[yl(1), yl(1), yl(2) ,yl(2) ,yl(1)];
	a.density=1; 
	a.name='temp'; 
	expwrite(a,profile);

	plotmodel(md,'data','rifts','offset',500,'figposition','halfright','title','',...
	'figure',2,'scaleruler',[-2.051*1e6,1.155*1e6,10000,1000,3],...
	'axis','image off','expdisp',profile,'expstyle','r','linewidth',2);
	h2=gca;

	text(-2.063*1e6,1.15*1e6,15,'1','Color','g','FontSize',35)
	text(-2.07*1e6,1.13*1e6,15,'2','Color','g','FontSize',35)

	

	plotmodel(md,'data','rifts','offset',1500,'figposition','halfright','title','','figure',1,'axis','image off',...
	'scaleruler',[-2.32*1e6,1.01*1e6,60000,2000,3],'expdisp',profile,'expstyle','r','linewidth',2);

	text(-2.063*1e6,1.15*1e6,15,'1','Color','g','FontSize',35)
	text(-2.07*1e6,1.13*1e6,15,'2','Color','g','FontSize',35)
	text(-2.132*1e6,1.1*1e6,15,'3','Color','g','FontSize',35)
	text(-2.25*1e6,1.255*1e6,15,'4','Color','g','FontSize',35)
	
	plot([A(1) xl(1)],[A(2) yl(2)],'r.-.','LineWidth',2);
	plot([B(1) xl(2)],[B(2) yl(2)],'r.-.','LineWidth',2);

	ax=axes('position',[.6 .67 .3 .18]);  axis off;
	copyobj(allchild(h2),ax); 
	xlim(xl);
	ylim(yl);
	h1=gca;

	export_fig(['Article/figures/' org.steps(org.currentstep).string '.png']);

end% }}}
if perform(org,'Figure2c'),% {{{ Setup b)

	close all;
	md=loadmodel(org,'Parameterization');

	plotmodel(md,'data',md.geometry.thickness<250,'colorbar','off','figposition','halfright','axis','image off',...
	'shpdisp#all',{rift},'shpstyle#all',{'g-'},'linewidth#all',3);

	text(-2.063*1e6,1.15*1e6,15,'1','Color','g','FontSize',35)
	text(-2.07*1e6,1.13*1e6,15,'2','Color','g','FontSize',35)
	text(-2.132*1e6,1.1*1e6,15,'3','Color','g','FontSize',35)
	text(-2.25*1e6,1.255*1e6,15,'4','Color','g','FontSize',35)

	export_fig(['Article/figures/' org.steps(org.currentstep).string '.png']);


end% }}}
if perform(org,'Figure3a'),% {{{ Delta vel rotation

	close all;
	figure(1),clf,hold on;

	if loadresults,
		loaddata(org,'ThinningCurve'); models_std=models; 
		loaddata(org,'ThinningCurveIceShelfOnly'); models_ice=models; 
		loaddata(org,'ThinningCurveMelangeOnly'); models_mel=models; 

		%needed: 
		md=models_std{1};

		ind=[11 30 39];
		vels_std=zeros(md.mesh.numberofvertices,2,length(ind));
		vels_mel=zeros(md.mesh.numberofvertices,2,length(ind));
		vels_ice=zeros(md.mesh.numberofvertices,2,length(ind));
		for j=1:length(ind),
			i=ind(j);

			%std:  {{{
			md=models_std{i};
			rifti1=md.rifts.riftstruct(1);
			rifti2=md.rifts.riftstruct(4);
			penaltypairs=[rifti1.penaltypairs(:,[1 2]); fliplr(rifti2.penaltypairs(:,[1 2]))];
			segments=[rifti1.segments;fliplr(rifti2.segments)];

			%set as NaN all velocities not on rifts
			u=NaN*ones(md.mesh.numberofvertices,1);
			v=NaN*ones(md.mesh.numberofvertices,1);
			u(penaltypairs(:,1))=md.initialization.vx(penaltypairs(:,1))-md.initialization.vx(penaltypairs(:,2));
			v(penaltypairs(:,1))=md.initialization.vy(penaltypairs(:,1))-md.initialization.vy(penaltypairs(:,2));

			%process data and model
			[x y z elements is2d isplanet]=processmesh(md,[],pairoptions);
			[vel datatype]=processdata(md,[u v],pairoptions);

			%pos=find(md.mesh.x>-2.09*1e6); vel(pos,:)=0;
			vels_std(:,:,j)=vel;
			%}}}
			%mel:  {{{
			md=models_mel{i};
			rifti1=md.rifts.riftstruct(1);
			rifti2=md.rifts.riftstruct(4);
			penaltypairs=[rifti1.penaltypairs(:,[1 2]); fliplr(rifti2.penaltypairs(:,[1 2]))];
			segments=[rifti1.segments;fliplr(rifti2.segments)];

			%set as NaN all velocities not on rifts
			u=NaN*ones(md.mesh.numberofvertices,1);
			v=NaN*ones(md.mesh.numberofvertices,1);
			u(penaltypairs(:,1))=md.initialization.vx(penaltypairs(:,1))-md.initialization.vx(penaltypairs(:,2));
			v(penaltypairs(:,1))=md.initialization.vy(penaltypairs(:,1))-md.initialization.vy(penaltypairs(:,2));

			%process data and model
			[x y z elements is2d isplanet]=processmesh(md,[],pairoptions);
			[vel datatype]=processdata(md,[u v],pairoptions);

			%pos=find(md.mesh.x>-2.09*1e6); vel(pos,:)=0;
			vels_mel(:,:,j)=vel;
			%}}}
			%ice:  {{{
			md=models_ice{i};
			
			rifti1=md.rifts.riftstruct(1);
			rifti2=md.rifts.riftstruct(4);
			penaltypairs=[rifti1.penaltypairs(:,[1 2]); fliplr(rifti2.penaltypairs(:,[1 2]))];
			segments=[rifti1.segments;fliplr(rifti2.segments)];

			%set as NaN all velocities not on rifts
			u=NaN*ones(md.mesh.numberofvertices,1);
			v=NaN*ones(md.mesh.numberofvertices,1);
			u(penaltypairs(:,1))=md.initialization.vx(penaltypairs(:,1))-md.initialization.vx(penaltypairs(:,2));
			v(penaltypairs(:,1))=md.initialization.vy(penaltypairs(:,1))-md.initialization.vy(penaltypairs(:,2));

			%process data and model
			[x y z elements is2d isplanet]=processmesh(md,[],pairoptions);
			[vel datatype]=processdata(md,[u v],pairoptions);

			%pos=find(md.mesh.x>-2.09*1e6); vel(pos,:)=0;
			vels_ice(:,:,j)=vel;
			%}}}

		end

		vmax=0;
		for i=1:length(ind),
			vel=vels_std(:,:,i);
			vmax=max(vmax,max(sqrt(vel(:,1).^2+vel(:,2).^2)));
			vel=vels_ice(:,:,i);
			vmax=max(vmax,max(sqrt(vel(:,1).^2+vel(:,2).^2)));
			vel=vels_mel(:,:,i);
			vmax=max(vmax,max(sqrt(vel(:,1).^2+vel(:,2).^2)));
		end
		for i=1:length(ind),
			vels_std(1,1,i)=vmax; vels_std(1,2,i)=0;
			vels_mel(1,1,i)=vmax; vels_mel(1,2,i)=0;
			vels_ice(1,1,i)=vmax; vels_ice(1,2,i)=0;
		end 
	end


	plotmodel(md,'nlines',3,'ncols',3,...
	'data',vels_std(:,:,1),...
	'data',vels_std(:,:,2),...
	'data',vels_std(:,:,3),...
	'data',vels_mel(:,:,1),...
	'data',vels_mel(:,:,2),...
	'data',vels_mel(:,:,3),...
	'data',vels_ice(:,:,1),...
	'data',vels_ice(:,:,2),...
	'data',vels_ice(:,:,3),...
	'ylim#all',[1.135 1.15]*1e6,'xlim#all',[    -2120000    -2050000],...
	'shpdisp#all',{domain,rift},'shpstyle#all',{'k','k'},'axis#all','off image',...
	'colorbar#1,2,3,4,5,6,7,8','off','fontsize#all',fontsize-8,'figposition','halfright','colorbarpos#9',[.91 .685 .015 .1],'scaling#all',1.2);
	

	fontsizet=20;
	subplot(3,3,1); 
	text(-2.098*1e6,1.15*1e6,sprintf('\\DeltaV_{avg}=%2.2g (m/yr)',penetration_std(ind(1))),'fontsize',fontsizet-5); 
	text(-2.145*1e6,1.1425*1e6,{'Ice Shelf','+','M\''elange'},'FontSize',fontsizet,'interpreter','latex'); 
	text(-2.09*1e6,1.165*1e6,sprintf('\\DeltaH=%2.2g (m)',thinning(ind(1))),'FontSize',fontsizet);
	text(-2.049*1e6,1.14*1e6,15,'1','Color','g','FontSize',20); text(-2.07*1e6,1.135*1e6,15,'2','Color','g','FontSize',20);

	subplot(3,3,2); 
	text(-2.098*1e6,1.15*1e6,sprintf('\\DeltaV_{avg}=%2.2g (m/yr)',penetration_std(ind(2))),'FontSize',fontsizet-5);
	text(-2.09*1e6,1.165*1e6,sprintf('\\DeltaH=%2.2g (m)',thinning(ind(2))),'FontSize',fontsizet);
	text(-2.049*1e6,1.14*1e6,15,'1','Color','g','FontSize',20); text(-2.07*1e6,1.135*1e6,15,'2','Color','g','FontSize',20);

	subplot(3,3,3); 
	text(-2.098*1e6,1.15*1e6,sprintf('\\DeltaV_{avg}=%2.2g (m/yr)',penetration_std(ind(3))),'FontSize',fontsizet-5);
	text(-2.09*1e6,1.165*1e6,sprintf('\\DeltaH=%2.2g (m)',thinning(ind(3))),'FontSize',fontsizet);
	text(-2.049*1e6,1.14*1e6,15,'1','Color','g','FontSize',20); text(-2.07*1e6,1.135*1e6,15,'2','Color','g','FontSize',20);
	
	subplot(3,3,4); 
	text(-2.098*1e6,1.15*1e6,sprintf('\\DeltaV_{avg}=%2.2g (m/yr)',penetration_mel(ind(1))),'FontSize',fontsizet-5); 
	text(-2.145*1e6,1.1425*1e6,'M\''elange','FontSize',fontsizet,'interpreter','latex');
	text(-2.049*1e6,1.14*1e6,15,'1','Color','g','FontSize',20); text(-2.07*1e6,1.135*1e6,15,'2','Color','g','FontSize',20);
	
	subplot(3,3,5);
	text(-2.098*1e6,1.15*1e6,sprintf('\\DeltaV_{avg}=%2.2g (m/yr)',penetration_mel(ind(2))),'FontSize',fontsizet-5);
	text(-2.049*1e6,1.14*1e6,15,'1','Color','g','FontSize',20); text(-2.07*1e6,1.135*1e6,15,'2','Color','g','FontSize',20);
	
	subplot(3,3,6);
	text(-2.098*1e6,1.15*1e6,sprintf('\\DeltaV_{avg}=%2.3g (m/yr)',penetration_mel(ind(3))),'FontSize',fontsizet-5);
	text(-2.049*1e6,1.14*1e6,15,'1','Color','g','FontSize',20); text(-2.07*1e6,1.135*1e6,15,'2','Color','g','FontSize',20);
	
	subplot(3,3,7); 
	text(-2.098*1e6,1.15*1e6,sprintf('\\DeltaV_{avg}=%2.2g (m/yr)',penetration_ice(ind(1))),'FontSize',fontsizet-5);
	text(-2.145*1e6,1.1425*1e6,'Ice Shelf','FontSize',fontsizet,'interpreter','latex');
	text(-2.049*1e6,1.14*1e6,15,'1','Color','g','FontSize',20); text(-2.07*1e6,1.135*1e6,15,'2','Color','g','FontSize',20);
	
	subplot(3,3,8);
	text(-2.098*1e6,1.15*1e6,sprintf('\\DeltaV_{avg}=%2.2g (m/yr)',penetration_ice(ind(2))),'FontSize',fontsizet-5);
	text(-2.049*1e6,1.14*1e6,15,'1','Color','g','FontSize',20); text(-2.07*1e6,1.135*1e6,15,'2','Color','g','FontSize',20);
	
	subplot(3,3,9);
	text(-2.098*1e6,1.15*1e6,sprintf('\\DeltaV_{avg}=%2.2g (m/yr)',penetration_ice(ind(3))),'FontSize',fontsizet-5);
	text(-2.049*1e6,1.14*1e6,15,'1','Color','g','FontSize',20); text(-2.07*1e6,1.135*1e6,15,'2','Color','g','FontSize',20);

	dx=.05; dy=.22;
	subplot(3,3,1); a=get(gca,'Position'); a(1)=a(1)+dx; set(gca,'Position',a);
	subplot(3,3,3); a=get(gca,'Position'); a(1)=a(1)-dx; set(gca,'Position',a);
	subplot(3,3,4); a=get(gca,'Position'); a(1)=a(1)+dx; a(2)=a(2)+dy; set(gca,'Position',a); 
	subplot(3,3,5); a=get(gca,'Position'); a(2)=a(2)+dy; set(gca,'Position',a);
	subplot(3,3,6); a=get(gca,'Position'); a(1)=a(1)-dx; a(2)=a(2)+dy; set(gca,'Position',a); 
	subplot(3,3,7); a=get(gca,'Position'); a(1)=a(1)+dx; a(2)=a(2)+2*dy; set(gca,'Position',a); 
	subplot(3,3,8); a=get(gca,'Position'); a(2)=a(2)+2*dy; set(gca,'Position',a);
	subplot(3,3,9); a=get(gca,'Position'); a(1)=a(1)-dx; a(2)=a(2)+2*dy; set(gca,'Position',a); 

	children=get(gcf,'Children'); colorb=children(1); tit=get(colorb,'title'); set(tit,'String','\DeltaV(m/yr)','FontSize',fontsizet); set(colorb,'FontSize',fontsizet-3);

	export_fig(['Article/figures/' org.steps(org.currentstep).string '.pdf']);

end
% }}}
if perform(org,'Figure3b'),% {{{ %OpeningRate vs Thinning

	close all	
	if loadresults,
		loaddata(org,'ThinningCurve');
		loaddata(org,'ThinningCurveIceShelfOnly');
		loaddata(org,'ThinningCurveMelangeOnly');
	end

	figure(1),clf,hold on;
	set(gcf,'Position',[584         603        1435         582]);
	p=plot(thinning,penetration_std,'r',thinning,penetration_mel,'b',thinning,penetration_ice,'k');
	set(p,'LineWidth',linewidth);
	ylabel('\DeltaV_{avg} (m/yr)');
	xlabel('Thinning (m)');
	l=legend('Shelf + M\''elange','M\''elange','Shelf','interpreter','latex');
	set(l,'Location','NorthWest');
	grid on;set(gca,'FontSize',fontsize+15); axis tight;
	set(gcf,'color','w');

	export_fig(['Article/figures/' org.steps(org.currentstep).string '.pdf'],'-nocrop');
	
end% }}}
