% MakeISSMMesh_v1

% Mike Wolovick, 8/2/2021

% This script is my own wrapper for the ISSM commands for making a refined
% model mesh.  It is based on runme.m from the PIG example.

% v2:  uses MultiWavelengthInterpolator instead of slowly doing all of that
% in here.  Also, extends ice thickness beyone the calving front beforehand
% in order to prevent artifacts with too-small ice thickness near the
% front.

% v3:  adds the option of expanding the domain outline into the ocean.  If
% so, no longer assumes that one boundary is the ice front.  Uses laplacian
% interpolation to expand ice thickness into the ocean and to fill velocity
% holes.  In addition, the construction of the control field is different.
% It proceeds by adding the log ov various fields in stages, so it no
% longer needs a text key at the start.  

% Each component of the control field is added like so:
% Component=sqrt(Component.^2+componentfloor^2);
% ControlField=ControlField+componentweight*log10(Component/componentfloor);

% The exception is the mask, which does not require a stabilizing floor or
% taking the logarithm.  

% v4:  allows you to separately set max mesh sizes in floating,
% grounded, and open ocean regions.  In addition, the analysis figures
% include area-weighted histograms in addition to basic histograms.  This
% script also allows you to loop through a series with variable err, hmin,
% and hmax, but the same controlfield.

clearvars
close all
tic

%% Parameters:

% File names and paths:
boundaryfile='~/Dropbox/FRISModeling/Boundaries/BoundaryPickingTest_v3_50km.exp';
gridinputfile='~/Dropbox/CombinedGrids/AntarcticaCombinedGrids_v4.mat';
velinputfile='~/Dropbox/CombinedGrids/AntarcticaMergedVel_MEASURESv2withphase_plus_FRmosaic.mat';  % can be commented
issmscriptpath='~/Dropbox/ISSM/';
generaloutputfile='~/Dropbox/FRISModeling/Meshes/FRISDomainMesh_SizeSeries_v1';  % no trailing .mat

% Mesh parameters:
hinit=2e3;                             % m   element size for the initial mesh
allhmax=[(1e4:1e4:1e5)',(5e3:5e3:5e4)',1e5*ones(10,1)];   % [Nx3] m   maximum element size of the final mesh [grounded,floating,ocean]
allhmin=(1e2:1e2:1e3)';                % [Nx1] m   minimum element size of the final mesh 
gradation=1.5;                         % unitless>1   maximum size ratio between two neighboring elements
allerr=(.01:.01:.1)';                  % [Nx1] controlfield units.   maximum error between interpolated and control field
numiterations_mesh=5;                  % integer

% Distance to expand into ocean:
dist2expand=5e5;                       % m (0 deactivates)
expandangle=105;                       % degrees counterclockwise relative to mean calving front orientation

% Velocity processing parameters:
velwavelength=2e3;                     % m (0 deactivates smoothing)
mingoodvalue4velfill=2;                % m/yr

% Control field parameters:
velweight=0;                           % unitless
velfloor=1;                            % m/yr
fluxweight=0;                          % unitless
fluxfloor=1e4;                         % m^2/yr
dissweight=1;                          % unitless
dissfloor=2e-2;                        % W/m^2
strainweight=4;                        % unitless (note: strain rate only applied to floating shelf)
strainfloor=3e-3;                      % 1/yr
maskweight=1;                          % unitless
controlfieldwavelength=2e3;            % m (0 deactivates smoothing)
frontmaskshift=-.75;                   % unitless (everything other than open ocean shifted by this amount)

% Include firn correction?
dofirncorrection=0;                    % logical

% Densities:
rho_i=917;                             % kg/m^3
rho_sw=1027;                           % kg/m^3

% Show progress in mesh construction?
showintermediateiterations=0;          % logical

% Seconds per year:
secondsperyear=60*60*24*365.25;        % s/yr
	
% Map cropping limits:
xlims_crop=[-1800,1000]*1e3;           % [1x2] m
ylims_crop=[-300,1750]*1e3;            % [1x2] m

% Threshold for analyzing fast-flow areas:
velanalysisthrehsold=10;               % m/yr (smoothing weighting btwn 1/2 and 2x this value)

% Axes limits for figures:
histogramlengthlims=[1e-1,1e2];        % [1x2] km
histogramarealims=[1e-2,1e4];          % [1x2] km^2
histogramanilims=[0,6];                % [1x2] unitless
meshsizeclims=[1e-2,1e2];              % [1x2] km

% Output size and resolution:
pagesize1=[12,10];                     % [1x2] inches (single map figures)
pagesize2=[20,10];                     % [1x2] inches (2x3 histogram figures)
pagesize3=[16,16];                     % [1x2] inches (3x2 map figures)
resolution=300;                        % dpi

% Record parameters:
MeshingParameters=struct('hinit',hinit,'hmax',allhmax,'hmin',allhmin,'gradation',gradation,'err',allerr,'numiterations_mesh',numiterations_mesh,'dist2expand',dist2expand,'expandangle',expandangle,...
    'velwavelength',velwavelength,'mingoodvalue4velfill',mingoodvalue4velfill,...
    'velweight',velweight,'velfloor',velfloor,'fluxweight',fluxweight,'fluxfloor',fluxfloor,'dissweight',dissweight,'dissfloor',dissfloor,...
    'strainweight',strainweight,'strainfloor',strainfloor,'maskweight',maskweight,'controlfieldwavelength',controlfieldwavelength,'frontmaskshift',frontmaskshift,...
    'dofirncorrection',dofirncorrection,'rho_i',rho_i,'rho_sw',rho_sw,...
    'xlims_crop',xlims_crop,'ylims_crop',ylims_crop,'boundaryfile',boundaryfile,'gridinputfile',gridinputfile);
if exist('velinputfile','var')
    MeshingParameters.velinputfile=velinputfile;
end

%% Load Grud Inputs:

% Check number of samples in series:
numsamples=length(allerr);
if size(allhmax,1)~=numsamples
    error('Parameter "allhmax" must have the same number of rows as the length of parameter "allerr".')
end
if length(allhmin)~=numsamples
    error('Parameter "allhmin" must have the same length as parameter "allerr".')
end

% Communicate:
disp('Loading inputs...')

% Add the ISSM script path:
addpath(genpath(issmscriptpath))

% Load gridded input file:
load(gridinputfile,'X','Y','Mask','IceThick','BedElev','SurfElev')
if dissweight~=0
    load(gridinputfile,"DrivingStress_alongflow")
end
if ~exist('velinputfile','var')
    load(gridinputfile,'U_obs','V_obs','VelMag_obs')
end
if dofirncorrection
    load(gridinputfile,'FirnCorrection')
end

% Crop input grids:
ind1=find(Y>=ylims_crop(1),1,'first');
ind2=find(Y<=ylims_crop(2),1,'last');
ind3=find(X>=xlims_crop(1),1,'first');
ind4=find(X<=xlims_crop(2),1,'last');
X_grid=X(ind3:ind4);
Y_grid=Y(ind1:ind2);
Mask_grid=Mask(ind1:ind2,ind3:ind4);
BedElev_grid=BedElev(ind1:ind2,ind3:ind4);
SurfElev_grid=SurfElev(ind1:ind2,ind3:ind4);
IceThick_grid=IceThick(ind1:ind2,ind3:ind4);
if ~exist('velinputfile','var')
    U_obs_grid=U_obs(ind1:ind2,ind3:ind4);
    V_obs_grid=V_obs(ind1:ind2,ind3:ind4);
    VelMag_obs_grid=VelMag_obs(ind1:ind2,ind3:ind4);
end
if dofirncorrection
    Firn_grid=FirnCorrection(ind1:ind2,ind3:ind4);
end
if dissweight~=0
    DrivingStress_grid=DrivingStress_alongflow(ind1:ind2,ind3:ind4);
end
clear X Y Mask U_obs V_obs VelMag_obs FirnCorrection BedElev IceThick SurfElev DrivingStress_alongflow

% Evaluate grid size:
[ysize,xsize]=size(Mask_grid);
dx_grid=X_grid(2)-X_grid(1);

% Reorder the mask:
ReorderedMask_grid=Mask_grid;  % 0 (ocean) and 2 (grounded ice) remain the same
ReorderedMask_grid(Mask_grid==3)=1; % shelf
ReorderedMask_grid(Mask_grid==1)=3; % rock

% Load separate velocity input file:
if exist('velinputfile','var')
    % Load:
    load(velinputfile,'X','Y','VelMag_merged','U_merged','V_merged')
    % Crop:
    ind1=find(Y>=ylims_crop(1),1,'first');
    ind2=find(Y<=ylims_crop(2),1,'last');
    ind3=find(X>=xlims_crop(1),1,'first');
    ind4=find(X<=xlims_crop(2),1,'last');
    X=X(ind3:ind4);
    Y=Y(ind1:ind2);
    U_merged=U_merged(ind1:ind2,ind3:ind4);
    V_merged=V_merged(ind1:ind2,ind3:ind4);
    VelMag_merged=VelMag_merged(ind1:ind2,ind3:ind4);
    % Smooth before interpolation:
    dx_vel=X(2)-X(1);
    if dx_grid/dx_vel>sqrt(2)
        U_merged=ThreeByThree(U_merged,round((.5*dx_grid/dx_vel)^2));
        V_merged=ThreeByThree(V_merged,round((.5*dx_grid/dx_vel)^2));
        VelMag_merged=ThreeByThree(VelMag_merged,round((.5*dx_grid/dx_vel)^2));
    end
    % Interpolate:
    U_obs_grid=interp2(X,Y,U_merged,repmat(X_grid,[ysize,1]),repmat(Y_grid,[1,xsize]));
    V_obs_grid=interp2(X,Y,V_merged,repmat(X_grid,[ysize,1]),repmat(Y_grid,[1,xsize]));
    VelMag_obs_grid=interp2(X,Y,VelMag_merged,repmat(X_grid,[ysize,1]),repmat(Y_grid,[1,xsize]));
    % Clear extraneous variables:
    clear X Y U_merged V_merged VelMag_merged
end


%% Velocity Processing:

% Communicate:
disp('Processing velocity field...')

% Smoothly fill in holes in observed velocity with neighboring values:
% identify bad values:
BadValues_grid=isnan(VelMag_obs_grid)|VelMag_obs_grid<=mingoodvalue4velfill;
% expand the bad values region by 2 in all directions:
BadValues_grid(2:end-1,2:end-1)=BadValues_grid(1:end-2,1:end-2)|BadValues_grid(2:end-1,1:end-2)|BadValues_grid(3:end,1:end-2)|...
    BadValues_grid(1:end-2,2:end-1)|BadValues_grid(2:end-1,2:end-1)|BadValues_grid(3:end,2:end-1)|...
    BadValues_grid(1:end-2,3:end)|BadValues_grid(2:end-1,3:end)|BadValues_grid(3:end,3:end);
BadValues_grid(2:end-1,2:end-1)=BadValues_grid(1:end-2,1:end-2)|BadValues_grid(2:end-1,1:end-2)|BadValues_grid(3:end,1:end-2)|...
    BadValues_grid(1:end-2,2:end-1)|BadValues_grid(2:end-1,2:end-1)|BadValues_grid(3:end,2:end-1)|...
    BadValues_grid(1:end-2,3:end)|BadValues_grid(2:end-1,3:end)|BadValues_grid(3:end,3:end);
% Flag bad values and interpolate:
U_obs_grid(BadValues_grid)=NaN;
U_obs_grid=LaplacianInterpolation(U_obs_grid,9);
V_obs_grid(BadValues_grid)=NaN;
V_obs_grid=LaplacianInterpolation(V_obs_grid,9);
VelMag_obs_grid(BadValues_grid)=NaN;
VelMag_obs_grid=LaplacianInterpolation(VelMag_obs_grid,9);

% Smooth out velocity:
if velwavelength>0
    U_obs_grid=ThreeByThree(U_obs_grid,round((.5*velwavelength/dx_grid)^2));
    V_obs_grid=ThreeByThree(V_obs_grid,round((.5*velwavelength/dx_grid)^2));
    VelMag_obs_grid=ThreeByThree(VelMag_obs_grid,round((.5*velwavelength/dx_grid)^2));
end

%% Generate Control Field:

% Communicate:
disp('Generating control field...')

% Initialize control field:
ControlField_grid=zeros(ysize,xsize);

% Add flux:
if fluxweight~=0
    % Compute component:
    Flux_grid=VelMag_obs_grid.*IceThick_grid;
    % Apply floor:
    Flux_grid=sqrt(Flux_grid.^2+fluxfloor^2);
    % Apply mask constraints:
    Flux_grid(ReorderedMask_grid<.5|ReorderedMask_grid>2.5)=0; % 0 outside the ice sheet or shelf
    % Add to control field:
    ControlField_grid=ControlField_grid+fluxweight*log10(Flux_grid/fluxfloor);
end

% Add dissipation:
if dissweight~=0
    % Compute component:
    Dissipation_grid=VelMag_obs_grid.*DrivingStress_grid/secondsperyear;
    % Apply mask constraints:
    Dissipation_grid(ReorderedMask_grid<1.5|ReorderedMask_grid>2.5)=0; % zero outside of grounded ice sheet
    % Apply floor:
    Dissipation_grid=sqrt(Dissipation_grid.^2+dissfloor^2);
    % Add to control field:
    ControlField_grid=ControlField_grid+dissweight*log10(Dissipation_grid/dissfloor);
end

% Add strain rate:
if strainweight~=0
    % Compute component:
    [StrainRateXX_grid,DUDY_grid]=gradient(U_obs_grid,dx_grid,dx_grid);
    [DVDX_grid,StrainRateYY_grid]=gradient(V_obs_grid,dx_grid,dx_grid);
    StrainRateXY_grid=.5*(DUDY_grid+DVDX_grid);
    StrainRateMag_grid=sqrt(StrainRateXY_grid.^2+.5*(StrainRateXX_grid.^2+StrainRateYY_grid.^2));
    % Apply mask constraints:
    StrainRateMag_grid(ReorderedMask_grid<.5|ReorderedMask_grid>2.5)=0; % 0 outside the ice sheet or shelf
    % Apply floor:
    StrainRateMag_grid=sqrt(StrainRateMag_grid.^2+strainfloor^2);
    % Add to control field:
    ControlField_grid(ReorderedMask_grid==1)=ControlField_grid(ReorderedMask_grid==1)+strainweight*log10(StrainRateMag_grid(ReorderedMask_grid==1)/strainfloor);
end

% Add velocity:
if velweight~=0
    % Apply floor:
    VelMag_obs_grid=sqrt(VelMag_obs_grid.^2+velfloor^2);
    % Apply mask constraints:
    VelMag_obs_grid(ReorderedMask_grid<.5|ReorderedMask_grid>2.5)=0; % 0 outside the ice sheet or shelf
    % Add to control field:
    ControlField_grid=ControlField_grid+velweight*log10(VelMag_obs_grid/velfloor);
end

% Add mask:
if maskweight~=0
    ControlField_grid(ReorderedMask_grid>0)=ControlField_grid(ReorderedMask_grid>0)+maskweight*(ReorderedMask_grid(ReorderedMask_grid>0)+frontmaskshift);
end

% Smooth out control field:
if controlfieldwavelength>0
    ControlField_grid=ThreeByThree(ControlField_grid,round((.5*controlfieldwavelength/dx_grid)^2));
end

% figure
% imagesc(X_grid/1000,Y_grid/1000,ControlField_grid)
% set(gca,'ydir','normal')
% colormap('jet')
% error('end here')

%% Load Boundary, modify if necessary:

% Communicate:
disp('Getting domain boundary, making initial mesh...')

% Load matlab version of the boundary file:
load([boundaryfile(1:end-4),'.mat'],'X_path','Y_path','BoundaryNum_path')

% Check whether to make new boundary:
if dist2expand>0
    
    % Identify first and last points of the ice front:
    ind1=find(BoundaryNum_path==2,1,'first');
    ind2=find(BoundaryNum_path==2,1,'last');
    
    % Compute angle of the front:
    frontangle=atan2(Y_path(ind2)-Y_path(ind1),X_path(ind2)-X_path(ind1));
    
    % Generate corners of new mask boundary:
    xnew=X_path([ind1,ind2])+dist2expand*cos(frontangle+pi*expandangle/180);
    ynew=Y_path([ind1,ind2])+dist2expand*sin(frontangle+pi*expandangle/180);
    
    % Make new boundary path:
    X_path_new=[X_path(1:ind1);xnew;X_path(ind2:end)];
    Y_path_new=[Y_path(1:ind1);ynew;Y_path(ind2:end)];
    BoundaryNum_path_new=[BoundaryNum_path(1:ind1);3;3;BoundaryNum_path(ind2:end)];
    
    % Save new argus file:
    CoordStruct=struct('x',X_path_new,'y',Y_path_new);
    newboundaryfile=[boundaryfile(1:end-4),'_expanded.exp'];
    expwrite(CoordStruct,newboundaryfile)
    
    % Assign this boundary file:
    thisboundaryfile=newboundaryfile;
    
else
    % Assign this boundary file:
    thisboundaryfile=boundaryfile;
end

% Generate an initial uniform mesh (resolution = hinit m)
initialmd=bamg(model,'domain',thisboundaryfile,'hmax',hinit);

%% Loop through samples:

% Communicate:
disp('Looping through samples')

for thissample=1:numsamples

    %% Prepare this sample:

    % Communicate:
    disp('...')
    disp(['Generating mesh number ',num2str(thissample),'/',num2str(numsamples)])
    disp('...')

    % Assign initial mesh to the main mesh:
    md=initialmd;

    % Assign parameter values:
    err=allerr(thissample);
    hmax=allhmax(thissample,:);
    hmin=allhmin(thissample);
    MeshingParameters.err=err;
    MeshingParameters.hmax=hmax;
    MeshingParameters.hmin=hmin;

    % Set output file:
    thisoutputfile=[generaloutputfile,'_mesh',num2str(thissample)];

    % Close figures:
    close all

    %% Iteratively Refine Mesh:

    % Communicate:
    disp('Iteratively refining mesh...')

    % Plot:
    if showintermediateiterations
        figure(1)
        patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','k','FaceColor','none')
        title('Iteration 0')
    end

    % Iterate on mesh refinement:
    for ii=1:numiterations_mesh

        % Interpolate control field onto mesh:
        ControlField_mesh=interp2(X_grid',Y_grid,ControlField_grid,md.mesh.x,md.mesh.y,'linear',0);

        % Interpolate mask to the mesh:
        Mask_mesh=interp2(X_grid',Y_grid,ReorderedMask_grid,md.mesh.x,md.mesh.y,'linear',0);

        % Assign max sizes to the mesh:
        MaxSize_mesh=zeros(size(md.mesh.x));
        MaxSize_mesh(Mask_mesh>=1.5&Mask_mesh<=2.5)=hmax(1);
        MaxSize_mesh(Mask_mesh>=0.5&Mask_mesh<1.5)=hmax(2);
        MaxSize_mesh(Mask_mesh<0.5|Mask_mesh>2.5)=hmax(3);

        % Adapt the mesh:
        md=bamg(md,'gradation',gradation,'field',ControlField_mesh,'err',err,'hmin',hmin,'hmaxVertices',MaxSize_mesh);

        % Plot:
        if showintermediateiterations
            figure(ii+1)
            patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','k','FaceColor','none')
            title(['Iteration ',num2str(ii)])
        end

        % Pause to allow bamg output to reach command line:
        pause(.01)

    end

    % Interpolate control field onto mesh one last time:
    ControlField_mesh=interp2(X_grid',Y_grid,ControlField_grid,md.mesh.x,md.mesh.y,'linear',0);

    % Plot mesh:
    %plotmodel(md,'data','mesh')
    %patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','k','FaceColor','none')

    % Interpolate control field onto refined mesh:
    %ControlField_mesh=InterpFromGridToMesh(X_grid',Y_grid,ControlField_grid,md.mesh.x,md.mesh.y,0);
    % Plot control field on refined mesh:
    %patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'FaceVertexCData',ControlField_mesh,'EdgeColor','none','FaceColor','interp')

    %% Interpolate Grids onto the Mesh:

    % Communicate:
    disp('Interpolating grids onto the mesh...')

    % Check if we"re doing the firn correction:
    if dofirncorrection
        ThisIceThick_grid=IceThick_grid-Firn_grid;
    else
        ThisIceThick_grid=IceThick_grid;
    end

    % Extrapolate ice thickness into the ocean:
    % identify bad values:
    BadValues_grid=ReorderedMask_grid<.5;
    % expand the bad values region by 2 in all directions:
    BadValues_grid(2:end-1,2:end-1)=BadValues_grid(1:end-2,1:end-2)|BadValues_grid(2:end-1,1:end-2)|BadValues_grid(3:end,1:end-2)|...
        BadValues_grid(1:end-2,2:end-1)|BadValues_grid(2:end-1,2:end-1)|BadValues_grid(3:end,2:end-1)|...
        BadValues_grid(1:end-2,3:end)|BadValues_grid(2:end-1,3:end)|BadValues_grid(3:end,3:end);
    BadValues_grid(2:end-1,2:end-1)=BadValues_grid(1:end-2,1:end-2)|BadValues_grid(2:end-1,1:end-2)|BadValues_grid(3:end,1:end-2)|...
        BadValues_grid(1:end-2,2:end-1)|BadValues_grid(2:end-1,2:end-1)|BadValues_grid(3:end,2:end-1)|...
        BadValues_grid(1:end-2,3:end)|BadValues_grid(2:end-1,3:end)|BadValues_grid(3:end,3:end);
    % Flag bad values and interpolate:
    ThisIceThick_grid(BadValues_grid)=NaN;
    ThisIceThick_grid=LaplacianInterpolation(ThisIceThick_grid,9);

    % Interpolate grids onto the mesh:
    BedElev_vertices=MultiWavelengthInterpolator(X_grid,Y_grid,BedElev_grid,md.mesh.x,md.mesh.y,md.mesh.elements);
    IceThick_vertices=MultiWavelengthInterpolator(X_grid,Y_grid,ThisIceThick_grid,md.mesh.x,md.mesh.y,md.mesh.elements);
    Mask_vertices=MultiWavelengthInterpolator(X_grid,Y_grid,ReorderedMask_grid,md.mesh.x,md.mesh.y,md.mesh.elements);

    % Compute level sets:
    IceLevelSet_vertices=1-2*(Mask_vertices<=2.5&Mask_vertices>=0.5);
    OceanLevelSet_vertices=1-2*(Mask_vertices<=1.5);

    % Set thickness to 0 outside of the ice sheet:
    IceThick_vertices(IceLevelSet_vertices>0)=0;

    % Force grounded areas to be at least 1m above hydrostatic equilibrium:
    BedElev_vertices(IceLevelSet_vertices<=0&OceanLevelSet_vertices>=0)=max(BedElev_vertices(IceLevelSet_vertices<=0&OceanLevelSet_vertices>=0),-(rho_i/rho_sw)*IceThick_vertices(IceLevelSet_vertices<=0&OceanLevelSet_vertices>=0)+1);

    % Force floating areas to be at least 1m below hydrostatic equilibrium:
    BedElev_vertices(IceLevelSet_vertices<=0&OceanLevelSet_vertices<0)=min(BedElev_vertices(IceLevelSet_vertices<=0&OceanLevelSet_vertices<0),-(rho_i/rho_sw)*IceThick_vertices(IceLevelSet_vertices<=0&OceanLevelSet_vertices<0)-1);

    % Compute surface elevation:
    SurfElev_vertices=max(BedElev_vertices+IceThick_vertices,(1-rho_i/rho_sw)*IceThick_vertices);

    % Compute ice base:
    Base_vertices=SurfElev_vertices-IceThick_vertices;

    % Ensure exact grounding:
    Base_vertices(OceanLevelSet_vertices>=0&IceLevelSet_vertices<=0)=BedElev_vertices(OceanLevelSet_vertices>=0&IceLevelSet_vertices<=0);

    % Re-Ensure exact ice thickness consistency:
    IceThick_vertices=SurfElev_vertices-Base_vertices;

    % Assign geometry and mask variables:
    md.geometry.surface=SurfElev_vertices;
    md.geometry.bed=BedElev_vertices;
    md.geometry.thickness=IceThick_vertices;
    md.geometry.base=Base_vertices;
    md.mask.ocean_levelset=OceanLevelSet_vertices;
    md.mask.ice_levelset=IceLevelSet_vertices;

    % Assign densities:
    md.materials.rho_ice=rho_i;
    md.materials.rho_water=rho_sw;

    % Save mesh:
    save([thisoutputfile,'.mat'],'md','MeshingParameters')

    %% Analyze Mesh:

    % Compute mesh properties:
    % Compute edge lengths:
    a=sqrt((md.mesh.x(md.mesh.elements(:,1))-md.mesh.x(md.mesh.elements(:,2))).^2+(md.mesh.y(md.mesh.elements(:,1))-md.mesh.y(md.mesh.elements(:,2))).^2);
    b=sqrt((md.mesh.x(md.mesh.elements(:,2))-md.mesh.x(md.mesh.elements(:,3))).^2+(md.mesh.y(md.mesh.elements(:,2))-md.mesh.y(md.mesh.elements(:,3))).^2);
    c=sqrt((md.mesh.x(md.mesh.elements(:,3))-md.mesh.x(md.mesh.elements(:,1))).^2+(md.mesh.y(md.mesh.elements(:,3))-md.mesh.y(md.mesh.elements(:,1))).^2);
    % Compute circumradius:
    CircumRadius=a.*b.*c./sqrt((a+b+c).*(b+c-a).*(c+a-b).*(a+b-c));
    % Compute minimum, maximum, and mean edge lengths:
    MinEdgeLength=min([a,b,c],[],2);
    MaxEdgeLength=max([a,b,c],[],2);
    MeanEdgeLength=mean([a,b,c],2);
    % Compute area:
    s=(a+b+c)/2;
    Area=sqrt(s.*(s-a).*(s-b).*(s-c));
    % Compute hydraulic diameter:
    HydraulicDiameter=2*Area./s;
    % Compute altitude lengths:
    Altitude_a=2*Area./a;
    Altitude_b=2*Area./b;
    Altitude_c=2*Area./c;
    % Compute anisotropy:
    %Anisotropy=MaxEdgeLength./MinEdgeLength;
    Anisotropy=(sqrt(3)/2)*max([a./Altitude_a,b./Altitude_b,c./Altitude_c],[],2);
    % Compute area fractions:
    IceAreaFraction=mean(IceLevelSet_vertices(md.mesh.elements)<0,2);
    OceanAreaFraction=mean(OceanLevelSet_vertices(md.mesh.elements)<0,2);

    % Make histograms of the mesh properties (histograms wrt number of elements):
    [Counts_hd,Edges_hd]=histcounts(log(HydraulicDiameter/1e3),100,'Normalization','probability');  % log scale, units of km
    [Counts_cr,Edges_cr]=histcounts(log(CircumRadius/1e3),100,'Normalization','probability');  % log scale, units of km
    [Counts_min,Edges_min]=histcounts(log(MinEdgeLength/1e3),100,'Normalization','probability');  % log scale, units of km
    [Counts_max,Edges_max]=histcounts(log(MaxEdgeLength/1e3),100,'Normalization','probability');  % log scale, units of km
    [Counts_mean,Edges_mean]=histcounts(log(MeanEdgeLength/1e3),100,'Normalization','probability');  % log scale, units of km
    [Counts_area,Edges_area]=histcounts(log(Area/1e6),100,'Normalization','probability');  % log scale, units of km^2
    [Counts_ani,Edges_ani]=histcounts(Anisotropy,100,'Normalization','probability');  % linear scale, unitless

    % Make histograms of the mesh properties (histograms wrt domain area):
    % Interpolate velocity onto mesh:
    VelMag_vertices=interp2(X_grid,Y_grid,VelMag_obs_grid,md.mesh.x,md.mesh.y,'linear',0);
    VelMag_elements=mean(VelMag_vertices(md.mesh.elements),2);
    % Compute weighting function for fast-flowing areas:
    FastFlowWeight=ones(size(md.mesh.elements,1),1);
    FastFlowWeight(VelMag_elements<.5*velanalysisthrehsold)=0;
    FastFlowWeight(VelMag_elements>=.5*velanalysisthrehsold&VelMag_elements<2*velanalysisthrehsold)=...
        .5*(1+sin(.5*pi*log2(VelMag_elements(VelMag_elements>=.5*velanalysisthrehsold&VelMag_elements<2*velanalysisthrehsold)/velanalysisthrehsold)));
    FastFlowWeight=FastFlowWeight.*IceAreaFraction.*(1-OceanAreaFraction);
    % compute area totals:
    totalarea=sum(Area);
    totalicearea=sum(Area.*IceAreaFraction);
    totalfloatingarea=sum(Area.*OceanAreaFraction.*IceAreaFraction);
    totalgroundedarea=totalicearea-totalfloatingarea;
    totalfastarea=sum(Area.*FastFlowWeight);
    % Pre-allocate:
    Counts_hd_totalarea=zeros(100,1);
    Counts_hd_groundedarea=zeros(100,1);
    Counts_hd_floatingarea=zeros(100,1);
    Counts_hd_fastarea=zeros(100,1);
    LogHydraulicDiameter=log(HydraulicDiameter/1e3);
    % Loop through bins:
    for ii=1:100
        % Generate counts:
        Counts_hd_totalarea(ii)=sum(Area(LogHydraulicDiameter>=Edges_hd(ii)&LogHydraulicDiameter<Edges_hd(ii+1)));
        Counts_hd_groundedarea(ii)=sum(Area(LogHydraulicDiameter>=Edges_hd(ii)&LogHydraulicDiameter<Edges_hd(ii+1)).*...
            IceAreaFraction(LogHydraulicDiameter>=Edges_hd(ii)&LogHydraulicDiameter<Edges_hd(ii+1)).*(1-OceanAreaFraction(LogHydraulicDiameter>=Edges_hd(ii)&LogHydraulicDiameter<Edges_hd(ii+1))));
        Counts_hd_floatingarea(ii)=sum(Area(LogHydraulicDiameter>=Edges_hd(ii)&LogHydraulicDiameter<Edges_hd(ii+1)).*...
            IceAreaFraction(LogHydraulicDiameter>=Edges_hd(ii)&LogHydraulicDiameter<Edges_hd(ii+1)).*OceanAreaFraction(LogHydraulicDiameter>=Edges_hd(ii)&LogHydraulicDiameter<Edges_hd(ii+1)));
        Counts_hd_fastarea(ii)=sum(Area(LogHydraulicDiameter>=Edges_hd(ii)&LogHydraulicDiameter<Edges_hd(ii+1)).*...
            FastFlowWeight(LogHydraulicDiameter>=Edges_hd(ii)&LogHydraulicDiameter<Edges_hd(ii+1)));
    end
    % Normalize to probablity density: (on log scale)
    Counts_hd_totalarea=Counts_hd_totalarea/(totalarea*(Edges_hd(2)-Edges_hd(1)));
    Counts_hd_groundedarea=Counts_hd_groundedarea/(totalgroundedarea*(Edges_hd(2)-Edges_hd(1)));
    Counts_hd_floatingarea=Counts_hd_floatingarea/(totalfloatingarea*(Edges_hd(2)-Edges_hd(1)));
    Counts_hd_fastarea=Counts_hd_fastarea/(totalfastarea*(Edges_hd(2)-Edges_hd(1)));


    %% Make Figures:

    % Communicate:
    disp('Making figures...')

    % Plot final mesh figure, if necessary:
    if showintermediateiterations==0
        figure(1)
        patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','k','FaceColor','none')
    end

    % Overlay grounding line and ice front from grid:
    hold on
    contour(X_grid/1000,Y_grid/1000,ReorderedMask_grid,.5:1:3.5,'r')

    % Compute axes limits:
    xlims=[min(md.mesh.x),max(md.mesh.x)]/1000;
    ylims=[min(md.mesh.y),max(md.mesh.y)]/1000;

    % Fix up final figure:
    xlim(xlims)
    ylim(ylims)
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    slashnum=strfind(thisoutputfile,'/');
    title(['Final Mesh for ',thisoutputfile(slashnum(end)+1:end)],'interpreter','none')
    % Save figure:
    thisfigname=[thisoutputfile,'_meshmap.png'];
    set(gcf,'PaperSize',pagesize1)
    set(gcf,'PaperPosition',[0,0,pagesize1])
    print('-dpng',thisfigname,['-r',num2str(resolution)])


    %% Make a 6-part map figure:
    figure
    % Bed topography:
    subplot(3,2,1)
    patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','none','FaceColor','interp','FaceVertexCData',md.geometry.bed)
    hold on
    contour(X_grid/1000,Y_grid/1000,ReorderedMask_grid,.5:1:3.5,'r')
    xlim(xlims)
    ylim(ylims)
    caxis([-2500,2500])
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','Elevation (m)')
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    title('a) Bed Topography')
    % Ice thickness:
    subplot(3,2,2)
    patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','none','FaceColor','interp','FaceVertexCData',md.geometry.surface)
    hold on
    contour(X_grid/1000,Y_grid/1000,ReorderedMask_grid,.5:1:3.5,'r')
    xlim(xlims)
    ylim(ylims)
    caxis([0,4000])
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','Elevation (m)')
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    title('b) Surface Topography')
    % Ice bottom topography:
    subplot(3,2,3)
    patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','none','FaceColor','interp','FaceVertexCData',md.geometry.base)
    hold on
    contour(X_grid/1000,Y_grid/1000,ReorderedMask_grid,.5:1:3.5,'r')
    xlim(xlims)
    ylim(ylims)
    caxis([-2500,2500])
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','Elevation (m)')
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    title('c) Ice Base Topography')
    % Ice thickness:
    subplot(3,2,4)
    patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','none','FaceColor','interp','FaceVertexCData',md.geometry.thickness)
    hold on
    contour(X_grid/1000,Y_grid/1000,ReorderedMask_grid,.5:1:3.5,'r')
    xlim(xlims)
    ylim(ylims)
    caxis([0,4000])
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','Thickness (m)')
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    title('d) Ice Thickness')
    % Ocean Level set:
    subplot(3,2,5)
    patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','none','FaceColor','interp','FaceVertexCData',md.mask.ocean_levelset)
    hold on
    contour(X_grid/1000,Y_grid/1000,ReorderedMask_grid,.5:1:3.5,'r')
    xlim(xlims)
    ylim(ylims)
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','Level Set Value')
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    title('e) Ocean Level Set')
    % Ice Level set:
    subplot(3,2,6)
    patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','none','FaceColor','interp','FaceVertexCData',md.mask.ice_levelset)
    hold on
    contour(X_grid/1000,Y_grid/1000,ReorderedMask_grid,.5:1:3.5,'r')
    xlim(xlims)
    ylim(ylims)
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','Level Set Value')
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    title('f) Ice Level Set')
    % Save figure:
    thisfigname=[thisoutputfile,'_6maps.png'];
    set(gcf,'PaperSize',pagesize3)
    set(gcf,'PaperPosition',[0,0,pagesize3])
    print('-dpng',thisfigname,['-r',num2str(resolution)])


    %% Make a 6-part histogram figure:

    figure

    % Top row (prob wrt element number)
    % Everything with length-units:
    subplot(2,3,1)
    histogram('BinEdges',exp(Edges_hd),'BinCounts',Counts_hd,'DisplayStyle','stairs','EdgeColor','k')
    hold on
    histogram('BinEdges',exp(Edges_cr),'BinCounts',Counts_cr,'DisplayStyle','stairs','EdgeColor',[.6,.3,0])
    %histogram('BinEdges',exp(Edges_min),'BinCounts',Counts_min,'DisplayStyle','stairs','EdgeColor','b')
    histogram('BinEdges',exp(Edges_mean),'BinCounts',Counts_mean,'DisplayStyle','stairs','EdgeColor','g')
    %histogram('BinEdges',exp(Edges_max),'BinCounts',Counts_max,'DisplayStyle','stairs','EdgeColor','r')
    set(gca,'XScale','log')
    theseylims=get(gca,'Ylim');
    plot((min(hmin)/1e3)*[1,1],theseylims,'--k')
    plot((max(hmax)/1e3)*[1,1],theseylims,'--k')
    text(min(hmin)/1e3,theseylims(1)+.7*(theseylims(2)-theseylims(1)),'hmin','HorizontalAlignment','left','VerticalAlignment','middle','Color','k')
    text(max(hmax)/1e3,theseylims(1)+.7*(theseylims(2)-theseylims(1)),'hmax','HorizontalAlignment','right','VerticalAlignment','middle','Color','k')
    xlim(histogramlengthlims)
    ylim(theseylims)
    xlabel('Length (km)')
    ylabel('Probability Density (wrt element count)')
    set(gca,'TickDir','out')
    legend('Hydraulic Diameter','Circumradius','Mean Edge Length','location','NorthEast')
    title('a) Mesh Element Length Measures')
    % Area:
    subplot(2,3,2)
    histogram('BinEdges',exp(Edges_area),'BinCounts',Counts_area,'DisplayStyle','stairs','EdgeColor','k')
    hold on
    set(gca,'XScale','log')
    theseylims=get(gca,'Ylim');
    plot((min(hmin)^2)*[1,1]/1e6,theseylims,'--k')
    plot((max(hmax)^2)*[1,1]/1e6,theseylims,'--k')
    text((min(hmin)^2)/1e6,theseylims(1)+.75*(theseylims(2)-theseylims(1)),'hmin^2','HorizontalAlignment','left','VerticalAlignment','middle','Color','k')
    text((max(hmax)^2/1e6),theseylims(1)+.75*(theseylims(2)-theseylims(1)),'hmax^2','HorizontalAlignment','right','VerticalAlignment','middle','Color','k')
    xlim(histogramarealims)
    ylim(theseylims)
    xlabel('Area (km^2)')
    ylabel('Probability Density (wrt element count)')
    set(gca,'TickDir','out')
    title('b) Mesh Element Area')
    % Anisotropy:
    subplot(2,3,3)
    histogram('BinEdges',Edges_ani,'BinCounts',Counts_ani,'DisplayStyle','stairs','EdgeColor','k')
    xlim(histogramanilims)
    xlabel('Max(side length / altitude)*sqrt(3)/2 (unitless)')
    ylabel('Probability Density (wrt element count)')
    set(gca,'TickDir','out')
    title('c) Mesh Element Anisotropy')

    % Bottom row:  (prob wrt domain area)
    % Grounded mesh size:
    subplot(2,3,4)
    histogram('BinEdges',exp(Edges_hd),'BinCounts',Counts_hd_groundedarea,'DisplayStyle','stairs','EdgeColor','k')
    hold on
    set(gca,'XScale','log')
    theseylims=get(gca,'Ylim');
    plot((hmin(1)/1e3)*[1,1],theseylims,'--k')
    plot((hmax(1)/1e3)*[1,1],theseylims,'--k')
    text(hmin(1)/1e3,theseylims(1)+.7*(theseylims(2)-theseylims(1)),'hmin','HorizontalAlignment','left','VerticalAlignment','middle','Color','k')
    text(hmax(1)/1e3,theseylims(1)+.7*(theseylims(2)-theseylims(1)),'hmax','HorizontalAlignment','right','VerticalAlignment','middle','Color','k')
    thiscdf=[0;cumsum(Counts_hd_groundedarea)/sum(Counts_hd_groundedarea)];
    [thisuniquecdf,thisind]=unique(thiscdf);
    thismedian=exp(interp1(thisuniquecdf,Edges_hd(thisind),.5,'linear'));
    plot(thismedian*[1,1],theseylims,'--r')
    text(thismedian,theseylims(1),['median= ',num2str(.01*round(100*thismedian)),' km'],'Color','k','VerticalAlignment','bottom','HorizontalAlignment','left','Rotation',90)
    xlim(histogramlengthlims)
    ylim(theseylims)
    xlabel('Length (km)')
    ylabel('Probability Density (wrt domain area)')
    set(gca,'TickDir','out')
    title('d) Hydraulic Diamter, Grounded Domain')
    % Floating mesh size:
    subplot(2,3,5)
    histogram('BinEdges',exp(Edges_hd),'BinCounts',Counts_hd_floatingarea,'DisplayStyle','stairs','EdgeColor','k')
    hold on
    set(gca,'XScale','log')
    theseylims=get(gca,'Ylim');
    plot((hmin/1e3)*[1,1],theseylims,'--k')
    plot((hmax(2)/1e3)*[1,1],theseylims,'--k')
    text(hmin/1e3,theseylims(1)+.7*(theseylims(2)-theseylims(1)),'hmin','HorizontalAlignment','left','VerticalAlignment','middle','Color','k')
    text(hmax(2)/1e3,theseylims(1)+.7*(theseylims(2)-theseylims(1)),'hmax','HorizontalAlignment','right','VerticalAlignment','middle','Color','k')
    thiscdf=[0;cumsum(Counts_hd_floatingarea)/sum(Counts_hd_floatingarea)];
    [thisuniquecdf,thisind]=unique(thiscdf);
    thismedian=exp(interp1(thisuniquecdf,Edges_hd(thisind),.5,'linear'));
    plot(thismedian*[1,1],theseylims,'--r')
    text(thismedian,theseylims(1),['median= ',num2str(.01*round(100*thismedian)),' km'],'Color','k','VerticalAlignment','bottom','HorizontalAlignment','left','Rotation',90)
    xlim(histogramlengthlims)
    ylim(theseylims)
    xlabel('Length (km)')
    ylabel('Probability Density (wrt domain area)')
    set(gca,'TickDir','out')
    title('e) Hydraulic Diameter, Floating Domain')
    % Grounded mesh size:
    subplot(2,3,6)
    histogram('BinEdges',exp(Edges_hd),'BinCounts',Counts_hd_fastarea,'DisplayStyle','stairs','EdgeColor','k')
    hold on
    set(gca,'XScale','log')
    theseylims=get(gca,'Ylim');
    plot((hmin(1)/1e3)*[1,1],theseylims,'--k')
    plot((hmax(1)/1e3)*[1,1],theseylims,'--k')
    text(hmin(1)/1e3,theseylims(1)+.7*(theseylims(2)-theseylims(1)),'hmin','HorizontalAlignment','left','VerticalAlignment','middle','Color','k')
    text(hmax(1)/1e3,theseylims(1)+.7*(theseylims(2)-theseylims(1)),'hmax','HorizontalAlignment','right','VerticalAlignment','middle','Color','k')
    thiscdf=[0;cumsum(Counts_hd_fastarea)/sum(Counts_hd_fastarea)];
    [thisuniquecdf,thisind]=unique(thiscdf);
    thismedian=exp(interp1(thisuniquecdf,Edges_hd(thisind),.5,'linear'));
    plot(thismedian*[1,1],theseylims,'--r')
    text(thismedian,theseylims(1),['median= ',num2str(.01*round(100*thismedian)),' km'],'Color','k','VerticalAlignment','bottom','HorizontalAlignment','left','Rotation',90)
    xlim(histogramlengthlims)
    ylim(theseylims)
    xlabel('Length (km)')
    ylabel('Probability Density (wrt domain area)')
    set(gca,'TickDir','out')
    title(['f) Hydraulic Diameter, Grounded Fast-Flow Domain (v>',num2str(round(velanalysisthrehsold)),' m/yr)'])

    % Save figure:
    thisfigname=[thisoutputfile,'_histograms.png'];
    set(gcf,'PaperSize',pagesize2)
    set(gcf,'PaperPosition',[0,0,pagesize2])
    print('-dpng',thisfigname,['-r',num2str(resolution)])


    %% Make a figure for the control field:
    figure
    % Control Field:
    imagesc(X_grid/1000,Y_grid/1000,ControlField_grid)
    set(gca,'ydir','normal')
    hold on
    if dist2expand>0
        plot(X_path_new/1000,Y_path_new/1000,'k')
    else
        plot(X_path/1000,Y_path/1000,'k')
    end
    %contour(X_grid/1000,Y_grid/1000,ReorderedMask_grid,.5:1:3.5,'r')
    caxis([min(ControlField_mesh),max(ControlField_mesh)])
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','Control Field Units')
    xlim(xlims)
    ylim(ylims)
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    title(['Control Field for ',thisoutputfile(slashnum(end)+1:end)],'interpreter','none')

    % Save figure:
    thisfigname=[thisoutputfile,'_controlfieldmap.png'];
    set(gcf,'PaperSize',pagesize1)
    set(gcf,'PaperPosition',[0,0,pagesize1])
    print('-dpng',thisfigname,['-r',num2str(resolution)])

    %% Make a figure for the mesh size:
    figure
    % Control Field:
    patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','none','FaceColor','flat','FaceVertexCData',log10(HydraulicDiameter/1000))
    hold on
    %contour(X_grid/1000,Y_grid/1000,ReorderedMask_grid,.5:1:3.5,'r')
    theseclims=log10(meshsizeclims);
    caxis(theseclims)
    % Make colorbar:
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','km')
    % Compute first tick for color bar:
    thisbase=10^floor(theseclims(1));
    thistick=thisbase*ceil((10^theseclims(1))/thisbase);
    % Create all ticks and labels for color bar:
    cticks=zeros(100,1);
    clabels=cell(100,1);
    for ii=1:100
        % Record this tick:
        cticks(ii)=log10(thistick);
        % Check if we hit the next power of ten:
        if abs(thistick-10*thisbase)/thisbase<1e-12 || (ii==1 && abs(thistick-thisbase)/thisbase<1e-12)
            % Generate label:
            clabels{ii}=['10^{',num2str(round(log10(thistick))),'}'];
            % Advance base:
            thisbase=thistick;
        end
        % Advance tick:
        thistick=thistick+thisbase;
        % Break from loop and truncate:
        if thistick>(1+1e-12)*(10^theseclims(2))
            cticks=cticks(1:ii);
            clabels=clabels(1:ii);
            break
        end
    end
    % Set the colorbar ticks and labels:
    set(hc,'YTick',cticks)
    set(hc,'YTickLabel',clabels)
    set(hc,'TickDir','out')
    % Fix up the rest:
    xlim(xlims)
    ylim(ylims)
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    title(['Mesh Size (hydraulic diameter) ',thisoutputfile(slashnum(end)+1:end)],'interpreter','none')

    % Save figure:
    thisfigname=[thisoutputfile,'_meshsizemap.png'];
    set(gcf,'PaperSize',pagesize1)
    set(gcf,'PaperPosition',[0,0,pagesize1])
    print('-dpng',thisfigname,['-r',num2str(resolution)])

end

% Final display:
disp('Done!')
toc

%% Cutting Room Floor:

% % Smooth out mask:
% if maskwavelength>0
%     SmoothReorderedMask_grid=ThreeByThree(ReorderedMask_grid,round((.5*maskwavelength/dx_grid)^2));
% end
%
% % Interpolate mask onto the mesh:
% Mask_mesh=InterpFromGridToMesh(X_grid',Y_grid,SmoothReorderedMask_grid,md.mesh.x,md.mesh.y,0);
%
% % Turn the mask into levelset variables:
% md.mask.ocean_levelset=-double(Mask_mesh<1.5);
% md.mask.ocean_levelset(md.mask.ocean_levelset>=0)=1;
% md.mask.ice_levelset=-double(Mask_mesh>.5&Mask_mesh<2.5);
% md.mask.ice_levelset(md.mask.ice_levelset>=0)=1;
%
% % Flag the ice front:  (assumes that boundary #2 is the front)
% % Loop through boundary points:
% for ii=1:length(X_path)-1
%     % Check if this is a front segment:
%     if BoundaryNum_path(ii)==2 && BoundaryNum_path(ii+1)==2
%         % Find boundary vertices in between these points:
%         TheseVertices=md.mesh.vertexonboundary==1&md.mesh.x>=min(X_path(ii:ii+1))&md.mesh.x<=max(X_path(ii:ii+1))&md.mesh.y>=min(Y_path(ii:ii+1))&md.mesh.y<=max(Y_path(ii:ii+1));
%         % Assume those vertices to the calving front:
%         md.mask.ice_levelset(TheseVertices)=0;
%     end
% end
%
% % Interpolate geometry onto the mesh:
% md.geometry.surface=InterpFromGridToMesh(X_grid',Y_grid,SurfElev_grid,md.mesh.x,md.mesh.y,0);
% md.geometry.thickness=InterpFromGridToMesh(X_grid',Y_grid,IceThick_grid,md.mesh.x,md.mesh.y,0);
% md.geometry.base=InterpFromGridToMesh(X_grid',Y_grid,SurfElev_grid-IceThick_grid,md.mesh.x,md.mesh.y,0);
% md.geometry.bed=InterpFromGridToMesh(X_grid',Y_grid,BedElev_grid,md.mesh.x,md.mesh.y,0);





% % Compute element centers:
% CentX=mean(md.mesh.x(md.mesh.elements),2);
% CentY=mean(md.mesh.y(md.mesh.elements),2);
% 
% % Pre-allocate interpolated products:
% Mask_elements=zeros(md.mesh.numberofelements,1);
% BedElev_elements=zeros(md.mesh.numberofelements,1);
% SurfElev_elements=zeros(md.mesh.numberofelements,1);
% IceThick_elements=zeros(md.mesh.numberofelements,1);
% 
% % Pre-allocate intermediate products:
% CloseGridCells=false(ysize,xsize);
% Dist=zeros(ysize,xsize);
% Weights=zeros(ysize,xsize);
% 
% % Make full grid:
% [X_grid_full,Y_grid_full]=meshgrid(X_grid,Y_grid);
% 
% % Loop through elements:
% for ii=1:md.mesh.numberofelements
%     % Compute this radius:
%     thisradius=HydraulicDiameter(ii)/2;
%     % Identify close elements:
%     CloseGridCells=abs(X_grid_full-CentX(ii))<4*thisradius&abs(Y_grid_full-CentY(ii))<4*thisradius;
%     % Compute distance:
%     Dist(CloseGridCells)=sqrt((X_grid_full(CloseGridCells)-CentX(ii)).^2+(Y_grid_full(CloseGridCells)-CentY(ii)).^2);
%     % Compute weighting:
%     Weights(CloseGridCells)=exp(-.5*((Dist(CloseGridCells)/thisradius).^2));
%     % Compute weighted means:
%     IceThick_elements(ii)=sum(IceThick_grid(CloseGridCells).*Weights(CloseGridCells))/sum(Weights(CloseGridCells));
%     BedElev_elements(ii)=sum(BedElev_grid(CloseGridCells).*Weights(CloseGridCells))/sum(Weights(CloseGridCells));
%     SurfElev_elements(ii)=sum(SurfElev_grid(CloseGridCells).*Weights(CloseGridCells))/sum(Weights(CloseGridCells));
%     Mask_elements(ii)=sum(ReorderedMask_grid(CloseGridCells).*Weights(CloseGridCells))/sum(Weights(CloseGridCells));
% end


% % Compute which wavelength to use for each element:
% % This is computed by taking double the hydraulic radius and then going to
% % the nearest wavelength on a log scale
% % pre-allocate:
% Wavelength2Use=zeros(md.mesh.numberofelements,1);
% % Compute cutoff wavelengths:
% cutoffwavelengths=[0,sqrt(interpwavelengths(1:end-1).*interpwavelengths(2:end)),max(HydraulicDiameter*2)+1];
% % Loop through wavelengths:
% for ii=1:numwavelengths
%     % Assign elements to this wavelength:
%     Wavelength2Use(2*HydraulicDiameter>=cutoffwavelengths(ii)&2*HydraulicDiameter<cutoffwavelengths(ii+1))=ii;
% end
% 
% % Compute element centers:
% CentX=mean(md.mesh.x(md.mesh.elements),2);
% CentY=mean(md.mesh.y(md.mesh.elements),2);
% 
% % Pre-allocate interpolated products:
% Mask_elements=zeros(md.mesh.numberofelements,1);
% BedElev_elements=zeros(md.mesh.numberofelements,1);
% SurfElev_elements=zeros(md.mesh.numberofelements,1);
% IceThick_elements=zeros(md.mesh.numberofelements,1);
% 
% % Loop through wavelengths and interpolate:
% for ii=1:numwavelengths
%     Mask_elements(Wavelength2Use==ii)=interp2(X_grid,Y_grid,squeeze(Mask_smoothed_grid(:,:,ii)),CentX(Wavelength2Use==ii),CentY(Wavelength2Use==ii));
%     BedElev_elements(Wavelength2Use==ii)=interp2(X_grid,Y_grid,squeeze(BedElev_smoothed_grid(:,:,ii)),CentX(Wavelength2Use==ii),CentY(Wavelength2Use==ii));
%     SurfElev_elements(Wavelength2Use==ii)=interp2(X_grid,Y_grid,squeeze(SurfElev_smoothed_grid(:,:,ii)),CentX(Wavelength2Use==ii),CentY(Wavelength2Use==ii));
%     IceThick_elements(Wavelength2Use==ii)=interp2(X_grid,Y_grid,squeeze(IceThick_smoothed_grid(:,:,ii)),CentX(Wavelength2Use==ii),CentY(Wavelength2Use==ii));
% end
% 
% % Average from elements onto vertices:
% % pre-allocate:
% Mask_vertices=zeros(md.mesh.numberofvertices,1);
% BedElev_vertices=zeros(md.mesh.numberofvertices,1);
% SurfElev_vertices=zeros(md.mesh.numberofvertices,1);
% IceThick_vertices=zeros(md.mesh.numberofvertices,1);
% Denom_vertices=zeros(md.mesh.numberofvertices,1);
% % loop through elements:
% for ii=1:md.mesh.numberofelements
%     % Add to the running sums:
%     Mask_vertices(md.mesh.elements(ii,:))=Mask_vertices(md.mesh.elements(ii,:))+Mask_elements(ii);
%     BedElev_vertices(md.mesh.elements(ii,:))=BedElev_vertices(md.mesh.elements(ii,:))+BedElev_elements(ii);
%     SurfElev_vertices(md.mesh.elements(ii,:))=SurfElev_vertices(md.mesh.elements(ii,:))+SurfElev_elements(ii);
%     IceThick_vertices(md.mesh.elements(ii,:))=IceThick_vertices(md.mesh.elements(ii,:))+IceThick_elements(ii);
%     Denom_vertices(md.mesh.elements(ii,:))=Denom_vertices(md.mesh.elements(ii,:))+1;
% end
% % finish computing means:
% Mask_vertices=Mask_vertices./Denom_vertices;
% BedElev_vertices=BedElev_vertices./Denom_vertices;
% SurfElev_vertices=SurfElev_vertices./Denom_vertices;
% IceThick_vertices=IceThick_vertices./Denom_vertices;
% 
% % Turn the mask into levelset variables:
% md.mask.ocean_levelset=-double(Mask_vertices<1.5);
% md.mask.ocean_levelset(md.mask.ocean_levelset>=0)=1;
% md.mask.ice_levelset=-double(Mask_vertices>.5&Mask_vertices<2.5);
% md.mask.ice_levelset(md.mask.ice_levelset>=0)=1;
% 
% % Flag the ice front:  (assumes that boundary #2 is the front)
% % Loop through boundary points:
% for ii=1:length(X_path)-1
%     % Check if this is a front segment:
%     if BoundaryNum_path(ii)==2 && BoundaryNum_path(ii+1)==2
%         % Find boundary vertices in between these points:
%         TheseVertices=md.mesh.vertexonboundary==1&md.mesh.x>=min(X_path(ii:ii+1))&md.mesh.x<=max(X_path(ii:ii+1))&md.mesh.y>=min(Y_path(ii:ii+1))&md.mesh.y<=max(Y_path(ii:ii+1));
%         % Assume those vertices to the calving front:
%         md.mask.ice_levelset(TheseVertices)=0;
%     end
% end
% 
% % Assign geometry variables:
% md.geometry.surface=SurfElev_vertices;
% md.geometry.bed=BedElev_vertices;
% md.geometry.thickness=IceThick_vertices;
% md.geometry.base=SurfElev_vertices-IceThick_vertices;
% 
% % Ensure exact grounding:
% md.geometry.base(md.mask.ocean_levelset==1)=md.geometry.bed(md.mask.ocean_levelset==1);
% 
% % Ensure exact ice thickness consistency:
% md.geometry.thickness=md.geometry.surface-md.geometry.base;


% % Fill in holes in observed velocity with neighboring values:
% if numiterations_velfill>0
%     % Remember old values:
%     OldVelMag_obs_grid=VelMag_obs_grid;
%     % Identify area within the boundary:
%     InDomain_grid=reshape(inpolygon(reshape(repmat(X_grid,[ysize,1]),[],1),reshape(repmat(Y_grid,[1,xsize]),[],1),X_path,Y_path),[ysize,xsize]);
%     for ii=1:numiterations_velfill
%         % Identify good values of velocity:
%         GoodValues_grid=isnan(VelMag_obs_grid)==0&VelMag_obs_grid>mingoodvalue4velfill;
%         % Loop through all grid cells:
%         for jj=1:ysize*xsize
%             % Get subscript indices:
%             [d1,d2]=ind2sub([ysize,xsize],jj);
%             % Check whether to skip this grid cell:
%             if GoodValues_grid(d1,d2) || InDomain_grid(d1,d2)==0
%                 continue
%             end
%             % Check if any neighbors are good values:
%             NeighborGoodValues=GoodValues_grid(max(1,d1-1):min(ysize,d1+1),max(1,d2-1):min(xsize,d2+1));
%             if sum(NeighborGoodValues(:))==0
%                 continue
%             end
%             % Take average of neighbor values:
%             NeighborVelMag=VelMag_obs_grid(max(1,d1-1):min(ysize,d1+1),max(1,d2-1):min(xsize,d2+1));
%             VelMag_obs_grid(d1,d2)=mean(NeighborVelMag(NeighborGoodValues));
%         end
%     end
% end
% 
% % Fill in remaining holes in velocity:
% VelMag_obs_grid(isnan(VelMag_obs_grid)|VelMag_obs_grid<velfiller)=velfiller;
% 
% % Apply soft floor to velocity:
% if velsoftfloor>0
%     VelMag_obs_grid=sqrt(VelMag_obs_grid.^2+velsoftfloor^2);
% end
% 
% % Smooth out velocity:
% if velwavelength>0
%     VelMag_obs_grid=ThreeByThree(VelMag_obs_grid,round((.5*velwavelength/dx_grid)^2));
% end


% % Extrapolate ice thickness into the ocean:
% % Idenify good values (ice-covered cells):
% GoodValues_grid=ReorderedMask_grid==1|ReorderedMask_grid==2;
% % Compute number of fill iterations:
% numfilliterations=1.5*ceil(hmax/dx_grid);
% % Initialized HasFilled variable:
% HasFilled_grid=false(ysize,xsize);
% for ii=1:numfilliterations
%     % Identify good values at this iteration:
%     GoodValues_grid=GoodValues_grid|HasFilled_grid;
%     % Loop through all grid cells:
%     for jj=1:ysize*xsize
%         % Get subscript indices:
%         [d1,d2]=ind2sub([ysize,xsize],jj);
%         % Check whether to skip this grid cell:
%         if GoodValues_grid(d1,d2) || ReorderedMask_grid(d1,d2)~=0
%             continue
%         end
%         % Check if any neighbors are good values:
%         NeighborGoodValues=GoodValues_grid(max(1,d1-1):min(ysize,d1+1),max(1,d2-1):min(xsize,d2+1));
%         if sum(NeighborGoodValues(:))==0
%             continue
%         end
%         % Take average of neighbor values:
%         NeighborIcethick=IceThick_grid(max(1,d1-1):min(ysize,d1+1),max(1,d2-1):min(xsize,d2+1));
%         IceThick_grid(d1,d2)=mean(NeighborIcethick(NeighborGoodValues));
%         % Record that this cell was filled:
%         HasFilled_grid(d1,d2)=1;
%     end
% end

% % Compute new ice level set:
% % Start with all ice:
% IceLevelSet_vertices=-1*ones(md.mesh.numberofvertices,1);
% % Flag the ice front:  (assumes that boundary #2 is the front)
% % Loop through boundary points:
% for ii=1:length(X_path)-1
%     % Check if this is a front segment:
%     if BoundaryNum_path(ii)==2 && BoundaryNum_path(ii+1)==2
%         % Find boundary vertices in between these points:
%         TheseVertices=md.mesh.vertexonboundary==1&md.mesh.x>=min(X_path(ii:ii+1))&md.mesh.x<=max(X_path(ii:ii+1))&md.mesh.y>=min(Y_path(ii:ii+1))&md.mesh.y<=max(Y_path(ii:ii+1));
%         % Assume those vertices to the calving front:
%         IceLevelSet_vertices(TheseVertices)=0;
%     end
% end
% % Remove rock:
% IceLevelSet_vertices(Mask_vertices>2.5)=1;
% 
% % Compute new ocean level set:
% OceanLevelSet_vertices=1-2*((Base_vertices-BedElev_vertices>groundingtolerance)&IceLevelSet_vertices<=0);