steps   = [21];
glacier = 'Denman';
clustername = oshostname;
%clustername = 'pfe';
%clustername = 'greenplanet';
%clustername = 'discovery'; %needs to be uncommented to run parameter sweep
loadonly = 0;

%default settings
ss_sigma_max = 0.1e6; %UNIT: MPa.
ss_eigen=9e7; %UNIT: m*a
ss_min_thickness = 200; %UNIT: m
ss_rc = 0.6;
ss_migmax = 10000; %UNIT: m

sigma_max_loop = 250000:5e3:350000;
eigen_loop = 1e8:0.5e8:1e9;
min_thickness_loop = 350:5:450;
rc_loop = 0.1:0.1:0.9;
migmax_loop = 2000:2000:10000;

%EXAMPLE: to run a parameter sweep for Ross on Discovery, enter in the following commands in Matlab
%for ss_sigma_max = 50000:5e3:150000; runme; end
%for ss_eigen = 1e8:0.5e8:1e9; runme; end
%for ss_min_thickness = 200:5:300; runme; end
%for ss_rc = 0.1:0.1:0.9; for ss_migmax = 2000:2000:10000; runme; end; end


%Cluster parameters{{{
if strcmpi(clustername,'pfe'),
	interactive = 0;
	%cluster=pfe('numnodes',8,'time',30,'interactive',interactive,'processor','har','queue','devel');
	cluster=pfe('numnodes',8,'time',30,'processor','wes','queue','devel');
	%cluster=pfe('numnodes',8,'time',30,'processor','wes','queue','debug');
	%cluster=pfe('numnodes',8,'time',30,'interactive',interactive,'processor','har','queue','alphatst@pbspl233');
elseif strcmpi(clustername,'astrid'),
	cluster=generic('name',oshostname(),'np',12);
	%cluster=generic('name',oshostname(),'np',1);
elseif strcmpi(clustername,'discovery')
      cluster=discovery('numnodes',1,'cpuspernode',16);
      cluster.time = 1; %hr
else
	cluster=generic('name',oshostname(),'np',12);
end
clear interactive
%}}}
org=organizer('repository',['./Models'],'prefix',['Model_' glacier '_'],'steps',steps); clear steps;
%addpath([jplsvn() '/proj-morlighem/CODE/']);
fprintf(['\n  ========  ' upper(glacier) '  ========\n\n']);

%%Steps 1-5
if perform(org,'Mesh'),% {{{

	%hmin = 3e3; %highest allowed resolution - DEFAULT
	if strcmp(glacier,'Thwaites')||strcmp(glacier,'Denman')
		hmin = 1.0e3; %Subject to change
	else
		hmin = 2.0e3;
	end

	md=triangle(model,['Exp/' glacier '.exp'],hmin);
	
	%refine mesh based on speed
	[vx vy] = interpMouginotAnt2019(md.mesh.x,md.mesh.y);
	vel  = sqrt(vx.^2+vy.^2);

	%Also increase resolution around the ice front (2 km)
	M=interpBedmachineAntarctica(md.mesh.x,md.mesh.y,'mask');
	ice_levelset = -ones(md.mesh.numberofvertices,1);
	ice_levelset(M==0) = +1;
	ice_levelset = reinitializelevelset(md,ice_levelset);
	hVertices = NaN(md.mesh.numberofvertices,1);
	pos = find(abs(ice_levelset)<100e3);
	hVertices(pos) = hmin;
	%hVertices(abs(ice_levelset)<100e3) = hmin;

	md=bamg(md,'hmin',hmin,'hmax',100e3,'field',vel,'err',9,'hVertices',hVertices);

	[md.mesh.lat,md.mesh.long]  = xy2ll(md.mesh.x,md.mesh.y,-1,0,71);
	md.mesh.epsg=3031;

	if strcmp(glacier,'Thwaites')
		md.mesh.scale_factor=(1+sin(md.mesh.lat*pi/180))/(1+sin(-71*pi/180));
		md.miscellaneous.name='mesh';
	end

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

	md=loadmodel(org,'Mesh');
	md=setflowequation(md,'SSA','all');
	
	if strcmpi(glacier,'Thwaites')
		md=parameterize(md,'Antarctica_Thwaites.par');
	else
		md=parameterize(md,'Antarctica.par');
	end

	%adjust mask for some domains
	if strcmpi(glacier, 'Thwaites')
		pos = find(ContourToNodes(md.mesh.x,md.mesh.y,'Exp/ThwaitesRemove.exp',2));
		md.mask.ice_levelset(pos) = +1;
		md.mask.ice_levelset = reinitializelevelset(md,md.mask.ice_levelset);
		md.materials.rheology_law = 'Cuffey'; 
	elseif strcmpi(glacier,'Denman')
		pos = find(md.mask.ocean_levelset>0 & interpBedmachineAntarctica(md.mesh.x,md.mesh.y,'mask')==3); %find all points that are marked as land and also as floating ice
		md.mask.ocean_levelset(pos)= -1; %and set these points to be ocean
		%md.geometry.bed(pos) = md.geometry.bed(pos) -1000; %Decrease the bed underneath floating ice
		pos = find(ContourToNodes(md.mesh.x,md.mesh.y,'Exp/DenmanRemove.exp',2));
		md.mask.ocean_levelset(pos) = -1; %force the island as ocean
		md.mask.ocean_levelset = reinitializelevelset(md,md.mask.ocean_levelset);
	end
	savemodel(org,md);
end%}}}
if perform(org,'InversionB'),% {{{

	if strcmpi(glacier,'Thwaites')==0
		md=loadmodel(org,'Param');

		%Control general
		md.inversion=m1qn3inversion(md.inversion);

		md.inversion.iscontrol=1;
		md.verbose=verbose('solution',false,'control',true);

		% Cost functions
		md.inversion.cost_functions=[101 103 502];
		md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,3);
		md.inversion.cost_functions_coefficients(:,1)=1000;
		md.inversion.cost_functions_coefficients(:,2)=1;
		md.inversion.cost_functions_coefficients(:,3)=1.e-17;

		%Fill in blanks in velocity data
		pos = find(isnan(md.inversion.vx_obs) | isnan(md.inversion.vy_obs));
		md.inversion.vx_obs(pos) = 0;
		md.inversion.vy_obs(pos) = 0;
		md.inversion.cost_functions_coefficients(pos,[1 2]) = 0;
		md.initialization.vx = md.inversion.vx_obs;
		md.initialization.vy = md.inversion.vy_obs;

		% Controls
		md.inversion.control_parameters={'MaterialsRheologyBbar'};
		md.inversion.min_parameters=md.materials.rheology_B;
		md.inversion.max_parameters=md.materials.rheology_B;
		md.inversion.gttol = 1e-8;
		md.inversion.maxsteps = 40;
		md.inversion.maxiter  = 40;
		pos = find(md.mask.ocean_levelset<0);
		md.inversion.min_parameters(pos) = cuffey(273);
		md.inversion.max_parameters(pos) = cuffey(200);

		%Additional parameters
		md.stressbalance.restol=0.01;
		md.stressbalance.reltol=0.1;
		md.stressbalance.abstol=NaN;

		%Go solve
		md.cluster=cluster;
		mds=extract(md,md.mask.ocean_levelset<0);
		%mds.toolkits.DefaultAnalysis=bcgslbjacobioptions();
		mds=solve(mds,'Stressbalance');

		% Update model rheology_B accordingly
		md.materials.rheology_B(mds.mesh.extractedvertices)=mds.results.StressbalanceSolution.MaterialsRheologyBbar;

		%Extrapolate outwards
		flags = zeros(md.mesh.numberofvertices,1);
		flags(md.mask.ocean_levelset<0 & md.mask.ice_levelset<0) = 1;
		pos1 = find(flags); pos2 = find(~flags);
		md.materials.rheology_B(pos2) = griddata(md.mesh.x(pos1),md.mesh.y(pos1), md.materials.rheology_B(pos1), md.mesh.x(pos2),md.mesh.y(pos2), 'nearest');


	% For Thwaites, set InversionB to that of of the runme that works for Mathieu's MICI code. 
	elseif strcmpi(glacier, 'Thwaites')
		md=loadmodel(org,'Param');

		%Control general
		md.inversion=m1qn3inversion(md.inversion);

		md.inversion.iscontrol=1;
		md.verbose=verbose('solution',false,'control',true);


		% Set inversion data
		md.inversion.vx_obs=md.initialization.vx; % initialization was defined in last step (raw data, with NaN)
		md.inversion.vy_obs=md.initialization.vy; % initialization was defined in last step (raw data, with NaN)

		pos=find(isnan(md.inversion.vx_obs) | isnan(md.inversion.vy_obs));
		md.inversion.vx_obs(pos)=0;
		md.inversion.vy_obs(pos)=0;
		md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
		md.initialization.vx(pos)=0;
		md.initialization.vy(pos)=0;
		md.initialization.vel(pos)=0;

		% Control general
		md.inversion.iscontrol=1;
		md.inversion.maxsteps=40;
		md.inversion.maxiter=40;
		md.inversion.dxmin=0.1;
		md.inversion.gttol=1.0e-6;
		md.inversion.incomplete_adjoint=0; % 0: non linear viscosity, 1: linear viscosity 04/29/2019 changed to non linear

		% Cost functions
		md.inversion.cost_functions=[101 103 502];
		md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,length(md.inversion.cost_functions));
		md.inversion.cost_functions_coefficients(:,2)=1; % always 1
		md.inversion.cost_functions_coefficients(:,end)=1.e-18; %1.e-19;%1.e-18;; % this may be changed according to the L-curve analysis
		md.inversion.cost_functions_coefficients(pos,:)=0; % positions with NaN in the velocity data set

		if 1 % keep this here (not in inversion for C)
			% Force "mounds" as Dirichlet BC
			% This avoids problems during the stressbalance convergence solver
			incontour=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,expread('./Exp/Mound_region_1995_correction_rev00.exp'),'node',1);
			pos=find(incontour);
			md.stressbalance.spcvx(pos)=md.initialization.vx(pos);
			md.stressbalance.spcvy(pos)=md.initialization.vy(pos);
			md.stressbalance.spcvz(pos)=md.initialization.vz(pos);
			% second region
			incontour=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,expread('./Exp/Mound_region2_1995_correction_rev00.exp'),'node',1);
			pos=find(incontour);
			md.stressbalance.spcvx(pos)=md.initialization.vx(pos);
			md.stressbalance.spcvy(pos)=md.initialization.vy(pos);
			md.stressbalance.spcvz(pos)=md.initialization.vz(pos);
			% third region (floating)
			incontour=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,expread('./Exp/Floating_region_1995_correction_rev00.exp'),'node',1);
			pos=find(incontour);
			md.mask.ice_levelset(pos)=1; % forcing no ice
		end

		% Controls
		% setting initial guess for rheology B
		Dmax=0.95; % maximum amount of damage
		md.inversion.control_parameters={'MaterialsRheologyBbar'};
		%md.inversion.min_parameters=(1-Dmax)*md.materials.rheology_B; % original, firstly tested
		%md.inversion.max_parameters=15*md.materials.rheology_B; % original, firstly tested
		% testing new value for min and max
		md.inversion.min_parameters=cuffey(273.15)*ones(size(md.materials.rheology_B)); % from Seroussi et al, 2014
		md.inversion.max_parameters=cuffey(273.15-30)*ones(size(md.materials.rheology_B)); % from Seroussi et al, 2014

		% Additional parameters
		md.stressbalance.restol=0.0001; % 04/29/2019
		md.stressbalance.reltol=0.01; % 04/29/2019
		md.stressbalance.abstol=10; % 04/29/2019
		md.stressbalance.maxiter=40; % 08/21/2019

		% Prepare to solve
		md.cluster=cluster;
		md.verbose=verbose('solution',false,'control',true);
		md.miscellaneous.name='inversion_B';
		mds=extract(md,md.mask.ocean_levelset<0);
		mds.friction.coefficient(:)=0; % make sure there is no basal friction

		% Solve
		%mds.toolkits.DefaultAnalysis=bcgslbjacobioptions();% biconjugate gradient with block Jacobi preconditioner
		mds.settings.solver_residue_threshold=NaN; % 11/05/2019
		mds.stressbalance.maxiter=50; % 10/24/2019
		mds.stressbalance.reltol=NaN; % 11/05/2019
		mds.stressbalance.abstol=NaN; % 11/05/2019

		mds=solve(mds,'Stressbalance'); % only extracted model

		% Update model rheology_B accordingly
		md.materials.rheology_B(mds.mesh.extractedvertices)=mds.results.StressbalanceSolution.MaterialsRheologyBbar;
	end
	savemodel(org,md);
end%}}}
if perform(org,'DragInversion'),% {{{
	if strcmpi(glacier,'Thwaites')==0
		md = loadmodel(org, 'InversionB');
		
		% Cost functions
		md.inversion.cost_functions=[101 103 501];
		md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,3);
		md.inversion.cost_functions_coefficients(:,1)=2000;
		md.inversion.cost_functions_coefficients(:,2)=1;
		md.inversion.cost_functions_coefficients(:,3)=8e-7;

		% Controls
		md.inversion.control_parameters={'FrictionCoefficient'};
		md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
		md.inversion.max_parameters=300*ones(md.mesh.numberofvertices,1);

		%Change default solver to increase performance
		md.toolkits.DefaultAnalysis=bcgslbjacobioptions();
		md.settings.solver_residue_threshold = 1e-5;

		%TEMP
		md=sethydrostaticmask(md);

		% Solve
		md=solve(md,'Stressbalance');

		% Update model friction fields accordingly
		md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient;


		savemodel(org,md);
	elseif strcmpi(glacier,'Thwaites')
		md.friction.p=ones(md.mesh.numberofelements,1);
      md.friction.q=ones(md.mesh.numberofelements,1);
      md.inversion.cost_functions=[101 103 501];
      md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,3); % IMPORTANT: set as 0 again the vertices with no data (NaN)
      md.inversion.cost_functions_coefficients(:,1)=2.2452e+03;%2.7526e+03; % for budd change this such that the cost functions 101 and 103 have the same contribution
      md.inversion.cost_functions_coefficients(:,2)=1; % always 1
      md.inversion.cost_functions_coefficients(:,end)=1e-8;  % for budd
      md.inversion.cost_functions_coefficients(pos,:)=0; % positions with NaN in the velocity data set
      %Controls
      md.friction.coefficient=EstimateFric_Budd(md); % initial guess from Driving Stress (using r=1 and s=1)
      md.inversion.control_parameters={'FrictionCoefficient'};
      md.inversion.min_parameters=0*ones(md.mesh.numberofvertices,1);
      %md.inversion.max_parameters=3000*ones(md.mesh.numberofvertices,1);
      md.inversion.max_parameters=300*ones(md.mesh.numberofvertices,1);
      md.miscellaneous.name='inversion_drag_budd';

		md.inversion.iscontrol=1;
		md.inversion.maxsteps=80;
		md.inversion.maxiter=80;
		md.inversion.dxmin=0.01;
		md.inversion.gttol=1.0e-6;
		md.inversion.incomplete_adjoint=0; % 0: non linear viscosity, 1: linear viscosity 04/29/2019 changed to non linear

		% Keep friction 0 in the floating part
		pos=find(md.mask.ocean_levelset<0);
		md.friction.coefficient(pos)=0;
		
		md.stressbalance.restol=0.002; % 04/29/2019
		md.stressbalance.reltol=0.01; % 04/29/2019
		md.stressbalance.abstol=10; % 04/29/2019
		md.stressbalance.maxiter=30; % 04/29/2019

		% Set cluster
		md.cluster=cluster;
		md.verbose=verbose('solution',false,'control',true);

		% Define solver
		%md.toolkits.DefaultAnalysis=bcgslbjacobioptions();% biconjugate gradient with block Jacobi preconditioner
		%md.settings.solver_residue_threshold=1000*1.0000e-06; %relaxing the threshold
		md.settings.solver_residue_threshold=NaN; % 11/05/2019
		md.stressbalance.maxiter=50; % 10/24/2019
		md.stressbalance.reltol=NaN; % 11/05/2019
		md.stressbalance.abstol=NaN; % 11/05/2019

		md.settings.solver_residue_threshold = 1e-5;
		md=solve(md,'Stressbalance');

		md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient;
		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
end%}}}
if perform(org,'TransientPrep'),% {{{

	md=loadmodel(org,'DragInversion');

	md.initialization.pressure = zeros(md.mesh.numberofvertices,1); %FIXME
	md.initialization.temperature = 250*ones(md.mesh.numberofvertices,1); %FIXME
	%set up mass transport BC
	md.masstransport.spcthickness = NaN(md.mesh.numberofvertices,1);
	pos = find(md.mesh.vertexonboundary);
	md.masstransport.spcthickness(pos) = md.geometry.thickness(pos);

	% Set parameters
	md.inversion.iscontrol=0;
	if 0
		md.timestepping.start_time = 0;
		md.timestepping.time_step  = .5;
		md.timestepping.final_time = 200; %200 years
	else
		md.timestepping = timesteppingadaptive();
		md.timestepping.start_time = 0;
		md.timestepping.final_time = 200; %200 years
		if strcmpi(glacier,'Thwaites')||strcmpi(glacier,'Denman')
			md.timestepping.time_step_min = .01/3; %Edit 11/21/22
		else
			md.timestepping.time_step_min = .01/2; %Edit 11/28/22
		end
		md.timestepping.time_step_max = 1;
	end

	if strcmpi(glacier,'Thwaites')||strcmpi(glacier,'Denman')
		md.settings.output_frequency=8; %Edit 11/21/22
	end

	% We set the transient parameters
	md.transient.ismovingfront   = 1;
	md.transient.isthermal       = 0;
	md.transient.isstressbalance = 1;
	md.transient.ismasstransport = 1;
	md.transient.isgroundingline = 1;
	md.groundingline.migration   = 'SubelementMigration';

	%Levelset parameters
	if strcmpi(glacier,'Denman')
		pos = find(md.mask.ocean_levelset>0 & interpBedmachineAntarctica(md.mesh.x,md.mesh.y,'mask')==3); %find all point that are marked as land and also as floating ice
		md.mask.ocean_levelset(pos)= -1; %and set these points to be ocean
	end
	md.mask.ice_levelset = reinitializelevelset(md,md.mask.ice_levelset);
	md.levelset.stabilization = 2;
	md.levelset.spclevelset= NaN(md.mesh.numberofvertices,1);
   pos = find(md.mesh.vertexonboundary);
   md.levelset.spclevelset(pos) = md.mask.ice_levelset(pos);
	
	%less melt along the sides
	md.frontalforcings.meltingrate = zeros(md.mesh.numberofvertices,1);

	% Set the requested outputs
	
	%md.transient.requested_outputs={'default','IceVolume','Calvingratex','Calvingratey','CalvingCalvingrate'}; %For rate-based calving laws
	%md.transient.requested_outputs={'default','IceVolume'} %For position-based calving laws

	md.toolkits.DefaultAnalysis=bcgslbjacobioptions();
	%if strcmpi(glacier,'Thwaites')
	%	md.toolkits.DefaultAnalysis=issmmumpssolver(); %EDIT 11/18/22
	%end
	md.verbose.solution=1;
	md.cluster = cluster;
	if ~isa(cluster,'pfe'),
		np = ceil(md.mesh.numberofelements/2000);
		if(md.cluster.np>np), md.cluster.np = np; end
	end

	savemodel(org,md);
end%}}}

%%Steps 6-10
if perform(org,['CalvingVM' num2str(ss_sigma_max)]),% {{{

	md=loadmodel(org,'TransientPrep');

	%Tuning parameter
	md.calving=calvingvonmises();
	md.calving.stress_threshold_floatingice = ss_sigma_max; 
	md.calving.stress_threshold_groundedice = 1e6; %1MPa
	%md.calving.min_thickness= 0;

	md.stressbalance.requested_outputs = {'Vel'};
	md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation','Calvingratex','Calvingratey','CalvingCalvingrate'};

	md.settings.solver_residue_threshold = 1e-5;
	md.levelset.stabilization = 5;
	md.levelset.reinit_frequency=10;
	md.levelset.migration_max = 10e3;
	md.cluster = cluster;

	%Make sure jobs are submitted without MATLAB waiting for the results, and job name is unique
   md.settings.waitonlock = 0;
   md.cluster.interactive = 0; %only needed if you are using the generic cluster
   md.miscellaneous.name = [glacier '_VM_' num2str(ss_sigma_max) ];


  %Submit job or download results, make sure that there is no runtime name (that includes the date)
   md=solve(md,'Transient','runtimename',false,'loadonly',loadonly);

	if loadonly
		savemodel(org,md);
	end
end%}}}
if perform(org,['CalvingLevermann' num2str(ss_eigen)]),% {{{
	md = loadmodel(org,'TransientPrep');

	md.calving=calvinglevermann();
	%Coefficient K2 in ms (need to divide by md.constants.yts to compare to Levermann's paper)
	%md.calving.coeff = 2e15*ones(md.mesh.numberofvertices,1); %2e13,2e15 advance; 2e17 retreat 

	%Add more output
	md.transient.requested_outputs={'default','StrainRateparallel','StrainRateperpendicular','IceVolume','IceVolumeAboveFloatation','Calvingratex','Calvingratey','CalvingCalvingrate'};
	md.settings.solver_residue_threshold = 1e-5;
		
	md.calving.coeff =ss_eigen*md.constants.yts *ones(md.mesh.numberofvertices,1); %2e13,2e15 advance; 2e17 retreat 
	md.levelset.stabilization = 5;
   md.levelset.reinit_frequency=10;
	md.levelset.migration_max = 10e3;
	md.cluster = cluster;

	%Make sure jobs are submitted without MATLAB waiting for the results, and job name is unique
   md.settings.waitonlock = 0;
	md.cluster.interactive = 0; %only needed if you are using the generic cluster
	md.miscellaneous.name = [glacier '_EC_' num2str(ss_eigen)];
	
	%Submit job or download results, make sure that there is no runtime name (that includes the date)
   md=solve(md,'Transient','runtimename',false,'loadonly',loadonly);
	%md=solve(md,'tr');
	%savemodel(org,md);

   if loadonly
      savemodel(org,md);
   end

end%}}}
if perform(org,'CalvingCrevasseDepth'),% {{{
   md = loadmodel(org,'TransientPrep')
	%This step is the generic Crevasse Depth law, NOT the Pollard crevasse depth law for floating ice shelves.

   md.calving=calvingcrevassedepth();
   md.calving.crevasse_opening_stress = 0;
	md.calving.crevasse_threshold = 1;
	md.calving.water_height = zeros(md.mesh.numberofvertices,1);

	md.levelset.stabilization = 5;
   md.levelset.reinit_frequency=10;
   md.cluster = cluster;
	
	md=solve(md,'tr');
   savemodel(org,md);
end%}}}
if perform(org,['CalvingMinThickness' num2str(ss_min_thickness)]),% {{{
	md = loadmodel(org,'TransientPrep');
	
	md.levelset.migration_max = 10e3;

	md.calving=calvingminthickness();
	md.calving.min_thickness = ss_min_thickness;

	md.levelset.stabilization = 5;
   md.levelset.reinit_frequency=10;
   md.cluster = cluster;

   %Make sure jobs are submitted without MATLAB waiting for the results, and job name is unique
   md.settings.waitonlock = 0;
   md.cluster.interactive = 0; %only needed if you are using the generic cluster
   md.miscellaneous.name = [glacier '_MT_' num2str(ss_min_thickness) ];

	md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation'};

   %Submit job or download results, make sure that there is no runtime name (that includes the date)
   md=solve(md,'Transient','runtimename',false,'loadonly',loadonly);
   %md=solve(md,'tr');
   %savemodel(org,md);

   if loadonly
      savemodel(org,md);
   end
   %savemodel(org,md);
end%}}}
if perform(org,'CalvingMinThickness_local'),% {{{
   md = loadmodel(org,'TransientPrep');

	%Minimum thickness step for testing on local machine (not running a parameter sweep on supercomputer)

	md.levelset.migration_max = 10e3; %restrict the amount that the leveset can migrate per year.

   md.calving=calvingminthickness();
   md.calving.min_thickness = 395;  %minimum thickness below which no ice is allowed [m]
	
   md.levelset.stabilization = 5;
   md.levelset.reinit_frequency=10;
   md.cluster = cluster;

   md=solve(md,'tr');
   savemodel(org,md);
end%}}}

%%Steps 11-15
if perform(org,'CalvingVMMT_local'),% {{{
   md = loadmodel(org,'TransientPrep');
	
	%Testing with VM and MT calving laws combined. For local machine use, not parameter sweeps for supercomputer.

   md.levelset.migration_max = 10e3; %restrict the amount that the leveset can migrate per year.
	   %Tuning parameter
   md.calving=calvingvonmises();
   md.calving.stress_threshold_floatingice = 110e3; %110 kPa
   md.calving.stress_threshold_groundedice = 1e6; %1MPa
   md.calving.min_thickness= 245; %m

   md.settings.solver_residue_threshold = 1e-5;

   md.levelset.stabilization = 2;
   md.levelset.reinit_frequency=10;
   md.cluster = cluster;

   md=solve(md,'tr');
   savemodel(org,md);
end%}}}
if perform(org,'CalvingVM_local'),% {{{
   md = loadmodel(org,'TransientPrep');

	%This step is for testing the VM calving law on the local machine, not on a supercomputer for parameter sweeps.

	md.transient.requested_outputs={'default','IceVolume','Calvingratex','Calvingratey','CalvingCalvingrate'};

      %Tuning parameter
   md.calving=calvingvonmises();
   md.calving.stress_threshold_floatingice = 190000; %110 kPa
   md.calving.stress_threshold_groundedice = 1e6; %1MPa

   md.settings.solver_residue_threshold = 1e-5;

   md.levelset.stabilization = 5;
   md.levelset.reinit_frequency=10;
   md.levelset.migration_max = 10e3;
	md.cluster = cluster;

   md=solve(md,'tr');
   savemodel(org,md);
end%}}}
if perform(org,'CalvingLevermann_local'),% {{{
   md = loadmodel(org,'TransientPrep');

	%This step is to test eigencalving on the local machine, not for the parameter sweeps on a supercomputer.

   md.calving=calvinglevermann();
   %Coefficient K2 in ms (need to divide by md.constants.yts to compare to Levermann's paper)
	calvingcoeff = 3e8;

   %Add more output
	md.transient.requested_outputs={md.transient.requested_outputs{:},'StrainRateparallel','StrainRateperpendicular','IceVolume','IceVolumeAboveFloatation','Calvingratex','Calvingratey','CalvingCalvingrate'};
   md.settings.solver_residue_threshold = 1e-5;

   md.calving.coeff =calvingcoeff*md.constants.yts *ones(md.mesh.numberofvertices,1); %2e13,2e15 advance; 2e17 retreat
   md.levelset.stabilization = 5;
   md.levelset.reinit_frequency=10;
   md.levelset.migration_max = 10e3;
   md.cluster = cluster;

	md=solve(md,'tr');
   savemodel(org,md);
end%}}}
if perform(org,'ExtractResults'),% {{{
	%This step extracts results from each run of the parameter sweep and creates a new "results" file to store the results
	for ss_sigma_max = sigma_max_loop;
		modelname = ['Models/Model_' glacier '_CalvingVM' num2str(ss_sigma_max)];
		disp(['Loading ' modelname]);
		load(modelname);
		levelset = md.results.TransientSolution(end).MaskIceLevelset;
		lsend =double(md.results.TransientSolution(end).MaskIceLevelset<0)*2-1;
		ls1 =double(md.mask.ice_levelset<0)*2-1;
		lsdiff = lsend-ls1;
		areabetween_ind = find(abs(lsdiff)>0);
		areabetween_indices = md.mesh.elements(areabetween_ind,:,:);
		areasum = sum(GetAreas(areabetween_indices,md.mesh.x,md.mesh.y));
		iso=isoline(md, md.mask.ice_levelset,'value',0);
		isolength = 0; %Change 5/3/23
		for i=1:length(iso) %Change 5/3/23
			isolength = max(isolength, sum(sqrt(diff(iso(i).x).^2 + diff(iso(i).y).^2))); %Find the largest calving front (there may be small spurious calving fronts produced)
		end
		
		firstfront=isoline(md, md.mask.ice_levelset,'output','matrix');
		lastfront=isoline(md, md.results.TransientSolution(end).MaskIceLevelset,'output','matrix');
		%frechetdist = frechet(sort(rmmissing(firstfront(:,1))),sort(rmmissing(firstfront(:,2))),sort(rmmissing(lastfront(:,1))),sort(rmmissing(lastfront(:,2)))); 
		frechetdist = DiscreteFrechetDist(sortrows(firstfront),sortrows(lastfront));
		save([modelname '_results'],'levelset','areasum','isolength','frechetdist');
	end
	for ss_eigen = eigen_loop;
	   modelname = ['Models/Model_' glacier '_CalvingLevermann' num2str(ss_eigen)];
      disp(['Loading ' modelname]);
      load(modelname);
      levelset = md.results.TransientSolution(end).MaskIceLevelset;
      lsend =double(md.results.TransientSolution(end).MaskIceLevelset<0)*2-1;
      ls1 =double(md.mask.ice_levelset<0)*2-1;
		lsdiff = lsend-ls1;
      areabetween_ind = find(abs(lsdiff)>0);
      areabetween_indices = md.mesh.elements(areabetween_ind,:,:);
      areasum = sum(GetAreas(areabetween_indices,md.mesh.x,md.mesh.y));
      iso=isoline(md, md.mask.ice_levelset,'value',0); %original
      %isolength=sum([1;sqrt(diff(iso(1).x).^2 + diff(iso(1).y).^2)]); %original
		isolength = 0; %Change 5/3/23
      for i=1:length(iso) %Change 5/3/23
         isolength = max(isolength, sum(sqrt(diff(iso(i).x).^2 + diff(iso(i).y).^2)));
      end
		firstfront=isoline(md, md.mask.ice_levelset,'output','matrix');
      lastfront=isoline(md, md.results.TransientSolution(end).MaskIceLevelset,'output','matrix');
		%frechetdist = frechet(sort(rmmissing(firstfront(:,1))),sort(rmmissing(firstfront(:,2))),sort(rmmissing(lastfront(:,1))),sort(rmmissing(lastfront(:,2))));
		frechetdist = DiscreteFrechetDist(sortrows(firstfront),sortrows(lastfront));
      save([modelname '_results'],'levelset','areasum','isolength','frechetdist');
	end
	for ss_min_thickness = min_thickness_loop;
	   modelname = ['Models/Model_' glacier '_CalvingMinThickness' num2str(ss_min_thickness)];
      disp(['Loading ' modelname]);
      load(modelname);
      levelset = md.results.TransientSolution(end).MaskIceLevelset;
      lsend =double(md.results.TransientSolution(end).MaskIceLevelset<0)*2-1;
      ls1 =double(md.mask.ice_levelset<0)*2-1;
		lsdiff = lsend-ls1;
      areabetween_ind = find(abs(lsdiff)>0);
      areabetween_indices = md.mesh.elements(areabetween_ind,:,:);
      areasum = sum(GetAreas(areabetween_indices,md.mesh.x,md.mesh.y));
      iso=isoline(md, md.mask.ice_levelset,'value',0); %original
      %isolength=sum([1;sqrt(diff(iso(1).x).^2 + diff(iso(1).y).^2)]); %original
      isolength = 0; %Change 5/3/23
      for i=1:length(iso) %Change 5/3/23
			isolength = max(isolength, sum(sqrt(diff(iso(i).x).^2 + diff(iso(i).y).^2)));
		end

      firstfront=isoline(md, md.mask.ice_levelset,'output','matrix');
      lastfront=isoline(md, md.results.TransientSolution(end).MaskIceLevelset,'output','matrix');
		%frechetdist = frechet(sort(rmmissing(firstfront(:,1))),sort(rmmissing(firstfront(:,2))),sort(rmmissing(lastfront(:,1))),sort(rmmissing(lastfront(:,2))));
		frechetdist = DiscreteFrechetDist(sortrows(firstfront),sortrows(lastfront)); 
      save([modelname '_results'],'levelset','areasum','isolength','frechetdist');
	end
	for ss_rc = rc_loop;
		for ss_migmax = migmax_loop;
			modelname = ['Models/Model_' glacier '_CalvingPollard' num2str(ss_rc*100) '_' num2str(ss_migmax)];
			disp(['Loading ' modelname]);
			load(modelname);
			levelset = md.results.TransientSolution(end).MaskIceLevelset;
			lsend =double(md.results.TransientSolution(end).MaskIceLevelset<0)*2-1;
			ls1 =double(md.mask.ice_levelset<0)*2-1;
			lsdiff = lsend-ls1;
			areabetween_ind = find(abs(lsdiff)>0);
			areabetween_indices = md.mesh.elements(areabetween_ind,:,:);
			areasum = sum(GetAreas(areabetween_indices,md.mesh.x,md.mesh.y));
			iso=isoline(md, md.mask.ice_levelset,'value',0); %original
			%isolength=sum([1;sqrt(diff(iso(1).x).^2 + diff(iso(1).y).^2)]); original
			isolength = 0; %Change 5/3/23
			for i=1:length(iso) %Change 5/3/23
				isolength = max(isolength, sum(sqrt(diff(iso(i).x).^2 + diff(iso(i).y).^2)));
			end

			firstfront=isoline(md, md.mask.ice_levelset,'output','matrix');
			lastfront=isoline(md, md.results.TransientSolution(end).MaskIceLevelset,'output','matrix');
			%frechetdist = frechet(sort(rmmissing(firstfront(:,1))),sort(rmmissing(firstfront(:,2))),sort(rmmissing(lastfront(:,1))),sort(rmmissing(lastfront(:,2))));
			if length(lastfront)>0
				frechetdist = DiscreteFrechetDist(sortrows(firstfront),sortrows(lastfront));
				save([modelname '_results'],'levelset','areasum','isolength','frechetdist');
			else
				save([modelname '_results'],'levelset','areasum','isolength');
			end
		end
   end

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

	% This step is for testing purposes only.

	for ss_sigma_max=0.1e6-5*1e4:5e3:0.1e6+5*1e4;
		resultsfile = ['Models/Model_Ross_CalvingVM' num2str(ss_sigma_max) '_results'];
		load(resultsfile);
	end

end%}}}

%%%Steps 16-20
if perform(org,'AreaCalculation'),% {{{
	%This step is not used in the production of the manuscript.

	lsend =double(md.results.TransientSolution(end).MaskIceLevelset<0)*2-1;
	ls1 =double(md.mask.ice_levelset<0)*2-1;
	lsdiff = lsend-ls1;
	areabetween_ind = find(abs(lsdiff)>0);
	areabetween_indices = md.mesh.elements(areabetween_ind,:);
	areasum = sum(GetAreas(areabetween_indices,md.mesh.x,md.mesh.y));
	fprintf('The area between the fronts is %g\n', areasum);
	iso=isoline(md, md.mask.ice_levelset,'value',0);
	isolength=sum([1;sqrt(diff(iso(1).x).^2 + diff(iso(1).y).^2)]);
	fprintf('The length of the first front is %g\n', isolength);
	fprintf('The average distance between the predicted fronts in km is  %g\n', areasum/isolength/1e3);
	plotmodel(md,'data',lsdiff);


	ls1   = md.mask.ice_levelset;
	lsend = md.results.TransientSolution(end).MaskIceLevelset;
	areas = GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
	pos   = find(ls1.*lsend<0);
	areasum = sum(areas(pos));
	fprintf('The area between the fronts is %g km^2\n', areasum/1e6);
	iso=isoline(md, ls1,'value',0);
	isolength=sum([1;sqrt(diff(iso(1).x).^2 + diff(iso(1).y).^2)]);
	fprintf('The length of the first front is %g km\n', isolength/1e3);
	fprintf('The average distance between the predicted fronts in km is  %g\n', areasum/isolength/1e3);
	plotmodel(md,'data',lsdiff);
end%}}}
if perform(org,'Plotting'),% {{{
	
	%This step is superceded by "Revised Plotting"

	calvinglaws = {'Levermann','VM','MinThickness','Pollard'};
	%calvinglaws = {'Levermann'}; %% TEST LINE. COMMENT OUT IF NOT TESTING
	results = 'results';
	directory = './Models';
	filesAndFolders =  dir(directory);
	filesInDir = filesAndFolders(~([filesAndFolders.isdir])); % Returns only the files in the directory

	sgtitle([glacier ' calving law comparison'])
	figure(1)
	for calvinglawnum = 1:length(calvinglaws)
		calvinglaw = calvinglaws{calvinglawnum};

		ContainsGlacier = zeros(numel(filesInDir),1);
		ContainsLaw = zeros(numel(filesInDir),1);
		ContainsResults = zeros(numel(filesInDir),1);
		ContainsAll = zeros(numel(filesInDir),1);
		ContainsLocal = zeros(numel(filesInDir),1);

		for i=1:numel(filesInDir)
			Name = filesInDir(i).name;
			%disp(Name)
			if strcmp(glacier,'Ronne') || strcmp(glacier,'Filchner')
				ContainsGlacier(i) = contains(Name,glacier)-contains(Name,'Filchner-Ronne');
				ContainsGlacier(ContainsGlacier(i)<0)=0;
			else
				ContainsGlacier(i) = contains(Name,glacier);
			end
			ContainsLaw(i) = contains(Name,calvinglaw);
			ContainsResults(i) = contains(Name,results);
			ContainsLocal(i) = contains(Name,'local');
			ContainsMesh(i) = contains(Name,'mesh');
			if (ContainsGlacier(i) && ContainsLaw(i) && ContainsResults(i)  == 1) && ContainsLocal(i) == 0 && ContainsMesh(i) == 0
			%if (ContainsGlacier(i) && ContainsLaw(i) && ContainsResults(i) && ContainsMesh(i)  == 1) && ContainsLocal(i) == 0
				ContainsAll(i) = 1;
			else
				ContainsAll(i) = 0;
			end
		end

		idx = find(ContainsAll); %Indices of filenames that contain both glacier and calving law, and are 'results' files

		resultsMat = zeros(numel(filesInDir),6);

		for k = idx'
			format longG
			idxName = filesInDir(k).name;
			%disp(idxName)
			load([directory '/' idxName])
			numvalue = str2num( regexprep( idxName, {'\D*([\d\.]+\d)[^\d]*', '[^\d\.]*'}, {'$1 ', ' '} ) ); %If not looking at decreased mesh
			%numvalue = str2num( regexprep( idxName(1:end-15), {'\D*([\d\.]+\d)[^\d]*', '[^\d\.]*'}, {'$1 ', ' '} ) );
			resultsMat(k,1) = numvalue(1);
			resultsMat(k,2) = areasum;
			resultsMat(k,3) = isolength;
			resultsMat(k,4) = areasum/isolength/1e3;
			resultsMat(k,5) = frechetdist;
			if calvinglawnum==4
				resultsMat(k,6) = numvalue(2);
			end
		end

		resultsMat = resultsMat(any(resultsMat,2),:); %remove rows with all zeroes
		resultsMat = sortrows(resultsMat,1); %Sort matrix along first column
		mastermax = 160;
		masterMax_frechet = 3e5;

		set(gcf, 'Color', 'w')
		if strcmp(calvinglaw,'VM')
			subplot(length(calvinglaws),1,2);
			plot(resultsMat(:,4),resultsMat(:,1),'-ok', 'MarkerFaceColor','k','MarkerSize',5); xlim([0 mastermax]);
			title('Von Mises');
			ylabel('\sigma_{max} [Pa]')
		elseif strcmp(calvinglaw,'Levermann')
			subplot(length(calvinglaws),1,1);
			plot(resultsMat(:,4),resultsMat(:,1),'-ok', 'MarkerFaceColor','k','MarkerSize',5); xlim([0 mastermax])
			title('Eigencalving');
			ylabel('K [m\timesyr]')
			eigenresults_x = resultsMat(:,4);
			eigenresults_y = resultsMat(:,1);
		elseif strcmp(calvinglaw,'MinThickness')
			subplot(length(calvinglaws),1,3);
			plot(resultsMat(:,4),resultsMat(:,1),'-ok', 'MarkerFaceColor','k','MarkerSize',5); xlim([0 mastermax])
			title('Minimum Thickness');
			ylabel('Min. thickness [m]')
		elseif strcmp(calvinglaw,'Pollard')
			subplot(length(calvinglaws),1,4);
			h=gscatter(resultsMat(:,4),resultsMat(:,1)/100,resultsMat(:,6)/1000); xlim([0 mastermax]);
         [leg,att] = legend('show');
			title(leg, '$\dot{M}_{\mathrm{\cdot max}}$ [km]','Interpreter','latex')
			set(h, 'linestyle', '-');
			title('Crevasse Depth');
         xlabel({'Avg. misfit between observed',' and modeled front [km]'});
         ylabel('r_c')
		end

		%%%%% MAP PLOTS %%%%

		[minvalue,minind] = min(resultsMat(:,4));
		[minvalue_frechet,minind_frechet] = min(resultsMat(:,5));
		if calvinglawnum==4
			min_calving_val =[resultsMat(minind,1),resultsMat(minind,6)];
			min_dist = resultsMat(minind,4);
			min_calving_val_frechet = [resultsMat(minind_frechet,1),resultsMat(minind_frechet,6)];
			min_dist_frechet = resultsMat(minind_frechet,5);
         display(['Minimum distance for ' calvinglaw ' is ' num2str(min_dist)])
			display(['Minimizing metric for ' calvinglaw ' is ' num2str(min_calving_val)])
			display(['Minimum frechet distance for ' calvinglaw ' is ' num2str(min_dist_frechet)])
         display(['Minimizing frechet metric for ' calvinglaw ' is ' num2str(min_calving_val_frechet)])
			md=loadmodel([directory '/Model_' glacier '_Calving' calvinglaw num2str(min_calving_val(1)) '_' num2str(min_calving_val(2)) '.mat']);
		else
			min_calving_val = resultsMat(minind,1); %This is the calving law value that minimizes the difference metric.
			min_dist = resultsMat(minind,4);
			min_calving_val_frechet = [resultsMat(minind_frechet,1),resultsMat(minind_frechet,6)];
         min_dist_frechet = resultsMat(minind_frechet,5);
			display(['Minimum distance for ' calvinglaw ' is ' num2str(min_dist)])
			display(['Minimizing metric for ' calvinglaw ' is ' num2str(min_calving_val)])
			display(['Minimum frechet distance for ' calvinglaw ' is ' num2str(min_dist_frechet)])
         display(['Minimizing frechet metric for ' calvinglaw ' is ' num2str(min_calving_val_frechet)])
			md=loadmodel([directory '/Model_' glacier '_Calving' calvinglaw num2str(min_calving_val) '.mat']);
		end
		md=radarpower(md);
		if calvinglawnum==1
			md1 = md;
		elseif calvinglawnum==2
			md2 = md;
		elseif calvinglawnum==3
			md3 = md;
		elseif calvinglawnum==4
			md4 = md;
		end
	
	end
	
%%%%% Beginning of Mathieu Code %%%%%

	%if ~exist('md1','var')
	%	md1=loadmodel('Models/Model_Ronne_CalvingVM150000.mat'); %for example
	%	md1 = radarpower(md1);
	%
	%	md2=loadmodel('Models/Model_Ronne_CalvingMinThickness240.mat'); %for example
	%	md2 = radarpower(md2);
	%end

	%Process model
	[X Y]=meshgrid(md1.radaroverlay.x,md1.radaroverlay.y);
	vmin=1; vmax = min(3000,max(md1.inversion.vel_obs));
	firstfront=isoline(md1, md1.mask.ice_levelset,'output','matrix');

	%Plot options
	figure(2)
	options = pairoptions('tightsubplot',2,'gap',[.08 .02],'vmargin',[.15 .05],'hmargin',[.05 .05]);
	%options = pairoptions('tightsubplot',1,'vmargin',[.15 .05]);
	left=.2;
	bottom=.3;
	width=0.8;
	height=0.6;
	set(gcf,'color','w','units','normalized','Position',[left bottom width height]);clf;
	cm = getcolormap(pairoptions('colormap',parula()));
	
	for i=1:4
		switch(i)
			case 1; md=md1; titlestr='Eigencalving (EC)';
			case 2; md=md2; titlestr='Von Mises (VM)';
			case 3; md=md3; titlestr='Min. Thickness (MT)';
			case 4; md=md4; titlestr='Crevasse Depth (CD)';
			otherwise error('not supported');
		end
		subplotmodel(2,2,i,options);
		data_in = md.results.TransientSolution(end).Vel; %data_in(data_in==0)=NaN;
		mask_in = md.results.TransientSolution(end).MaskIceLevelset; %data_in(data_in==0)=NaN;
		data = InterpFromMeshToGrid(md.mesh.elements,md.mesh.x,md.mesh.y,[data_in],X(1,1:end),Y(1:end,1),NaN);
		mask = InterpFromMeshToGrid(md.mesh.elements,md.mesh.x,md.mesh.y,[mask_in],X(1,1:end),Y(1:end,1),NaN);
		data(mask>0) = NaN;
		overlay(X(1,:),Y(:,1), data,'caxis',[vmin vmax],'colormap',cm,'log',10,'hemisphere',-1);
		hold on
		lastfront=isoline(md, md.results.TransientSolution(end).MaskIceLevelset,'output','matrix');
		h1=plot(firstfront(:,1),firstfront(:,2),'-k','LineWidth',2.5);
		h2=plot(lastfront(:,1),lastfront(:,2),'-m','LineWidth',2.5);
		title(titlestr,'FontSize',14);

		scaleruler(pairoptions('scaleruler',[min(md.mesh.x)+(max(md.mesh.x)-min(md.mesh.x))/21 max(md.mesh.y)-(max(md.mesh.y)-min(md.mesh.y))/1.05 50e3 (max(md.mesh.y)-min(md.mesh.y))/100 1],'scaleruler_fontsize',10,'FontColor','k'));

	end
	
	legendx = left+0.1;
	legendy = bottom-0.23;
	legendwidth = .15;
	legendheight = .02;

	h=legend([h1 h2],{'Initial front','Front after 200 years'});
   pos = get(h,'position');
   set(h, 'position',[legendx legendy legendwidth legendheight])

	%Add colorbar
	manualcb(vmin,vmax,cm,'Position',[legendx+0.2 legendy .3 .02],'log',10,'orientation','horizontal','title','Velocity (m/yr)');


	hreg = showregion(min(md.mesh.x),min(md.mesh.y),max(md.mesh.x),max(md.mesh.y),[legendx-1.1 legendy-0.4 .3 .3],'color',[1 0 0],'offset2',0*10^3,'LineWidth',1,'hemisphere',-1,'marker',1);
   axes(hreg)


end%}}}
if perform(org,'CalvingPollard_local'),% {{{
   md = loadmodel(org,'TransientPrep');

	%This step is to test the Pollard crevasse depth law for floating ice shelves on the local machine, not to run a paremeter sweep on a supercomputer.
   md.calving=calvingpollard();
   
	md.calving.rc = 0.45;

	%md.transient.requested_outputs={'default','IceVolume','Calvingratex','Calvingratey','CalvingCalvingrate','Gradient1'};
	md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation','Calvingratex','Calvingratey','CalvingCalvingrate'};

   %Add more output
   md.settings.solver_residue_threshold = 1e-5;

   md.levelset.stabilization = 5;
   md.levelset.reinit_frequency=10;
   md.levelset.migration_max = 6e3;
   md.cluster = cluster;

   md=solve(md,'tr');
   savemodel(org,md);
end%}}}
if perform(org,['CalvingPollard' num2str(ss_rc*100) '_' num2str(ss_migmax)]),% {{{
   md = loadmodel(org,'TransientPrep');

   md.calving=calvingpollard();

   md.calving.rc = ss_rc;

   %md.transient.requested_outputs={'default','IceVolume','Calvingratex','Calvingratey','CalvingCalvingrate','Gradient1'};
   md.transient.requested_outputs={'default','IceVolume','IceVolumeAboveFloatation','Calvingratex','Calvingratey','CalvingCalvingrate'};

   %Add more output
   md.settings.solver_residue_threshold = 1e-5;

   md.levelset.stabilization = 5;
   md.levelset.migration_max = ss_migmax;
   md.cluster = cluster;

	%Make sure jobs are submitted without MATLAB waiting for the results, and job name is unique
   md.settings.waitonlock = 0;
   md.cluster.interactive = 0; %only needed if you are using the generic cluster
	md.miscellaneous.name = [glacier '_CD_' num2str(ss_rc*100) '_' num2str(ss_migmax)];

   %Submit job or download results, make sure that there is no runtime name (that includes the date)
   md=solve(md,'Transient','runtimename',false,'loadonly',loadonly);
   %md=solve(md,'tr');
   %savemodel(org,md);

   if loadonly
      savemodel(org,md);
   end

end%}}}
if perform(org,'CalibrationPlotting'),% {{{
   %This step reproduced Figure 1 map plots from the manuscript.
	%Load Ross advance, stable, and retreat
   results = 'results';
   directory = './Models/';
   filesAndFolders =  dir(directory);
   filesInDir = filesAndFolders(~([filesAndFolders.isdir]));
   RossRetreatFile_EC = 'Model_Ross_CalvingLevermann1000000000.mat';
   RossAdvanceFile_EC = 'Model_Ross_CalvingLevermann100000000.mat';
   RossSteadyFile_EC = 'Model_Ross_CalvingLevermann300000000.mat';
   md1 = radarpower(loadmodel([directory RossRetreatFile_EC]));
   md3 = radarpower(loadmodel([directory RossAdvanceFile_EC]));
   md2 = radarpower(loadmodel([directory RossSteadyFile_EC]));
   [X Y]=meshgrid(md1.radaroverlay.x,md1.radaroverlay.y);
   vmin=1; vmax = min(3000,max(md1.inversion.vel_obs));
   firstfront=isoline(md1, md1.mask.ice_levelset,'output','matrix');

	options = pairoptions('tightsubplot',1,'gap',[.01 .05],'vmargin',[.12 .05],'hmargin',[.02 .02]);
   set(gcf,'color','w','Position',[600 200 1100 700]);clf;
   cm = getcolormap(pairoptions('colormap',parula()));

   for i=1:3
      switch(i)
         case 1; md=md1;
         case 2; md=md2;
         case 3; md=md3;
         otherwise error('not supported');
      end
      subplotmodel(3,1,i,options);
      data_in = md.results.TransientSolution(end).Vel; %data_in(data_in==0)=NaN;
      mask_in = md.results.TransientSolution(end).MaskIceLevelset; %data_in(data_in==0)=NaN;
      data = InterpFromMeshToGrid(md.mesh.elements,md.mesh.x,md.mesh.y,[data_in],X(1,1:end),Y(1:end,1),NaN);
      mask = InterpFromMeshToGrid(md.mesh.elements,md.mesh.x,md.mesh.y,[mask_in],X(1,1:end),Y(1:end,1),NaN);
      data(mask>0) = NaN;
      overlay(X(1,:),Y(:,1), data,'caxis',[vmin vmax],'colormap',cm,'log',10,'hemisphere',-1);
      hold on
      lastfront=isoline(md, md.results.TransientSolution(end).MaskIceLevelset,'output','matrix');
      h1=plot(firstfront(:,1),firstfront(:,2),'-k','LineWidth',2.5);
      h2=plot(lastfront(:,1),lastfront(:,2),'-m','LineWidth',2.5);
		scaleruler(pairoptions('scaleruler',[min(md.mesh.x)+(max(md.mesh.x)-min(md.mesh.x))/21 max(md.mesh.y)-(max(md.mesh.y)-min(md.mesh.y))/1.05 200e3 (max(md.mesh.y)-min(md.mesh.y))/100 1],'scaleruler_fontsize',10,'FontColor','k'));
	end
	hold off
	legendx=0.65; legendy=0.4; legendwidth=.1; legendheight=.05;
	legend('Initial front','Front after 200 years','FontSize',8,'Position',[legendx legendy legendwidth legendheight]);
	manualcb(vmin,vmax,cm,'Position',[legendx+legendwidth/3 legendy+1.5*legendheight legendwidth/6 legendheight*7],'log',10,'orientation','vertical','title','Velocity (m/yr)');
   hreg = showregion(min(md.mesh.x),min(md.mesh.y),max(md.mesh.x),max(md.mesh.y),[legendx/1.3 legendy*1.4 .4 .4],'color',[1 0 0],'offset2',0*10^3,'LineWidth',1,'hemisphere',-1,'marker',1);
   axes(hreg);	


end%}}}

%%%Steps 21+
if perform(org,'RevisedPlotting'),% {{{
	%This step reproduces the calving law comparison figures, with multiple fronts superimposed on the same figure.	
   calvinglaws = {'Levermann','VM','MinThickness','Pollard'};
   %calvinglaws = {'Levermann'}; %% TEST LINE. COMMENT OUT IF NOT TESTING
   results = 'results';
   directory = './Models';
   filesAndFolders =  dir(directory);
   filesInDir = filesAndFolders(~([filesAndFolders.isdir])); % Returns only the files in the directory

   sgtitle([glacier ' calving law comparison'])
   figure(1)
   for calvinglawnum = 1:length(calvinglaws)
      calvinglaw = calvinglaws{calvinglawnum};

      ContainsGlacier = zeros(numel(filesInDir),1);
      ContainsLaw = zeros(numel(filesInDir),1);
      ContainsResults = zeros(numel(filesInDir),1);
      ContainsAll = zeros(numel(filesInDir),1);
      ContainsLocal = zeros(numel(filesInDir),1);

      for i=1:numel(filesInDir)
         Name = filesInDir(i).name;
         %disp(Name)
         if strcmp(glacier,'Ronne') || strcmp(glacier,'Filchner')
            ContainsGlacier(i) = contains(Name,glacier)-contains(Name,'Filchner-Ronne');
            ContainsGlacier(ContainsGlacier(i)<0)=0;
         else
            ContainsGlacier(i) = contains(Name,glacier);
         end
         ContainsLaw(i) = contains(Name,calvinglaw);
         ContainsResults(i) = contains(Name,results);
         ContainsLocal(i) = contains(Name,'local');
         ContainsMesh(i) = contains(Name,'mesh');
         if (ContainsGlacier(i) && ContainsLaw(i) && ContainsResults(i)  == 1) && ContainsLocal(i) == 0 && ContainsMesh(i) == 0
            ContainsAll(i) = 1;
         else
            ContainsAll(i) = 0;
         end
      end

      idx = find(ContainsAll); %Indices of filenames that contain both glacier and calving law, and are 'results' files

      resultsMat = zeros(numel(filesInDir),6);

      for k = idx'
         format longG
         idxName = filesInDir(k).name;
         %disp(idxName)
         load([directory '/' idxName])
         numvalue = str2num( regexprep( idxName, {'\D*([\d\.]+\d)[^\d]*', '[^\d\.]*'}, {'$1 ', ' '} ) ); %If not looking at decreased mesh
         %numvalue = str2num( regexprep( idxName(1:end-15), {'\D*([\d\.]+\d)[^\d]*', '[^\d\.]*'}, {'$1 ', ' '} ) );
         resultsMat(k,1) = numvalue(1);
         resultsMat(k,2) = areasum;
         resultsMat(k,3) = isolength;
         resultsMat(k,4) = areasum/isolength/1e3;
         resultsMat(k,5) = frechetdist;
         if calvinglawnum==4
            resultsMat(k,6) = numvalue(2);
         end
      end

      resultsMat = resultsMat(any(resultsMat,2),:); %remove rows with all zeroes
      resultsMat = sortrows(resultsMat,1); %Sort matrix along first column
      mastermax = 160;
      masterMax_frechet = 3e5;

      set(gcf, 'Color', 'w')
      if strcmp(calvinglaw,'VM')
         subplot(length(calvinglaws),1,2);
         plot(resultsMat(:,4),resultsMat(:,1),'-ok', 'MarkerFaceColor','k','MarkerSize',5); xlim([0 mastermax]);
         title('Von Mises');
         %xlabel({'Avg. misfit between observed',' and modeled front [km]'});
         ylabel('\sigma_{max} [Pa]')
      elseif strcmp(calvinglaw,'Levermann')
         subplot(length(calvinglaws),1,1);
         plot(resultsMat(:,4),resultsMat(:,1),'-ok', 'MarkerFaceColor','k','MarkerSize',5); xlim([0 mastermax])
         title('Eigencalving');
         %xlabel({'Avg. misfit between observed',' and modeled front [km]'});
         ylabel('K [m\timesyr]')
         eigenresults_x = resultsMat(:,4);
         eigenresults_y = resultsMat(:,1);
      elseif strcmp(calvinglaw,'MinThickness')
         subplot(length(calvinglaws),1,3);
         plot(resultsMat(:,4),resultsMat(:,1),'-ok', 'MarkerFaceColor','k','MarkerSize',5); xlim([0 mastermax])
         title('Minimum Thickness');
         %xlabel({'Avg. misfit between observed',' and modeled front [km]'});
         ylabel('Min. thickness [m]')
      elseif strcmp(calvinglaw,'Pollard')
         subplot(length(calvinglaws),1,4);
         h=gscatter(resultsMat(:,4),resultsMat(:,1)/100,resultsMat(:,6)/1000); xlim([0 mastermax]);
         [leg,att] = legend('show');
         title(leg, '$\dot{M}_{\mathrm{\cdot max}}$ [km]','Interpreter','latex')
         set(h, 'linestyle', '-');
         title('Crevasse Depth');
         xlabel({'Avg. misfit between observed',' and modeled front [km]'});
         ylabel('r_c')
      end

      %%%%% MAP PLOTS %%%%
		[minvalue,minind] = min(resultsMat(:,4));
      [minvalue_frechet,minind_frechet] = min(resultsMat(:,5));
      if calvinglawnum==4
         min_calving_val =[resultsMat(minind,1),resultsMat(minind,6)];
         min_dist = resultsMat(minind,4);
         min_calving_val_frechet = [resultsMat(minind_frechet,1),resultsMat(minind_frechet,6)];
         min_dist_frechet = resultsMat(minind_frechet,5);
         display(['Minimum distance for ' calvinglaw ' is ' num2str(min_dist)])
         display(['Minimizing metric for ' calvinglaw ' is ' num2str(min_calving_val)])
         display(['Minimum frechet distance for ' calvinglaw ' is ' num2str(min_dist_frechet)])
         display(['Minimizing frechet metric for ' calvinglaw ' is ' num2str(min_calving_val_frechet)])
         md=loadmodel([directory '/Model_' glacier '_Calving' calvinglaw num2str(min_calving_val(1)) '_' num2str(min_calving_val(2)) '.mat']);
      else
         min_calving_val = resultsMat(minind,1); %This is the calving law value that minimizes the difference metric.
         min_dist = resultsMat(minind,4);
         min_calving_val_frechet = resultsMat(minind_frechet,1);
         min_dist_frechet = resultsMat(minind_frechet,5);
         display(['Minimum distance for ' calvinglaw ' is ' num2str(min_dist)])
         display(['Minimizing metric for ' calvinglaw ' is ' num2str(min_calving_val)])
         display(['Minimum frechet distance for ' calvinglaw ' is ' num2str(min_dist_frechet)])
         display(['Minimizing frechet metric for ' calvinglaw ' is ' num2str(min_calving_val_frechet)])
         md=loadmodel([directory '/Model_' glacier '_Calving' calvinglaw num2str(min_calving_val) '.mat']);
      end
      md=radarpower(md);
      if calvinglawnum==1
         md1 = md;
      elseif calvinglawnum==2
         md2 = md;
      elseif calvinglawnum==3
         md3 = md;
      elseif calvinglawnum==4
         md4 = md;
      end

   end	

	%Process model
   [X Y]=meshgrid(md1.radaroverlay.x,md1.radaroverlay.y);
   vmin=1; 
	vmax=1000;
	%vmax = min(3000,max(md1.inversion.vel_obs));
   firstfront=isoline(md1, md1.mask.ice_levelset,'output','matrix');
	alpha = .55;

	%Plot options
   f=figure(2);

	%Determine domain outline percentages by ice shelf. Hardcoded
	switch glacier
      case 'Amery'
			percent_x_min = 0; percent_x_max = 80; percent_y_min = 0; percent_y_max = 100;
      case 'Denman'
			 percent_x_min = 0; percent_x_max = 90; percent_y_min = 0; percent_y_max = 100;
      case 'Filchner'
			percent_x_min = 20; percent_x_max = 100; percent_y_min = 20; percent_y_max = 80;
      case 'LarsenC'
			percent_x_min = 0; percent_x_max = 100; percent_y_min = 0; percent_y_max = 100;
      case 'PineIsland'
			percent_x_min = 0; percent_x_max = 68; percent_y_min = 0; percent_y_max = 60;
      case 'Ronne'
			percent_x_min = 0; percent_x_max = 100; percent_y_min = 0; percent_y_max = 85;
      case 'Ross'
			 percent_x_min = 0; percent_x_max = 100; percent_y_min = 0; percent_y_max = 100;
      case 'Shackleton'
			percent_x_min = 0; percent_x_max = 100; percent_y_min = 0; percent_y_max = 100;
      case 'Thwaites'
			percent_x_min = 25; percent_x_max = 65; percent_y_min = 0; percent_y_max = 100;
		case 'Totten'
			percent_x_min = 0; percent_x_max = 50; percent_y_min = 16; percent_y_max = 100;
      otherwise
         error('Not supported')
   end

   options = pairoptions('tightsubplot',2,'gap',[.08 .02],'vmargin',[.15 .05],'hmargin',[.05 .05]);
   %options = pairoptions('tightsubplot',1,'vmargin',[.15 .05]);
   %set(gcf,'color','w','Position',[600 200 1100 700]);clf;
   left=.2;
   bottom=.3;
   width=0.8;
   height=0.6;
   set(gcf,'color','w','units','normalized','Position',[left bottom width height]);clf;
   %set(gcf,'color','w');clf;
   cm = getcolormap(pairoptions('colormap',parula(20)));
	h1=plot(firstfront(:,1),firstfront(:,2),'-k','LineWidth',2.5);
	colorpalette=lines(5); colorpalette(2,:) = [1 0 0];
	hold on
	for i=1:5
		color = colorpalette(i,:);
      switch(i)
         case 1; md=md1; titlestr='Eigencalving (EC)'; 
         case 2; md=md2; titlestr='Von Mises (VM)'; 
         case 3; md=md3; titlestr='Min. Thickness (MT)';
         case 4; md=md4; titlestr='Crevasse Depth (CD)';
         case 5; md=md1; titlstr='Initial'; color = 'w';
			otherwise error('not supported');
      end
      %subplotmodel(2,2,i,options);
		if i==5
			data_in=md.results.TransientSolution(1).Vel; 
			mask_in = md.results.TransientSolution(1).MaskIceLevelset;
		else
			data_in = md.results.TransientSolution(end).Vel; %data_in(data_in==0)=NaN;
			mask_in = md.results.TransientSolution(end).MaskIceLevelset; %data_in(data_in==0)=NaN;
		end
		x_mesh_min = (percent_x_min * (max(md.mesh.x) - min(md.mesh.x)) / 100) + min(md.mesh.x);
		x_mesh_max = (percent_x_max * (max(md.mesh.x) - min(md.mesh.x)) / 100) + min(md.mesh.x);
		idx_x_max = md.mesh.x(:)<=x_mesh_max;
		idx_x_min = md.mesh.x(:)>=x_mesh_min;
		idx_X_min = X(1,:)>=x_mesh_min;
		idx_X_max = X(1,:)<=x_mesh_max;
		X_min = find(idx_X_min,1,'first');
		X_max = find(idx_X_max,1,'last');
		mesh_x_subset = intersect(md.mesh.x(idx_x_max),md.mesh.x(idx_x_min),'stable');
		data_in_subset = intersect(data_in(idx_x_max),data_in(idx_x_min),'stable');
      y_mesh_min = (percent_y_min * (max(md.mesh.y) - min(md.mesh.y)) / 100) + min(md.mesh.y);
      y_mesh_max = (percent_y_max * (max(md.mesh.y) - min(md.mesh.y)) / 100) + min(md.mesh.y);
		idx_y_max = md.mesh.y(:)<=y_mesh_max;
      idx_y_min = md.mesh.y(:)>=y_mesh_min;
		idx_Y_min = Y(:,1)>=y_mesh_min;
      idx_Y_max = Y(:,1)<=y_mesh_max;
      Y_min = find(idx_Y_min,1,'first');
      Y_max = find(idx_Y_max,1,'last');
      mesh_y_subset = intersect(md.mesh.y(idx_y_max),md.mesh.y(idx_y_min),'stable');
		data = InterpFromMeshToGrid(md.mesh.elements,md.mesh.x,md.mesh.y,[data_in],X(1,1:end),Y(1:end,1),NaN);
		mask = InterpFromMeshToGrid(md.mesh.elements,md.mesh.x,md.mesh.y,[mask_in],X(1,1:end),Y(1:end,1),NaN);
		
		
		data(mask>0) = NaN;
      if i==5
			%overlay(X(1,:),Y(:,1), data,'caxis',[vmin vmax],'colormap',cm,'log',10,'hemisphere',-1,'alpha',alpha);
			overlay(X(1,X_min:X_max),Y(Y_min:Y_max,1), data(Y_min:Y_max,X_min:X_max),'caxis',[vmin vmax],'colormap',cm,'log',10,'hemisphere',-1,'alpha',alpha,'frame',5)
		else
			lastfront=isoline(md, md.results.TransientSolution(end).MaskIceLevelset,'output','matrix');
			h2=plot(lastfront(:,1),lastfront(:,2),'-','Color', color,'LineWidth',2.5);
		end
   end
	switch glacier
		case {'Amery', 'PineIsland','LarsenC'}
			textpos='a';
		case {'Ronne', 'Thwaites', 'Shackleton'}
			textpos='b';
		case {'Filchner','Ross','Denman'}
			textpos='c';
		case {'Totten'}
			textpos='d';
		otherwise
			error('Not supported');
	end
	
	%Formatting values switch. Hardcoded figure format values.
	switch glacier
      case 'Amery'
			cbarlength = 0.2; cbarx = 0.4; cbary = 0.05; northx=0.02; northy=0.1; northlengthdiv=15; hregx = 0.85; hregy = 0.01; sr_length=100e3;
		case 'Denman'
			cbarlength = 0.2; cbarx = 0.4; cbary = 0.05; northx=percent_x_min+0.034; northy=percent_y_min+0.15; northlengthdiv=15; hregx = percent_x_max/100-0.05; hregy = 0.025; sr_length=20e3;
		case 'Filchner'
			cbarlength = 0.2; cbarx = 0.4; cbary = 0.05; northx=0.08; northy=0.075; northlengthdiv=15; hregx = 0.63; hregy = 0.01; sr_length=200e3;
		case 'LarsenC'
			cbarlength = 0.2; cbarx = 0.4; cbary = 0.05; northx=0.2; northy=0.1; northlengthdiv=8; hregx = 0.6; hregy = 0.01; sr_length=100e3;
		case 'PineIsland'
			cbarlength = 0.2; cbarx = 0.4; cbary = 0.05; northx=percent_x_min+0.14; northy=percent_y_min+0.1; northlengthdiv=15; hregx = percent_x_max/100-0.04; hregy = 0.01; sr_length=60e3;
		case 'Ronne'
			cbarlength = 0.2; cbarx = 0.4; cbary = 0.05; northx=0.17; northy=0.1; northlengthdiv=15; hregx = 0.68; hregy = 0.01; sr_length=200e3;
		case 'Ross'
			cbarlength = 0.2; cbarx = 0.4; cbary = 0.05; northx=0.14; northy=0.20; northlengthdiv=15; hregx = 0.635; hregy = 0.01; sr_length=200e3;
		case 'Shackleton'
         cbarlength = 0.2; cbarx = 0.4; cbary = 0.05; northx=percent_x_min+0.6; northy=percent_y_min+0.15; northlengthdiv=15; hregx = percent_x_max/100-0.3; hregy = 0.01; sr_length=60e3;
		case 'Thwaites'
			cbarlength = .2; cbarx = 0.4; cbary = 0.05; northx=0.08; northy=0.15; northlengthdiv=15; hregx = 0.55; hregy = 0.01; sr_length=40e3;
		case 'Totten'
			cbarlength = 0.2; cbarx = 0.4; cbary = 0.05; northx=0.03; northy=0.12; northlengthdiv=15; hregx = percent_x_max/100-0.01; hregy = 0.01; sr_length=40e3;
		otherwise
			error('Not supported')
	end
	if vmax<=1000
		vmax=1000; %Ensure that max velocity is 1000 m/yr
	else
		vmax=vmax;
	end
	text(0.01,0.9999,textpos,'units','normalized','FontSize',28,'Color','k','HorizontalAlignment','left','VerticalAlignment','top');
	hold off
	switch glacier
		case 'Shackleton'
			sbar_offset = 0.55;
		otherwise
			sbar_offset = 0.02;
	end
	scaleruler(pairoptions('scaleruler',[min(md.mesh.x)+(max(md.mesh.x)-min(md.mesh.x))*(percent_x_min/100+sbar_offset) min(md.mesh.y)+(max(md.mesh.y)-min(md.mesh.y))*(percent_y_min/100+0.02) sr_length (max(md.mesh.y)-min(md.mesh.y))/100 1],'scaleruler_fontsize',16,'FontColor','k'));
	north_x0 = min(md.mesh.x)+(max(md.mesh.x)-min(md.mesh.x))*(percent_x_min/100+northx);
   north_y0 = min(md.mesh.y)+(max(md.mesh.y)-min(md.mesh.y))*(percent_y_min/100+northy);
   north_length = (max(md.mesh.x)-min(md.mesh.x))/northlengthdiv;
	north_ratio = 0.5;
   north_width = north_length/10;
   north_fontsize = 16;
   northarrow([north_x0,north_y0,north_length,north_ratio,north_width,north_fontsize]);
	legend('Initial','EC','VM','MT','CD')
	
	%Add colorbar
	u1 = findall(f,'type','axes');
   u2 = manualcb(vmin,vmax,(1-alpha)*cm+alpha,'Position',[cbarx cbary cbarlength .02],'log',10,'orientation','horizontal','title','Velocity (m yr^{-1})');	
   hreg = showregion(min(md.mesh.x),min(md.mesh.y),max(md.mesh.x),max(md.mesh.y),[hregx hregy .14 .14],'color',[1 0 0],'offset2',0*10^3,'LineWidth',1,'hemisphere',-1,'marker',1);
   axes(hreg);
	align([u1 u2],'center','distribute');
end%}}}

return

%The following code is scratch paper for test plotting purposes.

%plot first and last ice front
firstfront=isoline(md, md.mask.ice_levelset,'output','matrix');
lastfront=isoline(md, md.results.TransientSolution(end).MaskIceLevelset,'output','matrix');
plotmodel(md,'data',md.results.TransientSolution(end).Vel,'gridded',1,'levelset',md.results.TransientSolution(end).MaskIceLevelset); hold on
h1=plot(firstfront(:,1),firstfront(:,2),'-k','LineWidth',2.5);
h2=plot(lastfront(:,1),lastfront(:,2),'-c','LineWidth',2.5);
legend([h1 h2],{'Initial front','Front after 200 years'});

%plot front through time
N = numel(md.results.TransientSolution);
cmap = parula(N);
plotmodel(md,'data',md.results.TransientSolution(end).Vel,'gridded',1,'levelset',md.results.TransientSolution(end).MaskIceLevelset); hold on
for i=1:N
	front=isoline(md, md.results.TransientSolution(i).MaskIceLevelset,'output','matrix');
	plot(front(:,1),front(:,2),'-','Color',cmap(i,:));
end

%Return average calving rate (not truly a calving rate)
levelset = md.results.TransientSolution(end).MaskIceLevelset;
lsend =double(md.results.TransientSolution(end).MaskIceLevelset<0)*2-1;
ls1 =double(md.mask.ice_levelset<0)*2-1;
lsdiff = lsend-ls1;
areabetween_ind = find(abs(lsdiff)>0);
areabetween_indices = md.mesh.elements(areabetween_ind,:,:);
areasum = sum(GetAreas(areabetween_indices,md.mesh.x,md.mesh.y));
iso=isoline(md, md.mask.ice_levelset,'value',0);
isolength=sum([1;sqrt(diff(iso(1).x).^2 + diff(iso(1).y).^2)]);
firstfront=isoline(md, md.mask.ice_levelset,'output','matrix');
lastfront=isoline(md, md.results.TransientSolution(end).MaskIceLevelset,'output','matrix');
misfit=areasum/isolength/1e3;
calvingrate_avg = misfit/md.timestepping.final_time; %[km/yr]

%plot strain rate
[dvxdx,dvxdy]=slope(md,md.results.TransientSolution(end).Vx);
[dvydx,dvydy]=slope(md,md.results.TransientSolution(end).Vy);
plotmodel(md,'data',dvxdx,'data',dvxdy,'data',dvydx,'data',dvydy,'title','dvx/dx','fontsize',10,'title','dvx/dy','fontsize',10,'title','dvy/dx','fontsize',10,'title','dvy/dy','fontsize',10,'caxis',[-0.01 0.01],'caxis',[-0.01 0.01],'caxis',[-0.01 0.01],'caxis',[-0.01 0.01])
%Backup work
!rclone -v --skip-links --drive-upload-cutoff=16384k sync /totten_1/jwilner/proj-wilner/Calving gdrive_joel:proj-wilner/Calving

%Functions
function fric=EstimateFric_Budd(md) % {{{

   disp('Estimating friction coefficient - Budd')
   % initial guess from driving stress
   min_velocity=0.1; % m/yr
   min_pressure=1.;% Pa
   min_c=1;

   asurf=averaging(md,md.geometry.surface,20); % maybe executing 20 L2 projection is ok
   [sx,sy,s]=slope(md,asurf); % slope 's' comes on elements
   sslope=averaging(md,s,1); % average the slope once on the vertices, because 's' comes on elements, we need this data on vertices

   vel=md.inversion.vel_obs;
   vel=max(vel,min_velocity); % setting minimum velocity value
   vel=vel/md.constants.yts; % S.I.

   Neff=(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.base)*md.constants.g;
   Neff(find(Neff<=0))=min_pressure; % setting minimum positve pressure

   driving_stress=md.materials.rho_ice*md.constants.g*md.geometry.thickness.*(sslope);
   c=sqrt(driving_stress./(Neff.*vel));
   c=max(c,min_c);

   fric=c; % initial guess

end%}}}
