% ISSMInversion

% Mike Wolovick, 8/15/2021

% This script runs an inversion using ISSM for a basal friction
% coefficient and/or column-average rheology.

% v2:  The script now uses the MultiWavelengthInterpolator function and
% cuts out a lot of the slow repetitive code that did multiwavelength
% smoothing before.

% v2a:  This script flattens all unneccessarily 3D variables in its
% internal processing and when saving results to disc.  The model structure
% is only expanded when calling "solve".

% v3:  Various updates and streamlining for the inversion paper.  In SSA,
% the script can interpolate or solve for a 3D temperature field and then
% collapse the rheology back to a vertically averaged value.  Some of the
% complexity of old options is removed.  In addition, the overall
% functionality of the script relative to setup and loading results from
% cluster is handled differently.  The script now has two modes: 1) setup
% and send model to cluster, without waiting for results; and 2) load
% results from cluster and analyze them.  The script also now has the
% ability to use effective pressure in the sliding law.

% ISSM cost function codes:
% 101: SurfaceAbsVelMisfit
% 102: SurfaceRelVelMisfit
% 103: SurfaceLogVelMisfit
% 104: SurfaceLogVxVyMisfit
% 105: SurfaceAverageVelMisfit
% 201: ThicknessAbsMisfit
% 501: DragCoefficientAbsGradient
% 502: RheologyBbarAbsGradient
% 503: ThicknessAbsGradient

% My added codes:
% 102.1: homebrew SurfaceRelVelMisfit  (uses observed velocity magnitude to
%        normalize both components, rather than normalizing both separately
%        as the built-in does).

clearvars
close all
tic

% outputfiles={'/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_cuas_drag_v1_m3_mesh4.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_m3_mesh4.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_m5_mesh4.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_m3_mesh5.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_cuas_drag_v1a_mesh1.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_opc_v1_mesh1.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_op_v1_mesh1.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_opc_v1_mesh2.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_op_v1_mesh2.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_opc_v1_mesh4.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_op_v1_mesh4.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_opc_v1_mesh6.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_op_v1_mesh6.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_opc_v1_mesh8.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_op_v1_mesh8.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_HO_drag_v2_mesh4.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_HO_drag_v2_mesh5.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_HO_drag_v2_mesh6.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_HO_drag_v2_mesh8.mat'};
% outputfiles={'/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh1a.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh2.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh3.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh4.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh5.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh6.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh7.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh8.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh9.mat';...
%     '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh10.mat'};
 
 %for thislcurve=1:length(outputfiles)

%% Parameters:

% Do the startup:
%run('/home/ollie/mrueckam/modelling/ISSM/ISSM_startups/ISSM_JPL_v4p18_matlab_startup.m')

% File paths and names:
%issmscriptpath='~/Dropbox/ISSM/';
issmscriptpath='/home/ollie/mwolovic/ISSM_installations/';
%inputfile='~/Dropbox/FRISModeling/Meshes/FRISDomainMesh_SizeSeries_v1_mesh5.mat';
inputfile='/work/ollie/mwolovic/FRISModeling/Meshes/FRISDomainMesh_SizeSeries_v1_mesh6.mat';
%inputfile='/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_SizeSeries1Mesh5_SSA_v1_Lambda_1.mat';
%gridfile='~/Dropbox/CombinedGrids/AntarcticaCombinedGrids_v4.mat';
gridfile='/work/ollie/mwolovic/Inputs/AntarcticaCombinedGrids_v4.mat';
%velfile='~/Dropbox/CombinedGrids/AntarcticaCombinedGrids_v4.mat';
%velsuffix='_obs';
velfile='/work/ollie/mwolovic/Inputs/AntarcticaMergedVel_MEASURESv2withphase_plus_FRmosaic_velonly_900m.mat';
velsuffix='_merged';   % suffix appended to "U" and "V" in velocity file variable names, can be empty string
tempfile=[];
%tempfile='/work/ollie/tkleiner/PISM1Pal_RECISL_heatflux/post-processing/tzyx/04_run_regrid_16km_bm2_recisl_icegrav_qgeoM17_-010000_-_0000000.nc';
nfile='/work/ollie/mwolovic/Inputs/Ant1DMelt_spinup1360days_1000m_Sy0.4_Tvostok0.1_FRIS_output_cropped_last.nc';
%outputfile='~/Dropbox/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v1_mesh4_Lambda_0.562.mat';  % include trailing ".mat"
outputfile='/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh1a.mat';  % include trailing ".mat"
%outputfile=outputfiles{thislcurve}


% Overall control:
overallmode='loadandanalyze';         % 'setupandsend' or 'loadandanalyze'

% Ice rheology controls:
tempsource='const';                   % '1D', 'file', 'md', or 'const'
defaulticetemp=273.15-25;          % K
n=3;                               % unitless
rheologystyle='cuffeytemperate';   % 'paterson' or 'cuffeytemperate'
includewaterfraction=0;            % logical
minpismicethick=100;               % m
numtempsmoothings=4;               % integer
numshearheataveragings=5;          % integer
surftempshift=-10;                   % K
surftempshiftz0=200;               % m

% Sliding controls:
m=1;                               % unitless (sliding exponent, ub~sigma^m, or sigma~ub^1/m)
r=1;                                 % unitless (N exponent, sigma~N^r, or ub~N^-r*m)
nsource='none';                 % source of effective pressure: 'none', 'oceanp', 'oceanpwithcutoff', 'bedelev', or 'file'
minn=1e2;                       % Pa
numnaveragings=0;         % integer
numnunatakiterations=3;  % integer (for cutting out cells adjacent to nunataks in cuas output)

% Flow equation:
flowequation='HO';                % string (SIA, SSA, H0, L1L2, MLHO, or FS)

% What are we solving for?
includedrag=1;                     % logical
includerheology=0;                 % logical

% Cost function controls:
costfunctions=[101,501];           % [1xN] vector of ISSM cost function codes
intrinsicweights=[1,1];            % [1xN] vector of intrinsic weights (ie, weight after accounting for characteristic scales)
issmmeanvel=1000;                  % m/yr (arbitrary scaling velocity used internally in ISSM cost functions)
velabsscale='rms';                 % m/yr, 'mean', 'median', 'stdev', or 'rms'
vellogscale=2;                     % unitless>1, 'logstdev', '1+stdev/mean', '1+rms/mean'  (this quantity represents ratio before log is taken)
velstabilizer=.25;                 % m/yr, used in relative and log misfits 

% First guess controls:
dragcoeffguessstyle='compute';     % 'file', 'compute', 'previousinversion', 'previousspinup', or 'const'
dragguessaveragings=0;                 % integer

% Active area controls:
statsactiveareaonly=0;             % logical (only use the active area to compute statistics for scales, or use the whole (ice-covered) domain?)
rheologyfloatingonly=0;            % logical (only solve for bbar in the shelf?)
nummaskaveragings=0;           % integer (number of times to average the ice mask on the mesh?)

% L-curve analysis controls:
dolcurve=1;                        % logical
lambdamin=1e-3;                    % unitless
lambdamax=1e3;                     % unitless
numlambda=25;                      % integer

% Minimum thickness controls:
minicethick=10;                    % m
cutoutminthick=1;                  % logical (set areas at minthick to non-ice in IC)

% Vertical layers:
numlayers=10;                      % integer
densifyfactor=2;                   % unitless
inlandbcpower=5;                   % unitless exponent for shape function on inland bc (0 uses constant velocity)

% Should the rheology be modified to ensure that deformation velocity is
% less than observed velocity?  (not applied for SSA)
reducedefvel=0;                     % logical
numudefstressaveragings=3; % integer
numratioaveragings=1;            % integer
minratio=1e-2;                          % unitless<1t
velstabilizer_def=1e-2;            % m/yr

% Add noise to the data?
randnoiseamp=0;                    % m/yr
stripenoiseamp=0;                  % m/yr
stripenoiselength=1e5;             % m
stripenoisewidth=1e4;              % m

% Solver numerical controls:
numsteps=200;                      % integer (for inversion iteration)
dxmin=1e-3;                        % unitless (inversion convergence criteria, max normalized difference between two iterations)
gttol=1e-5;                        % unitless (inversion convergence criteria, reduction in gradient norm)
useincompleteadjoint=0;            % logical
solver_residue_threshold=NaN;      % unitless (NaN deactivates)
solvertype='asmstokes';                 % 'ilua', 'bjacobi', 'asmstokes', or 'mumps' (default)
restol=1e-4;                       % unitless (stress balance convergence criteria, mechanical residual tolerance)
reltol=1e-3;                        % unitless (stress balance convergence criteria, norm(du)/norm(u) tolerance)
abstol=1;                          % m/yr (stress balance convergence criteria, max(du) tolerance)
maxiter=200;                        % integer (stress balance iteration)
numtaskspernode=24;                % integer
numelementspertask=4e3;            % integer
rtime=3;                           % hours
lowmem=1;                          % logical

% Reasonable values of velocity, stress, and strain rate:
% (these parameters are used to compute max/min constraints or reasonable
% guesses for the inversion targets)
maxreasonablevel=2e4;               % m/yr
middlereasonablevel=1e2;            % m/yr
minreasonablevel=.1;                % m/yr
maxreasonablestress=3e5;            % Pa
middlereasonablestress=1e5;         % Pa
minreasonablestress_drag=1e2;            % Pa  1e2 for drag or 3e3 for bbar
minreasonablestress_bbar=3e3;            % Pa  1e2 for drag or 3e3 for bbar
maxreasonablestrainrate=1e3/5000;   % 1/yr
minreasonablestrainrate=.1/4000;    % 1/yr
minreasonablen=3e2;                        % Pa
middlereasonablen=1e5;                    % Pa
maxreasonablen=4e7;                       % Pa

% Seconds per year:
secondsperyear=60*60*24*365.25;     % s/yr

% LCurve fitting parameters:
numlambdasubsamples=1000;                % integer
numlambdawavelengths=50;              % integer
maxlambdawavelength=log(100);         % dist in log(cost) vs log(cost) space
outlierthreshold=2;                 % unitless>1
excludenonmonojobs=0;                % logical
excludenonmonojreg=0;                % logical
num2excludefromstart=0;                   % integer
lambdainterpmethod='linear';           % valid interp1 method
curvenormstyle='max';                % 'max' or 'quad'
curvevarstyle='var';                % 'var' or 'rms'

% Output figures controls:
makelcurvefig=1;                    % logical (if 0 it also skips l-curve analysis, loads data on best tradeoff from file)
makelcurvemapcomparisonfigs=0;         % logical
makelcurvefftcomparisonfig=0;         % logical
makelcurvedivergencecomparisonfig=0;         % logical
numbins=100;                        % integer
rawmaglims=[-50,50];                % [1x2 m/yr
normmaglims=[-1,1];                 % [1x2] uniltess
%rawmaglims=[-100,100];                % [1x2 m/yr
%normmaglims=[-.5,.5];                 % [1x2] uniltess
draglims=[0,150];                   % [1x2] kPa
dragcoeffratiolims=[-2,2];          % [1x2] log10(ratio)
viscositylims=[13,17];              % [1x2] log10(Pas)
bbarratiolims=[-2,1];               % [1x2] log10(ratio)
minvelobs4hist=1;                   % m/yr
lcurvepagesize=[20,10];           % [1x2] inches
vertaxesbuffer_lcurve=.06;          % unitless
horzaxesbuffer_lcurve=.04;          % unitless
titlebuffer_lcurve=.05;             % unitless
horznotextbuffer_lcurve=.02;        % unitless
lcurvefrac=.5;                      % unitless<1
widthperhist=5;                     % inches
heightperhist=3;                    % inches
widthpermap=6;                      % inches
heightpermap=4;                     % inches
resolution=300;                     % dpi
vertaxesbuffer_hist=.03;            % unitless
vertaxeslabelbuffer_hist=.04;       % unitless
titlebuffer_hist=.04;               % unitless
horzaxesbuffer_hist=.04;            % unitless
horzaxeslabelbuffer_hist=.03;       % unitless
horznotextbuffer_hist=.02;          % unitless
vertaxesbuffer_map=.01;             % unitless
vertaxeslabelbuffer_map=.01;        % unitless
titlebuffer_map=.01;                % unitless
horzaxesbuffer_map=.01;             % unitless
horzaxeslabelbuffer_map=.02;        % unitless
horznotextbuffer_map=.005;          % unitless

% Make a structure with all of my inversion parameters:
MyInversionParameters=struct('meshfile',inputfile,'velfile',velfile,'tempfile',tempfile,'gridfile',gridfile,'nfile',nfile,...
    'tempsource',tempsource,'defaulticetemp',defaulticetemp,'n',n,'rheologystyle',rheologystyle,'includewaterfraction',includewaterfraction,...
    'minpismicethick',minpismicethick,'numtempsmoothings',numtempsmoothings,'numshearheataveragings',numshearheataveragings,'surftempshift',surftempshift,'surftempshiftz0',surftempshiftz0,...
    'm',m,'r',r,'nsource',nsource,'minn',minn,'numnaveragings',numnaveragings,'numnunatakiterations',numnunatakiterations,'flowequation',flowequation,...
    'includedrag',includedrag,'includerheology',includerheology,'costfunctions',costfunctions,'intrinsicweights',intrinsicweights,'issmmeanvel',issmmeanvel,...
    'velabsscale',velabsscale,'vellogscale',vellogscale,'velstabilizer',velstabilizer,'dragcoeffguessstyle',dragcoeffguessstyle,'dragguessaveragings',dragguessaveragings,...
    'statsactiveareaonly',statsactiveareaonly,'rheologyfloatingonly',rheologyfloatingonly,'nummaskaveragings',nummaskaveragings,...
    'dolcurve',dolcurve,'lambdamin',lambdamin,'lambdamax',lambdamax,'numlambda',numlambda,...
    'minicethick',minicethick,'cutoutminthick',cutoutminthick,...
    'numlayers',numlayers,'densifyfactor',densifyfactor,'inlandbcpower',inlandbcpower,...
    'reducedefvel',reducedefvel,'numudefstressaveragings',numudefstressaveragings,'numratioaveragings',numratioaveragings,'minratio',minratio,...
    'randnoiseamp',randnoiseamp,'stripenoiseamp',stripenoiseamp,'stripenoiselength',stripenoiselength,'stripenoisewidth',stripenoisewidth,...
    'numsteps',numsteps,'dxmin',dxmin,'gttol',gttol,'useincompleteadjoint',useincompleteadjoint,'solver_residue_threshold',solver_residue_threshold,'solvertype',solvertype,...
    'restol',restol,'reltol',reltol,'abstol',abstol,'maxiter',maxiter,'numtaskspernode',numtaskspernode,'numelementspertask',numelementspertask,'rtime',rtime,'lowmem',lowmem,...
    'maxreasonablevel',maxreasonablevel,'middlereasonablevel',middlereasonablevel,'minreasonablevel',minreasonablevel,...
    'maxreasonablestress',maxreasonablestress,'middlereasonablestress',middlereasonablestress,'minreasonablestress_drag',minreasonablestress_drag,...
    'minreasonablestress_bbar',minreasonablestress_bbar,'maxreasonablestrainrate',maxreasonablestrainrate,'minreasonablestrainrate',minreasonablestrainrate,...
    'minreasonablen',minreasonablen,'middlereasonablen',middlereasonablen,'maxreasonablen',maxreasonablen,'secondsperyear',secondsperyear);


%% Load Model:

% Communicate:
disp('Loading and analyzing model mesh...')

% Check overall mode:
if ~strcmp(overallmode,'setupandsend') && ~strcmp(overallmode,'loadandanalyze')
    error('Parameter "overallmode" must be a string equal to "setupandsend" or "loadandanalyze"')
end
% Check the drag coefficient guess style:
if strcmp(dragcoeffguessstyle,'file')
    if exist('dragcoefffile','var')==0
        error('If parameter "dragcoeffguessstyle" is set to "file", then variable "dragcoefffile" must exist.')
    end
elseif strcmp(dragcoeffguessstyle,'compute')==0 && strcmp(dragcoeffguessstyle,'previousinversion')==0 && strcmp(dragcoeffguessstyle,'previousspinup')==0 && strcmp(dragcoeffguessstyle,'const')==0 && strcmp(dragcoeffguessstyle,'file')==0
    error('Parameter "dragcoeffguessstyle" must be a string equal to "file", "compute", "previousinversion", "previousspinup", or "const".')
end
% Check temperature source:
if strcmp(tempsource,'file')
    if isempty(tempfile)
        error('If parameter "tempsource" is set to "file", then parameter "tempfile" cannot be empty.')
    end
elseif ~strcmp(tempsource,'1D') && ~strcmp(tempsource,'md')  && ~strcmp(tempsource,'const')
    error('Parameter "tempsource" must be a string equal to "file", "1D", "md", or "const".')
end
% Check rheology style:
if strcmp(rheologystyle,'cuffeytemperate')==0 && strcmp(rheologystyle,'paterson')==0
    error('Parameter "rheologystyle" must be a string equal to "cuffeytemperate" or "paterson".')
end
% Check regularization terms:
if isempty(find(costfunctions==501,1))==0 && includedrag==0
    error('If cost function 501 is active, drag must be included in the inversion.')
end
if isempty(find(costfunctions==502,1))==0 && includerheology==0
    error('If cost function 502 is active, rheology must be included in the inversion.')
end
% Check that weights and cost functions are the same length:
if length(costfunctions)~=length(intrinsicweights)
    error('Parameter "costfunctions" must be the same length as parameter "intrinsicweights".')
end
% Check sum of weights for L-curve analysis:
if dolcurve
    % Identify which elements of the cost function are regularization:
    AreReg=costfunctions>500;
    % Check sum of intrinsic weights:
    if abs(sum(intrinsicweights(AreReg==0))-1)>1e-12
        error('If doing an L-curve analysis, the sum of the intrinsic weights of data terms in the cost function must be 1.')
    end
    if abs(sum(intrinsicweights(AreReg))-1)>1e-12
        error('If doing an L-curve analysis, the sum of the intrinsic weights of regularization terms in the cost function must be 1.')
    end
end

% Add the script path:
if exist('issmscriptpath','var')
    addpath(genpath(issmscriptpath))
end

% Check whether to load the input or the output file:
if strcmp(overallmode,'loadandanalyze')
    % Load the output file:
    load(outputfile,'md','*Parameters','BasicModel','Lambda','DataCost','RegularizationCost','*best*','hassolved')
    if exist('BasicModel','var')
        md=BasicModel;
    end
    % Overwrite a few variables from parameter structure:
    if isfield(MyInversionParameters,'lambdasmoothing')
        MyInversionParameters=rmfield(MyInversionParameters,'lambdasmoothing');
    end
    if isfield(MyInversionParameters,'outlierthreshold')
        MyInversionParameters=rmfield(MyInversionParameters,'outlierthreshold');
    end
    % Unpack inversion parameter structure (overwrites existing parameters):
    unpack(MyInversionParameters)
    % Selectively unpack coefficient conversion parameters:
    if exist('CoeffConversionParameters','var')
        thesefields=fieldnames(CoeffConversionParameters);
        for ii=1:length(thesefields)
            eval(['thisvariable=CoeffConversionParameters.',thesefields{ii},';'])
            if strcmp(thesefields{ii}(1:3),'new') && ~isempty(thisvariable)
                eval([thesefields{ii}(4:end),'=thisvariable;'])
            end
        end
    end
else
    % Load the input file:
    load(inputfile,'md','hassolved')
    % Check if we're loading a spinup run:
    if isfield(md.results,'TransientSolution')
        % Flag that we're starting from a previous spinup and communicate:
        usepreviousspinup=1;
        MyInversionParameters.usepreviousspinup=usepreviousspinup;
        disp('...detected transient output, using final model geometry')
        % Check if the spinup model reached a solution:
        if exist('hassolved','var')
            if ~hassolved
                error('Previous spinup did not reach a solution')
            end
        end
        % Assign final geometry to new geometry:
        md.geometry.thickness=md.results.TransientSolution(end).Thickness;
        md.geometry.surface=md.results.TransientSolution(end).Surface;
        md.geometry.base=md.results.TransientSolution(end).Base;
        % Assign final mask to new mask:
        md.mask.ocean_levelset=max(-1,min(1,md.results.TransientSolution(end).MaskOceanLevelset));
        if md.calving.ismovingfront
            md.mask.ice_levelset=max(-1,min(1,md.results.TransientSolution(end).MaskIceLevelset));
        end
        % Force exact grounding:
        md.geometry.base(md.mask.ocean_levelset==1)=md.geometry.bed(md.mask.ocean_levelset==1);
    end
    % Force consistent ice thickness:
    md.geometry.thickness=md.geometry.surface-md.geometry.base;
    % Impose minimum ice thickness, adjust base if necessary:
    md.geometry.thickness=max(md.geometry.thickness,minicethick);
    md.geometry.base=md.geometry.surface-md.geometry.thickness;
    md.geometry.bed=min(md.geometry.bed,md.geometry.base);
    % Cut out areas at minthick:
    if cutoutminthick
        md.mask.ice_levelset(md.geometry.thickness<minicethick+1e-3)=1;
    end
    % Adjust vertical coordinate for consistency:
    try
        oldnumlayers=md.mesh.numberoflayers;
    catch
        oldnumlayers=1;
    end
    if oldnumlayers>1
        % Compute Zhat:
        Zhat=(md.mesh.z(1:md.mesh.numberofvertices2d:end)-md.mesh.z(1))/(md.mesh.z((md.mesh.numberoflayers-1)*md.mesh.numberofvertices2d+1)-md.mesh.z(1));
        % Adjust mesh vertical coordinate:
        md.mesh.z=md.geometry.base+(md.geometry.surface-md.geometry.base).*reshape(repmat(Zhat',[md.mesh.numberofvertices2d,1]),[],1);
        % Ensure top layer is exactly at the surface:
        md.mesh.z(end-md.mesh.numberofvertices2d+1:end)=md.geometry.surface(end-md.mesh.numberofvertices2d+1:end);
    end
end

% Compute guesses and limits:
if includedrag
    mindragcoeff=sqrt(minreasonablestress_drag/((maxreasonablevel/secondsperyear)^(1/m)));
    maxdragcoeff=sqrt(maxreasonablestress/((minreasonablevel/secondsperyear)^(1/m)));
    if ~strcmp(nsource,'none')
        mindragcoeff=mindragcoeff/(maxreasonablen^(r/2));
        maxdragcoeff=maxdragcoeff/(minreasonablen^(r/2));
    end
    MyInversionParameters.mindragcoeff=mindragcoeff;
    MyInversionParameters.maxdragcoeff=maxdragcoeff;
end
if includerheology
    minbbar=minreasonablestress_bbar/((maxreasonablestrainrate/secondsperyear)^(1/n));
    maxbbar=maxreasonablestress/((minreasonablestrainrate/secondsperyear)^(1/n));
    MyInversionParameters.minbbar=minbbar;
    MyInversionParameters.maxbbar=maxbbar;
end
if strcmp(dragcoeffguessstyle,'const')
    dragcoeffguess=sqrt(middlereasonablestress/((middlereasonablevel/secondsperyear)^(1/m)));
    if ~strcmp(nsource,'none')
        dragcoeffguess=dragcoeffguess/(middlereasonablen^(r/2));
    end
    MyInversionParameters.dragcoeffguess=dragcoeffguess;
end

% Compute mesh properties:
% Extract 2D mesh:
try
    x2d=md.mesh.x2d;
    y2d=md.mesh.y2d;
    elements2d=md.mesh.elements2d;
    oldnumlayers=md.mesh.numberoflayers;
catch
    x2d=md.mesh.x;
    y2d=md.mesh.y;
    elements2d=md.mesh.elements;
    oldnumlayers=1;
end
numelements2d=size(elements2d,1);
numvertices2d=length(x2d);
% Compute edge lengths:
a=sqrt((x2d(elements2d(:,1))-x2d(elements2d(:,2))).^2+(y2d(elements2d(:,1))-y2d(elements2d(:,2))).^2);
b=sqrt((x2d(elements2d(:,2))-x2d(elements2d(:,3))).^2+(y2d(elements2d(:,2))-y2d(elements2d(:,3))).^2);
c=sqrt((x2d(elements2d(:,3))-x2d(elements2d(:,1))).^2+(y2d(elements2d(:,3))-y2d(elements2d(:,1))).^2);
% Compute area:
s=(a+b+c)/2;
Area=sqrt(s.*(s-a).*(s-b).*(s-c));
% Compute hydraulic diameter:
HydraulicDiameter_elements=2*Area./s;
% Compute element centers:
CentX=mean(x2d(elements2d),2);
CentY=mean(y2d(elements2d),2);
% Find mesh limits:
xlims=[min(x2d),max(x2d)];
ylims=[min(y2d),max(y2d)];
% Extract masks on 2D mesh and truncate:
OceanLevelSet_vertices=max(-1,min(1,md.mask.ocean_levelset(1:numvertices2d)));
IceLevelSet_vertices=max(-1,min(1,md.mask.ice_levelset(1:numvertices2d)));
% Smooth the ice mask:
if nummaskaveragings>0 && strcmp(overallmode,'setupandsend')
    IceLevelSet_vertices=averaging(md,repmat(IceLevelSet_vertices,[oldnumlayers,1]),nummaskaveragings);
    IceLevelSet_vertices(IceLevelSet_vertices<0)=-1;
    IceLevelSet_vertices(IceLevelSet_vertices>=0)=1;
    md.mask.ice_levelset=IceLevelSet_vertices;
end
% Compute ocean and ice area fractions:
OceanAreaFraction_elements=(1-mean(OceanLevelSet_vertices(elements2d),2))/2;
IceAreaFraction_elements=(1-mean(IceLevelSet_vertices(elements2d),2))/2;
% Compute area totals:
totalarea=sum(Area);
totalicearea=sum(Area.*IceAreaFraction_elements);
totalgroundedarea=sum(Area.*IceAreaFraction_elements.*(1-OceanAreaFraction_elements));
% Compute gradient of linear basis functions:
BasisGradX_elements=[(y2d(elements2d(:,2))-y2d(elements2d(:,3)))./(x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))...
    +x2d(elements2d(:,3)).*y2d(elements2d(:,1))+x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))),...
    (y2d(elements2d(:,3))-y2d(elements2d(:,1)))./(x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))...
    +x2d(elements2d(:,1)).*y2d(elements2d(:,2))+x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))),...
    (y2d(elements2d(:,1))-y2d(elements2d(:,2)))./(x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))...
    +x2d(elements2d(:,2)).*y2d(elements2d(:,3))+x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1)))];
BasisGradY_elements=[-(x2d(elements2d(:,2))-x2d(elements2d(:,3)))./(x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))...
    +x2d(elements2d(:,3)).*y2d(elements2d(:,1))+x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))),...
    -(x2d(elements2d(:,3))-x2d(elements2d(:,1)))./(x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))...
    +x2d(elements2d(:,1)).*y2d(elements2d(:,2))+x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))),...
    -(x2d(elements2d(:,1))-x2d(elements2d(:,2)))./(x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))...
    +x2d(elements2d(:,2)).*y2d(elements2d(:,3))+x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1)))];
% Interpolate hydraulic diameter to vertices:
% pre-allocate:
HydraulicDiameter_vertices=zeros(numvertices2d,1);
Denom_vertices=zeros(numvertices2d,1);
% loop through elements:
for jj=1:numelements2d
    % Add to the running sums:
    HydraulicDiameter_vertices(elements2d(jj,:))=HydraulicDiameter_vertices(elements2d(jj,:))+HydraulicDiameter_elements(jj);
    Denom_vertices(elements2d(jj,:))=Denom_vertices(elements2d(jj,:))+1;
end
% finish computing means:
HydraulicDiameter_vertices=HydraulicDiameter_vertices./Denom_vertices;


% Compute driving stress:
% Compute surface gradients:
SurfGradX_elements=sum(BasisGradX_elements.*md.geometry.surface(elements2d),2);
SurfGradY_elements=sum(BasisGradY_elements.*md.geometry.surface(elements2d),2);
% Compute gradient magnitude:
SurfGradMag_elements=sqrt(SurfGradX_elements.^2+SurfGradY_elements.^2);
% Average ice thickness onto the elements:
Icethick_elements=mean(md.geometry.thickness(elements2d),2);
% Compute driving stress:
DrivingStressX_elements=-md.materials.rho_ice*md.constants.g*SurfGradX_elements.*Icethick_elements;
DrivingStressY_elements=-md.materials.rho_ice*md.constants.g*SurfGradY_elements.*Icethick_elements;
DrivingStressMag_elements=md.materials.rho_ice*md.constants.g*SurfGradMag_elements.*Icethick_elements;
% Compute (area-weighted) stats:
% mean:
drivingstressmean=sum(DrivingStressMag_elements.*Area.*IceAreaFraction_elements)/totalicearea;
drivingstressgroundedmean=sum(DrivingStressMag_elements.*Area.*(1-OceanAreaFraction_elements).*IceAreaFraction_elements)/totalgroundedarea;
% stdev:
drivingstressstdev=sqrt(sum(((DrivingStressMag_elements-drivingstressmean).^2).*Area.*IceAreaFraction_elements)/totalicearea);
drivingstressgroundedstdev=sqrt(sum(((DrivingStressMag_elements-drivingstressgroundedmean).^2).*Area.*(1-OceanAreaFraction_elements).*IceAreaFraction_elements)/totalgroundedarea);
% Average driving stress from elements onto vertices:
% pre-allocate:
DrivingStressX_vertices=zeros(numvertices2d,1);
DrivingStressY_vertices=zeros(numvertices2d,1);
Denom_vertices=zeros(numvertices2d,1);
% loop through elements:
for jj=1:numelements2d
    % Add to the running sums:
    DrivingStressX_vertices(elements2d(jj,:))=DrivingStressX_vertices(elements2d(jj,:))+DrivingStressX_elements(jj);
    DrivingStressY_vertices(elements2d(jj,:))=DrivingStressY_vertices(elements2d(jj,:))+DrivingStressY_elements(jj);
    Denom_vertices(elements2d(jj,:))=Denom_vertices(elements2d(jj,:))+1;
end
% finish computing means:
DrivingStressX_vertices=DrivingStressX_vertices./Denom_vertices;
DrivingStressY_vertices=DrivingStressY_vertices./Denom_vertices;
% COmpute magnitude on the vertices:
DrivingStressMag_vertices=sqrt(DrivingStressX_vertices.^2+DrivingStressY_vertices.^2);

% Flatten old mesh:  (without erasing information)
if oldnumlayers>1
    md=FlattenModelStructure(md,'flatten');
end

%toc 
%return

%% Load, Process, and Interpolate Grids:

% Check if we even need to do any of this:
if strcmp(overallmode,'setupandsend')
    
    % Communicate:
    disp('Loading observed velocity grid...')
    
    % Load observed velocity:
    load(velfile,'X','Y',['U',velsuffix],['V',velsuffix])
    
    % Rename variables:
    X_obs=X;
    Y_obs=Y;
    if ~strcmp(velsuffix,'_obs')
        eval(['U_obs=U',velsuffix,';'])
        eval(['V_obs=V',velsuffix,';'])
    end
    
    % Clear extranesous variables:
    if ~strcmp(velsuffix,'_obs')
        clear('X','Y',['U',velsuffix],['V',velsuffix])
    else
        clear X Y
    end
    
    % Get grid spacing:
    dx_obs=X_obs(2)-X_obs(1);
    
    % Crop grids:
    ind1=find(Y_obs>=ylims(1)-3*dx_obs,1,'first');
    ind2=find(Y_obs<=ylims(2)+3*dx_obs,1,'last');
    ind3=find(X_obs>=xlims(1)-3*dx_obs,1,'first');
    ind4=find(X_obs<=xlims(2)+3*dx_obs,1,'last');
    X_obs=X_obs(ind3:ind4);
    Y_obs=Y_obs(ind1:ind2);
    U_obs=U_obs(ind1:ind2,ind3:ind4);
    V_obs=V_obs(ind1:ind2,ind3:ind4);
    
    % Get grid size:
    [ysize_obs,xsize_obs]=size(U_obs);
    
    % Identify good data:
    GoodData_obs=U_obs~=0&V_obs~=0&~isnan(U_obs)&~isnan(V_obs)&~isinf(U_obs)&~isinf(V_obs);
    U_obs(~GoodData_obs)=0;
    V_obs(~GoodData_obs)=0;
    
    % Add striped noise to the observations:
    if stripenoiseamp>0
        % Compute number of columns and rows:
        numcolumns=ceil(dx_obs*xsize_obs/stripenoiselength);
        numrows=ceil(dx_obs*ysize_obs/stripenoisewidth);
        % Loop through columns and rows:
        for d1=1:numrows
            for d2=1:numcolumns
                % Compute indices of grid cells in the stripe:
                ind1=1+ceil((d1-1)*stripenoisewidth/dx_obs);
                ind2=min(ysize_obs,ceil(d1*stripenoisewidth/dx_obs));
                ind3=1+ceil((d2-1)*stripenoiselength/dx_obs);
                ind4=min(xsize_obs,ceil(d2*stripenoiselength/dx_obs));
                %                     % Identify grid cells in the stripe:
                %                     InStripe=repmat((X_obs>=X_obs(1)+stripenoiselength*(d2-1))&(X_obs<X_obs(1)+stripenoiselength*d2),[ysize_obs,1])&...
                %                         repmat((Y_obs>=Y_obs(1)+stripenoisewidth*(d1-1))&(Y_obs<Y_obs(1)+stripenoisewidth*d1),[1,xsize_obs]);
                % Add correlated noise to velocity data:
                U_obs(ind1:ind2,ind3:ind4)=U_obs(ind1:ind2,ind3:ind4)+(sqrt(2)/2)*stripenoiseamp*randn(1);
                V_obs(ind1:ind2,ind3:ind4)=V_obs(ind1:ind2,ind3:ind4)+(sqrt(2)/2)*stripenoiseamp*randn(1);
            end
        end
    end
    
    % Add random noise to the observations:
    if randnoiseamp>0
        U_obs=U_obs+(sqrt(2)/2)*randnoiseamp*randn(size(U_obs));
        V_obs=V_obs+(sqrt(2)/2)*randnoiseamp*randn(size(V_obs));
    end
    
    % Reset velocity to zero in no data areas:
    if randnoiseamp>0 || stripenoiseamp>0
        U_obs(~GoodData_obs)=0;
        V_obs(~GoodData_obs)=0;
    end
    
    % Check if a drag coefficient guess file has been specified:
    if strcmp(dragcoeffguessstyle,'file')
        
        % Communicate:
        disp('Loading drag coefficient guess grid...')
        
        % Load drag coefficient file:
        loadnetcdf(dragcoefffile)
        
        % Reshape and rename drag coefficient:
        DragCoeff=C';
        clear C
        
        % Rename drag coeff coordinate variables:
        X_drag=X;
        Y_drag=Y;
        clear X Y
        
        % Get grid spacing:
        dx_drag=X_drag(2)-X_drag(1);
        
        % Crop drag grid:
        ind1=find(Y_drag>=ylims(1)-3*dx_drag,1,'first');
        ind2=find(Y_drag<=ylims(2)+3*dx_drag,1,'last');
        ind3=find(X_drag>=xlims(1)-3*dx_drag,1,'first');
        ind4=find(X_drag<=xlims(2)+3*dx_drag,1,'last');
        X_drag=X_drag(ind3:ind4);
        Y_drag=Y_drag(ind1:ind2);
        DragCoeff=DragCoeff(ind1:ind2,ind3:ind4);
        
        % Get grid size:
        [ysize_drag,xsize_drag]=size(DragCoeff);
        
        % Convert drag coefficient from m/yr to m/s:
        DragCoeff=DragCoeff*(secondsperyear^(1/m));
        
        % Take square root:
        DragCoeff=sqrt(DragCoeff);
        
        % Replace NaNs and INFs:
        DragCoeff(isnan(DragCoeff)|isinf(DragCoeff))=mean(DragCoeff(~isnan(DragCoeff)&~isinf(DragCoeff)));
        
    end
    
    % Communicate:
    if strcmp(dragcoeffguessstyle,'file')
        disp('Interpolating velocity and drag coefficient onto the mesh...')
    else
        disp('Interpolating velocity onto the mesh...')
    end
    
    % Extrapolate ice velocity into the ocean and the data gaps:
    U_obs(~GoodData_obs)=NaN;
    V_obs(~GoodData_obs)=NaN;
    U_obs=LaplacianInterpolation(U_obs,9);
    V_obs=LaplacianInterpolation(V_obs,9);
    
    % Interpolate grids onto the mesh:
    U_obs_vertices=MultiWavelengthInterpolator(X_obs,Y_obs,U_obs,x2d,y2d,elements2d);
    V_obs_vertices=MultiWavelengthInterpolator(X_obs,Y_obs,V_obs,x2d,y2d,elements2d);
    GoodData_obs_vertices=MultiWavelengthInterpolator(X_obs,Y_obs,double(GoodData_obs),x2d,y2d,elements2d);
    if strcmp(dragcoeffguessstyle,'file')
        DragCoeff_vertices=MultiWavelengthInterpolator(X_drag,Y_drag,DragCoeff,x2d,y2d,elements2d);
        DragCoeff_elements=mean(DragCoeff_vertices(elements2d),2);
    end
    
    % Replace remaining bad values with zero:
    BadValues=isnan(U_obs_vertices)|isnan(V_obs_vertices)|isnan(GoodData_obs_vertices);
    U_obs_vertices(BadValues)=0;
    V_obs_vertices(BadValues)=0;
    GoodData_obs_vertices(BadValues)=0;
    
    % Record good data mask in the model structure:
    md.miscellaneous.dummy=GoodData_obs_vertices;
    
    % Average velocity onto the elements:
    U_obs_elements=mean(U_obs_vertices(elements2d),2);
    V_obs_elements=mean(V_obs_vertices(elements2d),2);
    
    % Compute velocity magnitude:
    VelMag_obs_vertices=sqrt(U_obs_vertices.^2+V_obs_vertices.^2);
    VelMag_obs_elements=sqrt(U_obs_elements.^2+V_obs_elements.^2);
    
    
    
    % Get effective pressure:
    if strcmp(nsource,'oceanp')
        % Communicate:
        disp('Computing effective pressure from ocean water pressure (no cutoff)')
        % Compute water pressure (can be negative):
        WaterPressure_vertices=-md.materials.rho_water*md.constants.g*md.geometry.bed;
        % Compute ice pressure:
        IcePressure_vertices=md.materials.rho_ice*md.constants.g*md.geometry.thickness;
        % Compute effective pressure:
        EffectivePressure_vertices=max(IcePressure_vertices-WaterPressure_vertices,minn);
        % Average to the elements:
        EffectivePressure_elements=mean(EffectivePressure_vertices(elements2d),2);
    elseif strcmp(nsource,'oceanpwithcutoff')
        % Communicate:
        disp('Computing effective pressure from ocean water pressure (with cutoff)')
        % Compute water pressure (cannot be negative):
        WaterPressure_vertices=max(0,-md.materials.rho_water*md.constants.g*md.geometry.bed);
        % Compute ice pressure:
        IcePressure_vertices=md.materials.rho_ice*md.constants.g*md.geometry.thickness;
        % Compute effective pressure:
        EffectivePressure_vertices=max(IcePressure_vertices-WaterPressure_vertices,minn);
        % Average to the elements:
        EffectivePressure_elements=mean(EffectivePressure_vertices(elements2d),2);
    elseif strcmp(nsource,'bedelev')
        % Communicate:
        disp('Computing effective pressure from directly from bed elevation alone')
        % Compute effective pressure:
        EffectivePressure_vertices=max(minn+md.materials.rho_ice*md.constants.g*(md.geometry.bed-min(md.geometry.bed(IceLevelSet_vertices<0&OceanLevelSet_vertices>0))),minn);
        % Average to the elements:
        EffectivePressure_elements=mean(EffectivePressure_vertices(elements2d),2);
    elseif strcmp(nsource,'file')
        % Communicate:
        disp('Loading and interpolating effective pressure from a file')
        % NOTE:  you may need to change variable names here to match the
        % variable names in the file.
        % Load netcdf:
        loadnetcdf(nfile,'x','y','peffective','bnd_mask','effective_storativity','thk','topg')
        % Transpose variables to get in matlab standard orientation:
        peffective=peffective';
        bnd_mask=bnd_mask';
        effective_storativity=effective_storativity';
        x=x';
        % Check if the file had geometry variables:
        if exist('thk','var') && exist('topg','var') 
            % Transpose:
            thk=thk';
            topg=topg';
        else
            % Load geometry data from grid:
            GridGeom=load(gridfile,'X','Y','IceThick','BedElev');
            % Interpolate geometry to hydrology grid:
            thk=interp2(GridGeom.X,GridGeom.Y,GridGeom.IceThick,repmat(x,[length(y),1]),repmat(y,[1,length(x)]));
            topg=interp2(GridGeom.X,GridGeom.Y,GridGeom.BedElev,repmat(x,[length(y),1]),repmat(y,[1,length(x)]));
            clear GridGeom
        end
        % Fill shelf islands with Pw=Pocean:
        peffective(bnd_mask==2&topg<0)=md.materials.rho_ice*md.constants.g*thk(bnd_mask==2&topg<0)+md.materials.rho_water*md.constants.g*topg(bnd_mask==2&topg<0);
        % Identify grounded nunataks:
        isnunatak=(bnd_mask==2&topg>=0)|(bnd_mask==0&effective_storativity>.3999); % the threshold of 0.4 was suggested by Thomas Kleiner in an email
        % Identify halo around the nunataks:
        for thisiter=numnunatakiterations
            isnunatak=[false(1,length(x));false(length(y)-2,1),isnunatak(1:end-2,1:end-2)|isnunatak(2:end-1,1:end-2)|isnunatak(3:end,1:end-2)|...
                isnunatak(1:end-2,2:end-1)|isnunatak(2:end-1,2:end-1)|isnunatak(3:end,2:end-1)|...
                isnunatak(1:end-2,3:end)|isnunatak(2:end-1,3:end)|isnunatak(3:end,3:end),false(length(y)-2,1);false(1,length(x))];
        end
        % Smoothly interpolate over remaining invalid entries:
        peffective(isnunatak|bnd_mask==1|bnd_mask==4|peffective<0)=NaN;
        peffective=LaplacianInterpolation(peffective,9);
        % Interpolate onto the mesh:
        EffectivePressure_vertices=max(MultiWavelengthInterpolator(x,y,peffective,x2d,y2d,elements2d),minn);  % minimum applied after interpolation
        % Average to the elements:
        EffectivePressure_elements=mean(EffectivePressure_vertices(elements2d),2);
        % Clear variables:
        clear x y peffective bnd_mask topg thk haf isnunatak effective_storativity
    elseif strcmp(nsource,'none')
        % Dummy variable:
        EffectivePressure_elements=ones(numelements2d,1);
    else
        error('Parameter "nsource" must be a string equal to "none", "oceanp", "oceanpwithcutoff", "bedelev", or "file".')
    end
    
    % Smooth effective pressure:
    if numnaveragings>0 && ~strcmp(nsource,'none')
        EffectivePressure_vertices=averaging(md,EffectivePressure_vertices,numnaveragings);
        EffectivePressure_elements=mean(EffectivePressure_vertices(elements2d),2);
    end
    
    % Communicate:
    if ~strcmp(dragcoeffguessstyle,'file')
        disp('Computing drag coefficient guess...')
    end
    
    % Compute drag coefficient:
    if strcmp(dragcoeffguessstyle,'compute')
        % Resolve driving stress in the flow direction:
        ResolvedDrivingStress_elements=(U_obs_elements.*DrivingStressX_elements+V_obs_elements.*DrivingStressY_elements)./VelMag_obs_elements;
        % Compute (Weertman) drag coefficient on the elements:
        DragCoeff_elements=sqrt(max(ResolvedDrivingStress_elements,minreasonablestress_drag)./((max(VelMag_obs_elements,minreasonablevel)/secondsperyear).^(1/m)));
        % Modify drag coefficient for effective pressure:
        DragCoeff_elements=DragCoeff_elements./(EffectivePressure_elements.^(r/2));
        % Average from elements onto vertices:
        % pre-allocate:
        DragCoeff_vertices=zeros(numvertices2d,1);
        Denom_vertices=zeros(numvertices2d,1);
        % loop through elements:
        for ii=1:numelements2d
            % Add to the running sums:
            DragCoeff_vertices(elements2d(ii,:))=DragCoeff_vertices(elements2d(ii,:))+DragCoeff_elements(ii);
            Denom_vertices(elements2d(ii,:))=Denom_vertices(elements2d(ii,:))+1;
        end
        % finish computing means:
        DragCoeff_vertices=DragCoeff_vertices./Denom_vertices;
        DragCoeff_vertices=max(mindragcoeff,min(maxdragcoeff,DragCoeff_vertices));
    elseif strcmp(dragcoeffguessstyle,'previousinversion')
        % Use previous result for the drag coefficient:
        DragCoeff_vertices=md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d);
        % Average drag coefficient to the 2D elements:
        DragCoeff_elements=mean(DragCoeff_vertices(elements2d),2);
    elseif strcmp(dragcoeffguessstyle,'const')
        DragCoeff_vertices=dragcoeffguess*ones(numvertices2d,1);
        DragCoeff_elements=dragcoeffguess*ones(numelements2d,1);
    elseif strcmp(dragcoeffguessstyle,'previousspinup')
        % Use previous result for the drag coefficient:
        DragCoeff_vertices=md.friction.coefficient(1:numvertices2d);
        % Average drag coefficient to the 2D elements:
        DragCoeff_elements=mean(DragCoeff_vertices(elements2d),2);
    end
    
    % Apply smoothing to the guess:
    if dragguessaveragings>0
        DragCoeff_vertices=averaging(md,DragCoeff_vertices,dragguessaveragings);
        DragCoeff_elements=mean(DragCoeff_vertices(elements2d),2);
    end
    
    
    % Communicate:
    disp('Generating active masks...')
    
    % Compute active mask for drag:
    if includedrag
        ActiveMask_drag=IceLevelSet_vertices<=0&OceanLevelSet_vertices>=0;
    end
    
    % Compute active mask for rheology:
    if includerheology
        % Generate first Bbar mask:
        if rheologyfloatingonly
            ActiveMask_bbar=IceLevelSet_vertices<=0&OceanLevelSet_vertices<=0;
        else
            ActiveMask_bbar=IceLevelSet_vertices<=0;
        end
    end
    
    % Communicate:
    disp('Computing statistics...')
    
    % Compute (area-weighted) stats:
    if statsactiveareaonly
        % Get total active area on the elements:
        if includedrag && includerheology
            ActiveMask_elements=mean(ActiveMask_bbar(elements2d)|ActiveMask_drag(elements2d),2);
        elseif includedrag
            ActiveMask_elements=mean(ActiveMask_drag(elements2d),2);
        else
            ActiveMask_elements=mean(ActiveMask_bbar(elements2d),2);
        end
        % Ice thickness:
        icethickmean=sum(Icethick_elements.*Area.*ActiveMask_elements)/sum(Area.*ActiveMask_elements);
        % Velocity:
        % linear mean:
        velmean=sum(VelMag_obs_elements.*Area.*ActiveMask_elements)/sum(Area.*ActiveMask_elements);
        % linear stdev:
        velstdev=sqrt(sum(((VelMag_obs_elements-velmean).^2).*Area.*ActiveMask_elements)/sum(Area.*ActiveMask_elements));
        % RMS:
        velrms=sqrt(sum((VelMag_obs_elements.^2).*Area.*ActiveMask_elements)/sum(Area.*ActiveMask_elements));
        % Log mean:
        vellogmean=sum(log(VelMag_obs_elements(VelMag_obs_elements>0)).*Area(VelMag_obs_elements>0).*ActiveMask_elements(VelMag_obs_elements>0))/sum(Area(VelMag_obs_elements>0).*ActiveMask_elements(VelMag_obs_elements>0));
        % log stdev:
        vellogstdev=sqrt(sum(((log(VelMag_obs_elements(VelMag_obs_elements>0))-vellogmean).^2).*Area(VelMag_obs_elements>0).*ActiveMask_elements(VelMag_obs_elements>0))/sum(Area(VelMag_obs_elements>0).*ActiveMask_elements(VelMag_obs_elements>0)));
        % Median:
        [SortedVel,Ind]=sort(VelMag_obs_elements);
        CumulativeSortedArea=cumsum(Area(Ind).*ActiveMask_elements(Ind));
        velmedian=SortedVel(find(CumulativeSortedArea>=.5*sum(Area.*ActiveMask_elements),1,'first'));
        % Drag coefficient:
        if strcmp(dragcoeffguessstyle,'const')==0
            dragcoeffmean=sum(DragCoeff_elements.*Area.*ActiveMask_elements.*(1-OceanAreaFraction_elements).*IceAreaFraction_elements)/...
                sum(Area.*ActiveMask_elements.*(1-OceanAreaFraction_elements).*IceAreaFraction_elements);
            dragcoeffstdev=sqrt(sum(((DragCoeff_elements-dragcoeffmean).^2).*Area.*ActiveMask_elements.*(1-OceanAreaFraction_elements).*IceAreaFraction_elements)/...
                sum(Area.*ActiveMask_elements.*(1-OceanAreaFraction_elements).*IceAreaFraction_elements));
        else
            dragcoeffmean=dragcoeffguess;
            dragcoeffstdev=NaN;
        end
    else
        % Ice thickness:
        icethickmean=sum(Icethick_elements.*Area.*IceAreaFraction_elements)/totalicearea;
        % Velocity:
        % linear mean:
        velmean=sum(VelMag_obs_elements.*Area.*IceAreaFraction_elements)/totalicearea;
        % linear stdev:
        velstdev=sqrt(sum(((VelMag_obs_elements-velmean).^2).*Area.*IceAreaFraction_elements)/totalicearea);
        % RMS:
        velrms=sqrt(sum((VelMag_obs_elements.^2).*Area.*IceAreaFraction_elements)/totalicearea);
        % Log mean:
        vellogmean=sum(log(VelMag_obs_elements(VelMag_obs_elements>0)).*Area(VelMag_obs_elements>0).*IceAreaFraction_elements(VelMag_obs_elements>0))/...
            sum(Area(VelMag_obs_elements>0).*IceAreaFraction_elements(VelMag_obs_elements>0));
        % log stdev:
        vellogstdev=sqrt(sum(((log(VelMag_obs_elements(VelMag_obs_elements>0))-vellogmean).^2).*Area(VelMag_obs_elements>0).*IceAreaFraction_elements(VelMag_obs_elements>0))/...
            sum(Area(VelMag_obs_elements>0).*IceAreaFraction_elements(VelMag_obs_elements>0)));
        % Median:
        [SortedVel,Ind]=sort(VelMag_obs_elements);
        CumulativeSortedArea=cumsum(Area(Ind).*IceAreaFraction_elements(Ind));
        velmedian=SortedVel(find(CumulativeSortedArea>=.5*totalicearea,1,'first'));
        % Drag coefficient:
        if strcmp(dragcoeffguessstyle,'const')==0
            dragcoeffmean=sum(DragCoeff_elements.*Area.*(1-OceanAreaFraction_elements).*IceAreaFraction_elements)/totalgroundedarea;
            dragcoeffstdev=sqrt(sum(((DragCoeff_elements-dragcoeffmean).^2).*Area.*(1-OceanAreaFraction_elements).*IceAreaFraction_elements)/totalgroundedarea);
        else
            dragcoeffmean=dragcoeffguess;
            dragcoeffstdev=NaN;
        end
    end
    % Record stats:
    MyInversionParameters.icethickmean=icethickmean;
    MyInversionParameters.velmean=velmean;
    MyInversionParameters.velstdev=velstdev;
    MyInversionParameters.velrms=velrms;
    MyInversionParameters.vellogmean=vellogmean;
    MyInversionParameters.vellogstdev=vellogstdev;
    MyInversionParameters.velmedian=velmedian;
    MyInversionParameters.dragcoeffmean=dragcoeffmean;
    MyInversionParameters.dragcoeffstdev=dragcoeffstdev;
    
    
else
    
    % Unpack the good data mask:
    GoodData_obs_vertices=md.miscellaneous.dummy;
    
    % Compute active mask for drag:
    if includedrag
        ActiveMask_drag=IceLevelSet_vertices<=0&OceanLevelSet_vertices>=0;
    end
    
    % Compute active mask for rheology:
    if includerheology
        if rheologyfloatingonly
            ActiveMask_bbar=IceLevelSet_vertices<=0&OceanLevelSet_vertices<=0;
        else
            ActiveMask_bbar=IceLevelSet_vertices<=0;
        end
    end
    
end


%% Set Up Model Structure:

% Check if we even need to do any of this:
if strcmp(overallmode,'setupandsend')
    
    % Communicate:
    disp('Setting up model structure...')
    
    % Check whether to extrude or collapse the mesh:
    if numlayers~=oldnumlayers
        % Collapse the old mesh:
        if oldnumlayers>1
            md=FlattenModelStructure(md,'expand');
            md=collapse(md);
        end
        % Extrude and flatten the new mesh:
        if numlayers>1
            md=extrude(md,numlayers,densifyfactor);
            md=FlattenModelStructure(md,'flatten');
        end
    end
    
    % Get mesh properties:
    numelements=md.mesh.numberofelements;
    numvertices=md.mesh.numberofvertices;
    
    % Set Dirichlet BC and initial velocity guess:
    md.stressbalance.spcvx=NaN*zeros(numvertices,1);
    md.stressbalance.spcvy=NaN*zeros(numvertices,1);
    md.stressbalance.spcvz=NaN*zeros(numvertices,1);
    if inlandbcpower>0 && numlayers>1 && ~strcmp(flowequation,'SSA')
        % Compute Zhat:
        Zhat=(md.mesh.z(1:numvertices2d:end)-md.mesh.z(1))/md.geometry.thickness(1);
        % Compute shape function:
        ShapeFunction=1-(1-Zhat).^inlandbcpower;
        % Set BC:
        md.stressbalance.spcvx(repmat(md.mesh.vertexonboundary==1&IceLevelSet_vertices<0,[numlayers,1]))=...
            repmat(U_obs_vertices(md.mesh.vertexonboundary==1&IceLevelSet_vertices<0),[numlayers,1]).*...
            reshape(repmat(ShapeFunction',[sum(md.mesh.vertexonboundary==1&IceLevelSet_vertices<0),1]),[],1);
        md.stressbalance.spcvy(repmat(md.mesh.vertexonboundary==1&IceLevelSet_vertices<0,[numlayers,1]))=...
            repmat(V_obs_vertices(md.mesh.vertexonboundary==1&IceLevelSet_vertices<0),[numlayers,1]).*...
            reshape(repmat(ShapeFunction',[sum(md.mesh.vertexonboundary==1&IceLevelSet_vertices<0),1]),[],1);
        % Set initialization velocity:
        md.initialization.vx=repmat(U_obs_vertices,[numlayers,1]).*reshape(repmat(1-.5*(1-ShapeFunction'),[numvertices2d,1]),[],1); % slip fraction=0.5 in first guess
        md.initialization.vy=repmat(V_obs_vertices,[numlayers,1]).*reshape(repmat(1-.5*(1-ShapeFunction'),[numvertices2d,1]),[],1);
        md.initialization.vel=repmat(VelMag_obs_vertices,[numlayers,1]).*reshape(repmat(1-.5*(1-ShapeFunction'),[numvertices2d,1]),[],1);
    else
        % Set BC:
        md.stressbalance.spcvx(repmat(md.mesh.vertexonboundary==1&IceLevelSet_vertices<0,[numlayers,1]))=repmat(U_obs_vertices(md.mesh.vertexonboundary==1&IceLevelSet_vertices<0),[numlayers,1]);
        md.stressbalance.spcvy(repmat(md.mesh.vertexonboundary==1&IceLevelSet_vertices<0,[numlayers,1]))=repmat(V_obs_vertices(md.mesh.vertexonboundary==1&IceLevelSet_vertices<0),[numlayers,1]);
        % Set initialization velocity:
        md.initialization.vx=repmat(U_obs_vertices,[numlayers,1]);
        md.initialization.vy=repmat(V_obs_vertices,[numlayers,1]);
        md.initialization.vel=repmat(VelMag_obs_vertices,[numlayers,1]);
    end
    
    % Selectively expand geometry variables:
    md.geometry.thickness=repmat(md.geometry.thickness,[numlayers,1]);
    md.geometry.surface=repmat(md.geometry.surface,[numlayers,1]);
    md.geometry.base=repmat(md.geometry.base,[numlayers,1]);
    md.geometry.bed=repmat(md.geometry.bed,[numlayers,1]);
    md.mesh.x=repmat(md.mesh.x,[numlayers,1]);
    md.mesh.y=repmat(md.mesh.y,[numlayers,1]);
    
    % Check temperature source:
    if strcmp(tempsource,'file')
        % Interpolate with my new script:
        md=pism2issm_mike(md,tempfile,minpismicethick,numtempsmoothings);
        % Set any NaNs to default ice temp:
        md.initialization.temperature(isnan(md.initialization.temperature))=defaulticetemp;
        md.initialization.waterfraction(isnan(md.initialization.waterfraction))=0;
    elseif strcmp(tempsource,'1D')
        % Load boundary conditions from grid file:
        load(gridfile,'X','Y','SurfTemp','Accum','Gflux_martos')
        % Interpolate BC onto model mesh:
        SurfTemp_vertices=interp2(X,Y,SurfTemp,x2d,y2d);
        Accum_vertices=interp2(X,Y,Accum,x2d,y2d);
        Gflux_vertices=interp2(X,Y,Gflux_martos,x2d,y2d);
        % Fill NaNs in BC:
        if sum(isnan(SurfTemp_vertices))>0
            SurfTemp_vertices=LaplacianInterpolation_mesh(x2d,y2d,elements2d,SurfTemp_vertices);
        end
        if sum(isnan(Accum_vertices))>0
            Accum_vertices=LaplacianInterpolation_mesh(x2d,y2d,elements2d,Accum_vertices);
        end
        if sum(isnan(Gflux_vertices))>0
            Gflux_vertices=LaplacianInterpolation_mesh(x2d,y2d,elements2d,Gflux_vertices);
        end
        % Shift surface temperature at low elevations:
        SurfTemp_vertices=SurfTemp_vertices+surftempshift*exp(-md.geometry.surface(1:numvertices2d)/surftempshiftz0);
        % Estimate shear heating:
        ShearHeat_vertices=max(0,DrivingStressX_vertices.*U_obs_vertices+DrivingStressY_vertices.*V_obs_vertices)/secondsperyear;
        % Smooth out shear heating:
        ShearHeat_vertices=averaging(md,repmat(ShearHeat_vertices,[numlayers,1]),numshearheataveragings,1);
        ShearHeat_vertices=ShearHeat_vertices(1:numvertices2d);
        % Create structures for the 1D column function:
        BoundaryConditions=struct('surftemp',0,'accum',0,'icethick',0,'gflux',0,'meltrate',0,'g',md.constants.g);
        MaterialParameters=struct('cond_i',md.materials.thermalconductivity,'rho_i',md.materials.rho_ice,'specheat_i',md.materials.heatcapacity,'latentheat',md.materials.latentheat,'beta',md.materials.beta);
        % Pre-allocate temperature:
        md.initialization.temperature=zeros(numvertices2d*numlayers,1);
        % Loop through vertices:
        for ii=1:numvertices2d
            % Check if we have ice here:
            if IceLevelSet_vertices(ii)>0
                % Fill in with default temp:
                md.initialization.temperature(ii:numvertices2d:end)=defaulticetemp;
                continue
            end
            % Fill in boundary conditions structure:
            BoundaryConditions.surftemp=SurfTemp_vertices(ii);
            BoundaryConditions.accum=Accum_vertices(ii);
            BoundaryConditions.gflux=Gflux_vertices(ii)+ShearHeat_vertices(ii);
            BoundaryConditions.icethick=md.geometry.thickness(ii);
            % Check if this is a floating shelf:
            if OceanLevelSet_vertices(ii)<=0
                % Compute tied temperature solution:
                [ThisDepth,ThisTemp,~,~]=TempProfile1D_simple_v2(BoundaryConditions,MaterialParameters,2,3*numlayers);
            else
                % Try untied temperature solution:
                [ThisDepth,ThisTemp,~,~]=TempProfile1D_simple_v2(BoundaryConditions,MaterialParameters,1,3*numlayers);
                % Check if basal temperature exceeds the melting point:
                if ThisTemp(end)>-ThisDepth(end)*MaterialParameters.rho_i*BoundaryConditions.g*MaterialParameters.beta
                    % Compute tied temperature solution:
                    [ThisDepth,ThisTemp,~,~]=TempProfile1D_simple_v2(BoundaryConditions,MaterialParameters,2,3*numlayers);
                end
            end
            % Get depth for this column:
            ThisISSMDepth=max(0,min(md.geometry.thickness(ii),md.geometry.surface(ii)-md.mesh.z(ii:numvertices2d:end)));
            % Interpolate temperature onto this column:
            md.initialization.temperature(ii:numvertices2d:end)=md.materials.meltingpoint+interp1(ThisDepth,ThisTemp,ThisISSMDepth); % also converts from C to K
        end
        % Set zero water fraction:
        md.initialization.waterfraction=zeros(numvertices,1);
    elseif strcmp(tempsource,'const')
        % Use default temperature:
        md.initialization.temperature=defaulticetemp*ones(numvertices,1);
        md.initialization.waterfraction=zeros(numvertices,1);
    else % md
        % Check that existing model temperature field is legit:
        if length(md.initialization.temperature)~=numvertices
            error('md.initialization.temperature must be length numvertices')
        end
        if length(md.initialization.waterfraction)~=numvertices
            error('md.initialization.waterfraction must be length numvertices')
        end
        if sum(isnan(md.initialization.temperature))>0 || sum(isinf(md.initialization.temperature))>0
            error('md.initialization.temperature contains NaN or Inf')
        end
        if sum(isnan(md.initialization.waterfraction))>0 || sum(isinf(md.initialization.waterfraction))>0
            error('md.initialization.waterfraction contains NaN or Inf')
        end
        if min(md.initialization.temperature)<md.materials.meltingpoint-100 || max(md.initialization.temperature)>md.materials.meltingpoint+1
            error('md.initialization.temperature contains values out of bounds')
        end
        if min(md.initialization.waterfraction)<0 || max(md.initialization.waterfraction)>1
            error('md.initialization.waterfraction contains values out of bounds')
        end
    end
    
    % Set rheology:
    md.materials.rheology_n=n*ones(numelements,1);  % rheology can be vertically variable
    if strcmp(rheologystyle,'cuffeytemperate')
        if includewaterfraction
            md.materials.rheology_B=cuffeytemperate(md.initialization.temperature,md.initialization.waterfraction,n);
        else
            md.materials.rheology_B=cuffeytemperate(md.initialization.temperature,zeros(numvertices,1),n);
        end
        % Set invalid entries to default ice temperature:
        md.materials.rheology_B(md.materials.rheology_B<=0|isnan(md.materials.rheology_B)|isinf(md.materials.rheology_B))=cuffeytemperate(defaulticetemp,0,n);
    else
        md.materials.rheology_B=paterson(md.initialization.temperature);
        % Set invalid entries to default ice temperature:
        md.materials.rheology_B(md.materials.rheology_B<=0|isnan(md.materials.rheology_B)|isinf(md.materials.rheology_B))=paterson(defaulticetemp);
    end
    
    % Modify rheology to try to prevent deformation velocity greater than
    % observed velocity:
    if reducedefvel && ~strcmp(flowequation,'SSA')
        % Compute zhat:
        Zhat=(md.mesh.z(1:numvertices2d:end)-md.geometry.base(1))/md.geometry.thickness(1);
        % Compute shear stress guess from smoothing of driving stress:
        if numudefstressaveragings>0
            ShearStressXGuess=averaging(md,repmat(DrivingStressX_vertices,[numlayers,1]),numudefstressaveragings);
            ShearStressYGuess=averaging(md,repmat(DrivingStressY_vertices,[numlayers,1]),numudefstressaveragings);
        else
            ShearStressXGuess=repmat(DrivingStressX_vertices,[numlayers,1]);
            ShearStressYGuess=repmat(DrivingStressY_vertices,[numlayers,1]);
        end
        % Resolve shear stress in the flow direction:
        ShearStressFlowDirGuess=max(0,(repmat(U_obs_vertices,[numlayers,1]).*ShearStressXGuess+repmat(V_obs_vertices,[numlayers,1]).*ShearStressYGuess)./repmat(VelMag_obs_vertices+velstabilizer_def,[numlayers,1]));
        % Modify shear stress guess to a linear profile between bed and
        % surface:
        ShearStressFlowDirGuess=ShearStressFlowDirGuess.*reshape(repmat(1-Zhat',[numvertices2d,1]),[],1);
        % Compute rheology A coefficient:
        A=md.materials.rheology_B.^(-n);
        % Compute SIA deformation velocity:
        Udef=secondsperyear*md.geometry.thickness(1:numvertices2d).*trapz(Zhat,reshape(A.*(ShearStressFlowDirGuess.^n),[numvertices2d,numlayers])')';  % uses trapezoidal integration
        % Compute ratio of observed vel to deformation vel:
        Ratio=Udef./(velstabilizer_def+VelMag_obs_vertices);
        % Apply mask bounds:
        Ratio(IceLevelSet_vertices>0|OceanLevelSet_vertices<0)=1;
        % Smooth the ratio:
        if numratioaveragings>0
            Ratio=averaging(md,repmat(Ratio,[numlayers,1]),numratioaveragings); 
        else
            Ratio=repmat(Ratio,[numlayers,1]);
        end
        % Modify rheology to make the ice stiffer where necessary:
        md.materials.rheology_B=md.materials.rheology_B.*(max(1,Ratio).^(1/n));
    end
    
    % Compute vertically averaged rheology:
    InitialBbar_vertices=DepthAverage(md,md.materials.rheology_B);
    
    % Smoothly fill non-ice areas:
    InitialBbar_vertices(IceLevelSet_vertices>0)=NaN;
    InitialBbar_vertices=LaplacianInterpolation_mesh(x2d,y2d,elements2d,InitialBbar_vertices);
    
    % Collapse model structure if in SSA:
    if strcmp(flowequation,'SSA')
        md=FlattenModelStructure(md,'expand');
        md=collapse(md);
        md.materials.rheology_B=InitialBbar_vertices;
        numelements=md.mesh.numberofelements;
        numvertices=md.mesh.numberofvertices;
    end
    
    % Get rheology stats:
    if includerheology
        % Average to elements:
        InitialBbar_elements=mean(InitialBbar_vertices(elements2d),2);
        % Get stats:
        if statsactiveareaonly
            ActiveMask_bar_elements=mean(ActiveMask_bbar(elements2d),2);
            bbarmean=sum(ActiveMask_bar_elements.*InitialBbar_elements.*Area.*IceAreaFraction_elements)/sum(Area.*ActiveMask_bar_elements.*IceAreaFraction_elements);
            bbarstdev=sqrt(sum(((InitialBbar_elements-bbarmean).^2).*Area.*ActiveMask_bar_elements.*IceAreaFraction_elements)/sum(Area.*ActiveMask_bar_elements.*IceAreaFraction_elements));
        else
            bbarmean=sum(InitialBbar_elements.*Area.*IceAreaFraction_elements)/totalicearea;
            bbarstdev=sqrt(sum(((InitialBbar_elements-bbarmean).^2).*Area.*IceAreaFraction_elements)/totalicearea);
        end
        % Record stats:
        MyInversionParameters.bbarmean=bbarmean;
        MyInversionParameters.bbarstdev=bbarstdev;
        % Replace vertically resolved rheology with initial guess of Bbar:
        if numlayers>1 && ~strcmp(flowequation,'SSA')
            md.materials.rheology_B=repmat(InitialBbar_vertices,[numlayers,1]);
        end
    end
    
    % Un-expand geometry:
    md.geometry.thickness=md.geometry.thickness(1:numvertices2d);
    md.geometry.surface=md.geometry.surface(1:numvertices2d);
    md.geometry.base=md.geometry.base(1:numvertices2d);
    md.geometry.bed=md.geometry.bed(1:numvertices2d);
    md.mesh.x=md.mesh.x(1:numvertices2d);
    md.mesh.y=md.mesh.y(1:numvertices2d);
    
    % Set the flow equation:
    md=setflowequation(md,flowequation,'all');
    
    % Set loading force and referential:
    %md.stressbalance.loadingforce=[zeros(numvertices,2),md.constants.g*ones(numvertices,1)];
    md.stressbalance.loadingforce=[zeros(numvertices,2),md.constants.g*md.materials.rho_ice*ones(numvertices,1)];
    md.stressbalance.referential=NaN*zeros(numvertices,6);
    
    % Set name:
    slashnum=strfind(outputfile,'/');
    md.miscellaneous.name=outputfile(slashnum(end)+1:end-4);
    
    % Set basal forcings:
    md.basalforcings.groundedice_melting_rate=zeros(numvertices2d,1);
    md.basalforcings.floatingice_melting_rate=zeros(numvertices2d,1);
    
    % Initialize inversion structure:
    md.inversion=m1qn3inversion();
    md.inversion.iscontrol=1;
    md.inversion.incomplete_adjoint=useincompleteadjoint;
    md.inversion.maxsteps=numsteps;
    md.inversion.maxiter=numsteps;
    md.inversion.dxmin=dxmin;
    md.inversion.gttol=gttol;
    
    % Assign velocity to inversion structure:
    md.inversion.vx_obs=U_obs_vertices;
    md.inversion.vy_obs=V_obs_vertices;
    md.inversion.vel_obs=VelMag_obs_vertices;
    
    % Set up friction structure:
    md.friction.coefficient=DragCoeff_vertices;
    md.friction.p=m*ones(numelements2d,1);
    if ~strcmp(nsource,'none')
        md.friction.q=r*m*ones(numelements2d,1);
        md.friction.effective_pressure=EffectivePressure_vertices;
        md.friction.coupling=3;
    else
        md.friction.q=zeros(numelements2d,1);
        md.friction.effective_pressure=zeros(numvertices2d,1);
        md.friction.coupling=0;
    end
    
    % What type of inversion?
    if includedrag && includerheology
        % Assign control parameters:
        md.inversion.control_parameters={'FrictionCoefficient','MaterialsRheologyBbar'};
        % Compute total area for normalizing cost functions:
        thisarea=totalicearea;
    elseif includedrag
        % Assign control parameters:
        md.inversion.control_parameters={'FrictionCoefficient'};
        % Compute total area for normalizing cost functions:
        thisarea=totalgroundedarea;
    elseif includerheology
        % Assign control parameters:
        md.inversion.control_parameters={'MaterialsRheologyBbar'};
        % Compute total area for normalizing cost functions:
        thisarea=totalicearea-totalgroundedarea;
    else
        error('So if you dont want me to include drag, and you dont want me to include rheology, what the fuck kind of inversion do you want me to run here?')
    end
    
    % Set up cost functions:
    % Assign codes:
    md.inversion.cost_functions=costfunctions;
    % Pre-allocate coefficients:
    md.inversion.cost_functions_coefficients=ones(numvertices2d,length(costfunctions));
    % Loop through each component:
    for ii=1:length(costfunctions)
        % Check component and compute characteristic scale:
        if costfunctions(ii)==101 % vel absolute misfit
            if strcmp(velabsscale,'mean')
                thisscale=.5*(velmean/secondsperyear)^2;
            elseif strcmp(velabsscale,'stdev')
                thisscale=.5*(velstdev/secondsperyear)^2;
            elseif strcmp(velabsscale,'rms')
                thisscale=.5*(velrms/secondsperyear)^2;
            elseif strcmp(velabsscale,'median')
                thisscale=.5*(velmedian/secondsperyear)^2;
            else
                thisscale=.5*(velabsscale/secondsperyear)^2;
            end
        elseif costfunctions(ii)==102  % vel relative misfit
            thisscale=.5;
            % Apply stabilizer:
            thisscale=thisscale*(VelMag_obs_vertices+velstabilizer)./VelMag_obs_vertices;
        elseif costfunctions(ii)==102.1  % homebrew vel relative misfit
            % Set the internal cost function code to abslute velocity misfit:
            md.inversion.cost_functions(ii)=101;
            % Construct a scaling coefficient that makes the absolute misfit
            % act like a relative misfit:
            thisscale=.5*((VelMag_obs_vertices/secondsperyear).^2+(velstabilizer/secondsperyear)^2);
        elseif costfunctions(ii)==103 % vel log misfit
            if strcmp(vellogscale,'logstdev')
                thisscale=4*((issmmeanvel/secondsperyear)^2)*(vellogstdev^2);
            elseif strcmp(vellogscale,'1+stdev/mean')
                thisscale=4*((issmmeanvel/secondsperyear)^2)*(log(1+velstdev/velmean)^2);
            elseif strcmp(vellogscale,'1+rms/mean')
                thisscale=4*((issmmeanvel/secondsperyear)^2)*(log(1+velrms/velmean)^2);
            else
                thisscale=4*((issmmeanvel/secondsperyear)^2)*(log(vellogscale)^2);
            end
            % Note that all of these account for an arbitrary scaling
            % factor used internally in ISSM, 4*((issmmeanvel/secondsperyear)^2)
            % Apply stabilizer:
            thisscale=thisscale*(VelMag_obs_vertices+velstabilizer)./VelMag_obs_vertices;
        elseif costfunctions(ii)==501 % drag coeff gradient misfit
            % Use ice thickness as the wavelength:
            thiswavelength=icethickmean;
            % Check if we should use mean or stdev for scale:
            if strcmp(dragcoeffguessstyle,'const')==0
                thisscale=(pi*dragcoeffstdev/thiswavelength)^2;
            else
                thisscale=(pi*dragcoeffmean/thiswavelength)^2;
            end
        elseif costfunctions(ii)==502 % Bbar gradient misfit
            % Use ice thickness as the wavelength:
            thiswavelength=icethickmean;
            % Check if we should use mean or stdev for scale:
            if ~strcmp(tempsource,'const')
                thisscale=(pi*bbarstdev/thiswavelength)^2;
            else
                thisscale=(pi*bbarmean/thiswavelength)^2;
            end
        else
            error('I only know how to compute the characteristic scales for cost function components 101, 102, 102.1, 103, 501, and 502.')
        end
        % Compute coefficients:
        md.inversion.cost_functions_coefficients(:,ii)=intrinsicweights(ii)./(thisscale*thisarea);
        % Modify coefficient for bad data:
        if costfunctions(ii)<500
            md.inversion.cost_functions_coefficients(:,ii)=md.inversion.cost_functions_coefficients(:,ii).*(GoodData_obs_vertices.^2);
        end
        % Apply active masks:
        if costfunctions(ii)==501
            % Drag regularization:
            md.inversion.cost_functions_coefficients(:,ii)=md.inversion.cost_functions_coefficients(:,ii).*ActiveMask_drag;
        elseif costfunctions(ii)==502
            % Rheology regularization:
            md.inversion.cost_functions_coefficients(:,ii)=md.inversion.cost_functions_coefficients(:,ii).*ActiveMask_bbar;
        else
            % Data regularization, must check type of inversion:
            if includedrag && includerheology
                md.inversion.cost_functions_coefficients(:,ii)=md.inversion.cost_functions_coefficients(:,ii).*(ActiveMask_drag|ActiveMask_bbar);
            elseif includedrag
                md.inversion.cost_functions_coefficients(:,ii)=md.inversion.cost_functions_coefficients(:,ii).*ActiveMask_drag;
            else
                md.inversion.cost_functions_coefficients(:,ii)=md.inversion.cost_functions_coefficients(:,ii).*ActiveMask_bbar;
            end
        end
        % Remove non-ice regions from the cost function:
        md.inversion.cost_functions_coefficients(:,ii)=md.inversion.cost_functions_coefficients(:,ii).*(IceLevelSet_vertices<=0);
        % Remove inland boundary from the cost function (for data):
        if costfunctions(ii)<500
            md.inversion.cost_functions_coefficients(:,ii)=md.inversion.cost_functions_coefficients(:,ii).*(~(md.mesh.vertexonboundary&IceLevelSet_vertices<0));
        end
    end
    
    % Assign max/min limits:
    if includedrag && includerheology
        % Constant limits:
        md.inversion.min_parameters=[mindragcoeff*ones(numvertices2d,1),minbbar*ones(numvertices2d,1)];
        md.inversion.max_parameters=[maxdragcoeff*ones(numvertices2d,1),maxbbar*ones(numvertices2d,1)];
        % Modify for active regions (drag):
        md.inversion.min_parameters(ActiveMask_drag==0,1)=0;
        md.inversion.max_parameters(ActiveMask_drag==0,1)=0;
        % Modify for active regions (rheology):
        md.inversion.min_parameters(ActiveMask_bbar==0,2)=md.materials.rheology_B(ActiveMask_bbar==0);
        md.inversion.max_parameters(ActiveMask_bbar==0,2)=md.materials.rheology_B(ActiveMask_bbar==0);
    elseif includedrag
        % Constant limits:
        md.inversion.min_parameters=mindragcoeff*ones(numvertices2d,1);
        md.inversion.max_parameters=maxdragcoeff*ones(numvertices2d,1);
        % Modify for active regions:
        md.inversion.min_parameters(ActiveMask_drag==0)=0;
        md.inversion.max_parameters(ActiveMask_drag==0)=0;
    elseif includerheology
        % Constant limits:
        md.inversion.min_parameters=minbbar*ones(numvertices2d,1);
        md.inversion.max_parameters=maxbbar*ones(numvertices2d,1);
        % Modify for active regions:
        md.inversion.min_parameters(ActiveMask_bbar==0)=md.materials.rheology_B(ActiveMask_bbar==0);
        md.inversion.max_parameters(ActiveMask_bbar==0)=md.materials.rheology_B(ActiveMask_bbar==0);
    end
    
    % Compute control scaling factors:
    if includedrag && includerheology
        md.inversion.control_scaling_factors=[dragcoeffmean,bbarmean];
    else
        md.inversion.control_scaling_factors=1;
        %     elseif includedrag
        %         md.inversion.control_scaling_factors=dragcoeffmean;
        %     elseif includerheology
        %         md.inversion.control_scaling_factors=bbarmean;
    end
    
    % Set solvers:
    md.toolkits=toolkits;
    if strcmp(solvertype,'ilua')
        md.toolkits=addoptions(md.toolkits,'DefaultAnalysis',iluasmoptions());
        md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis',iluasmoptions());
    elseif strcmp(solvertype,'bjacobi')
        md.toolkits=addoptions(md.toolkits,'DefaultAnalysis',gmresbjacobioptions());
        md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis',gmresbjacobioptions());
    elseif strcmp(solvertype,'asmstokes')
        md.toolkits=addoptions(md.toolkits,'StressbalanceAnalysis',asmstokesoptions());
        %md.toolkits=addoptions(md.toolkits,'DefaultAnalysis',asmstokesoptions());
        md.toolkits=addoptions(md.toolkits,'DefaultAnalysis',iluasmoptions());  % asmstokes only used for stress balance!
    else
        error('Parameter "solvertype" not recognized')
    end
    md.settings.solver_residue_threshold=solver_residue_threshold;
    md.settings.lowmem=lowmem;
    md.stressbalance.restol=restol;
    md.stressbalance.reltol=reltol;
    md.stressbalance.abstol=abstol;
    md.stressbalance.maxiter=maxiter;
    
    % Control verbosity:
    md.verbose=verbose('control',true,'convergence',true);
    
    % Control waitonlock:
    md.settings.waitonlock=0;
    
    % Compute number of nodes and tasks:
    %numtasks=max(1,floor(numelements/numelementspertask));
    %numnodes=ceil(numtasks/numtaskspernode);
    numnodes=ceil(numelements/(numtaskspernode*numelementspertask));
    numtasks=numnodes*numtaskspernode;
    
    % Set up cluster info:
    md.cluster=ollie('nodes',numnodes,'ntasks',numtasks,'ntaskspernode',numtaskspernode,'time',rtime);
    
    % Set to local cluster:
    %md.cluster=local('name',oshostname(),'np',3);
    %md.cluster.np=3;
    
    % Remove the transient solution from md:
    if isfield(md.results,'TransientSolution')
        md.results=rmfield(md.results,'TransientSolution');
    end
    
    % Remove the old stress balance solution from md:
    if isfield(md.results,'StressbalanceSolution')
        md.results=rmfield(md.results,'StressbalanceSolution');
    end
    
end


%% Solve Model:

% Check overall mode:
if strcmp(overallmode,'setupandsend')
    
    % Check if we're doing an l-curve analysis:
    if dolcurve
        
        % Communicate:
        disp('Running L-curve analysis...')
        
        % Make lambda sampling:  (log scale)
        Lambda=exp(linspace(log(lambdamin),log(lambdamax),numlambda)');
        
        % Identify which elements of the cost function are regularization:
        AreReg=costfunctions>500;
        
        % Save basic model structure:
        BasicModel=md;
        
        % Communicate:
        disp('...saving pre-solution model structure')
        
        % Save pre-solution model:
        hassolved=0;
        save(outputfile,'BasicModel','hassolved','Lambda','*Parameters','-v7.3')
        
        % Communicate:
        disp('...looping through lambda values')
        
        % Loop through lambda values:
        for ii=1:numlambda
            
            % Communicate:
            disp('...')
            disp(['......lambda=',num2str(Lambda(ii))])
            
            % Construct this file name:
            slashnum=strfind(outputfile,'/');
            thisbase=10^floor(log10(Lambda(ii)));
            thisname=[outputfile(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))];
            
            % Make this model structure:
            md=BasicModel;
            
            % Apply lambda values:
            md.inversion.cost_functions_coefficients(repmat(AreReg,[numvertices2d,1]))=md.inversion.cost_functions_coefficients(repmat(AreReg,[numvertices2d,1]))*Lambda(ii);
            
            % Assign name:
            md.miscellaneous.name=thisname;
            
            % Save pre-solution model:
            disp('......saving pre-solution model structure')
            hassolved=0;
            save([outputfile(1:slashnum(end)),thisname,'.mat'],'md','hassolved','MyInversionParameters','-v7.3')
            
            % Solve model:
            disp('......sending model to cluster')
            md=FlattenModelStructure(md,'expand');
            md=solve(md,'sb');
            md=FlattenModelStructure(md,'flatten');
            
        end
        
        % End here:
        disp('...')
        disp('All models submitted to server, my work here is done!')
        toc
        return
        
    else
        
        % Communicate:
        disp('Running Inversion...')
        
        % Save pre-solution model:
        disp('...saving pre-solution model structure')
        hassolved=0;
        save(outputfile,'md','hassolved','MyInversionParameters','-v7.3')
        
        % Communicate:
        disp('...sending model to cluster')
        
        % Solve model:
        md=FlattenModelStructure(md,'expand');
        md=solve(md,'sb');
        md=FlattenModelStructure(md,'flatten');
        
        % End here:
        disp('...')
        disp('Model submitted to server, my work here is done!')
        toc
        return
        
    end
    
else
    
    % Check if we're doing an L-Curve analysis:
    if dolcurve
        
        % Check if it's solved:
        if ~hassolved
            
            % Communicate:
            disp('Loading L-curve results...')
            
            % Pre-allocate cost values:
            DataCost=zeros(numlambda,1);
            RegularizationCost=zeros(numlambda,1);
            
            % Loop through lambda values:
            for ii=1:numlambda
                
                % Communicate:
                disp('...')
                disp(['...lambda=',num2str(Lambda(ii))])
                
                % Construct this file name:
                slashnum=strfind(outputfile,'/');
                thisbase=10^floor(log10(Lambda(ii)));
                thisname=[outputfile(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))];
                
                % Load this model:
                load([outputfile(1:slashnum(end)),thisname,'.mat'],'md','hassolved')
                
                % Check if it has been solved:
                if ~hassolved
                    disp('......fetching results from server')
                    thesepaths=dir(['/work/ollie/mwolovic/ISSM_execution/',md.miscellaneous.name,'-*']);
                    % Sort by date:
                    folderdates=zeros(size(thesepaths));
                    for jj=1:length(thesepaths)
                        %folderdates(jj)=datenum(thesepaths(jj).date);
                        folderdates(jj)=thesepaths(jj).datenum;
                    end
                    [~,sortedind]=sort(folderdates);
                    % Load from server:
                    md=FlattenModelStructure(md,'expand');
                    md=solve(md,'sb','loadonly','true','restart',thesepaths(sortedind(end)).name);
                    md=FlattenModelStructure(md,'flatten');
                    % Save results:
                    disp('......saving post-solution matfile')
                    hassolved=1;
                    save([outputfile(1:slashnum(end)),thisname,'.mat'],'md','hassolved','MyInversionParameters','-v7.3')
                end
                
                % Record cost values:
                if isfield(md.results,'StressbalanceSolution')
                    DataCost(ii)=sum(md.results.StressbalanceSolution.J(end,AreReg==0));
                    RegularizationCost(ii)=sum(md.results.StressbalanceSolution.J(end,AreReg))/Lambda(ii);
                else
                    DataCost(ii)=NaN;
                    RegularizationCost(ii)=NaN;
                end
                
                % Communicate some more:
                disp(['Data cost= ',num2str(DataCost(ii)),', Regularization cost= ',num2str(RegularizationCost(ii))])
                disp('...')
                
            end
            
            % Save all results:
            disp('Saving overall results...')
            hassolved=1;
            save(outputfile,'Lambda','DataCost','RegularizationCost','BasicModel','MyInversionParameters','hassolved','-v7.3')
            
        end
        
    else
        
        % Check if the model has been solved:
        if ~hassolved
            disp('Loading results from server')
            thesepaths=dir(['/work/ollie/mwolovic/ISSM_execution/',md.miscellaneous.name,'-*']);
            % Sort by date:
            folderdates=zeros(size(thesepaths));
            for ii=1:length(thesepaths)
                %folderdates(ii)=datenum(thesepaths(ii).date);
                folderdates(ii)=thesepaths(ii).datenum;
            end
            [~,sortedind]=sort(folderdates);
            % Load from server:
            md=FlattenModelStructure(md,'expand');
            md=solve(md,'sb','loadonly','true','restart',thesepaths(sortedind(end)).name);
            md=FlattenModelStructure(md,'flatten');
            % Save results:
            disp('Saving post-solution matfile')
            hassolved=1;
            save(outputfile,'md','hassolved','MyInversionParameters','-v7.3')
        end
        
    end
    
end
%toc
%return

%% Display Results:

% Check whether to make any figures:
if strcmp(overallmode,'loadandanalyze')
    
    % L-curve figures:
    if dolcurve && exist('Lambda','var')
        
        % Check whether to do lcurve analysis and summary figure:
        if makelcurvefig
            
            %% Do L-Curve Analysis:
            
            % Identify outlier models:
            % prep:
            IsOutlier=isnan(DataCost);
            done=0;
            numiterations_outlierid=1;
            % loop:
            while done==0
                % record last outlier list:
                IsOutlier_last=IsOutlier;
                % isolateremaining good points:
                TheseDataCost=DataCost(~IsOutlier_last);
                TheseRegCost=RegularizationCost(~IsOutlier_last);
                % check remaining good points for new outliers:
                IsOutlier(~IsOutlier_last)=[TheseDataCost(1)>outlierthreshold*TheseDataCost(2);TheseDataCost(2:end-1)>outlierthreshold*sqrt(TheseDataCost(1:end-2).*TheseDataCost(3:end));...
                    TheseDataCost(end)>outlierthreshold*TheseDataCost(end-1)]|...
                    [TheseRegCost(1)>outlierthreshold*TheseRegCost(2);TheseRegCost(2:end-1)>outlierthreshold*sqrt(TheseRegCost(1:end-2).*TheseRegCost(3:end));...
                    TheseRegCost(end)>outlierthreshold*TheseRegCost(end-1)];
                % check for change in outlier list:
                if isequal(IsOutlier,IsOutlier_last)
                    % break from loop:
                    done=1;
                elseif numiterations_outlierid>=numlambda
                    error('Unable to identify outlier models')
                else
                    % advance counter:
                    numiterations_outlierid=numiterations_outlierid+1;
                end
            end
            
            % Exclude non-monotonic points from the beggining of the analysis:
            if excludenonmonojobs
                mindatacostind=find(DataCost==min(DataCost),1,'first');
                if mindatacostind>1
                    IsOutlier(1:mindatacostind-1)=1;
                end
            end
            if excludenonmonojreg
                maxregcostind=find(RegularizationCost==max(RegularizationCost),1,'first');
                if maxregcostind>1
                    IsOutlier(1:maxregcostind-1)=1;
                end
            end
            
            % Exclude an arbitrary number of points from the start of the
            % analysis:
            if num2excludefromstart>0
                IsOutlier(1:num2excludefromstart)=1;
            end
            
            % Convert points to log scale:
            LogDataCost_points=log(DataCost);
            LogRegCost_points=log(RegularizationCost);
            LogLambda_points=log(Lambda);
            
            % Compute lambda subsampling:
            LogLambda_path=linspace(LogLambda_points(find(~IsOutlier,1,'first')),LogLambda_points(find(~IsOutlier,1,'last')),numlambdasubsamples)';
            
            % Interpolate log(variables) onto the new lambdasubsampling:
            LogDataCost_path=interp1(LogLambda_points(~IsOutlier),LogDataCost_points(~IsOutlier),LogLambda_path,lambdainterpmethod);
            LogRegCost_path=interp1(LogLambda_points(~IsOutlier),LogRegCost_points(~IsOutlier),LogLambda_path,lambdainterpmethod);
            
            % Compute dlambda:
            dloglambda_path=LogLambda_path(2)-LogLambda_path(1);
            
            % Compute minimum lambda wavelength:
            minlambdawavelength=.5*(LogLambda_points(2)-LogLambda_points(1));
            
            % Generate wavelength sampling:
            LambdaWavelengths=linspace(minlambdawavelength,maxlambdawavelength,numlambdawavelengths);
            
            % Compute characteristic scales:
            if strcmp(curvenormstyle,'quad')
                % Estimate scales from quadratic best-fits:
                A=[ones(sum(~IsOutlier),1),LogLambda_points(~IsOutlier),LogLambda_points(~IsOutlier).^2];
                equation_data=(A'*A)\A'*LogDataCost_points(~IsOutlier);
                equation_reg=(A'*A)\A'*LogRegCost_points(~IsOutlier);
                Fit_data=A*equation_data;
                Fit_reg=A*equation_reg;
                Misfit_data=LogDataCost_points(~IsOutlier)-Fit_data;
                Misfit_reg=LogRegCost_points(~IsOutlier)-Fit_reg;
                curvescale=sqrt(equation_data(3)^2+equation_reg(3)^2);
                scatterscale=sqrt(mean(Misfit_data.^2+Misfit_reg.^2));
            elseif strcmp(curvenormstyle,'max')
                curvescale=1;
                scatterscale=1;
            else
                error('Parameter "curvenormstyle" must be a string equal to "quad" or "max".')
            end
            
            % Pre-allocate:
            LogDataCost_curves=zeros(numlambdasubsamples,numlambdawavelengths);
            LogRegCost_curves=zeros(numlambdasubsamples,numlambdawavelengths);
            CurvatureVariance_curves=zeros(1,numlambdawavelengths);
            ScatterVariance_curves=zeros(1,numlambdawavelengths);
            
            % Loop through potential smoothings:
            for ii=1:numlambdawavelengths
                % Check wavelength:
                if LambdaWavelengths(ii)>0
                    % Smooth the variables:
                    LogDataCost_curves(:,ii)=intuitive_lowpass(LogDataCost_path,LambdaWavelengths(ii)/dloglambda_path,'quasilocal');
                    LogRegCost_curves(:,ii)=intuitive_lowpass(LogRegCost_path,LambdaWavelengths(ii)/dloglambda_path,'quasilocal');
                else
                    % No smoothing required:
                    LogDataCost_curves(:,ii)=LogDataCost_path;
                    LogRegCost_curves(:,ii)=LogRegCost_path;
                end
                % Compute curvature of this curve:
                ThisLogDataCurvature=(LogDataCost_curves(1:end-2,ii)+LogDataCost_curves(3:end,ii)-2*LogDataCost_curves(2:end-1,ii))/(dloglambda_path^2);
                ThisLogRegCurvature=(LogRegCost_curves(1:end-2,ii)+LogRegCost_curves(3:end,ii)-2*LogRegCost_curves(2:end-1,ii))/(dloglambda_path^2);
                % Compute curvature variance:
                if strcmp(curvevarstyle,'var')
                    CurvatureVariance_curves(ii)=mean(ThisLogDataCurvature.^2+ThisLogRegCurvature.^2)/(curvescale^2);
                elseif strcmp(curvevarstyle,'rms')
                    CurvatureVariance_curves(ii)=sqrt(mean(ThisLogDataCurvature.^2+ThisLogRegCurvature.^2))/curvescale;
                else
                    error('Parameter "curvevarstyle" must be a string equal to "var" or "rms".')
                end
                % Compute scatter about this curve:
                ThisDataScatter=LogDataCost_points(~IsOutlier)-interp1(LogLambda_path,LogDataCost_curves(:,ii),LogLambda_points(~IsOutlier),lambdainterpmethod);
                ThisRegScatter=LogRegCost_points(~IsOutlier)-interp1(LogLambda_path,LogRegCost_curves(:,ii),LogLambda_points(~IsOutlier),lambdainterpmethod);
                % Compute scatter variance:
                if strcmp(curvevarstyle,'var')
                    ScatterVariance_curves(ii)=mean(ThisDataScatter.^2+ThisRegScatter.^2)/(scatterscale^2);
                else % rms
                    ScatterVariance_curves(ii)=sqrt(mean(ThisDataScatter.^2+ThisRegScatter.^2))/scatterscale;
                end
            end
            
            % Finish normalizing curvature variances:
            if strcmp(curvenormstyle,'max')
                CurvatureVariance_curves=CurvatureVariance_curves/max(CurvatureVariance_curves);
                ScatterVariance_curves=ScatterVariance_curves/max(ScatterVariance_curves);
            end
            
            % Determine best wavelength:
            TotalVariance_curves=ScatterVariance_curves+CurvatureVariance_curves;
            bestwavelengthind=find(TotalVariance_curves==min(TotalVariance_curves),1,'first');
            bestwavelength=LambdaWavelengths(bestwavelengthind);
            
            % Isolate best curve:
            LogDataCost_bestcurve=LogDataCost_curves(:,bestwavelengthind);
            LogRegCost_bestcurve=LogRegCost_curves(:,bestwavelengthind);
            
            % Compute gradient and curvature of best curve:
            LogDataGradient_bestcurve=(LogDataCost_bestcurve(2:end)-LogDataCost_bestcurve(1:end-1))/dloglambda_path;
            LogRegGradient_bestcurve=(LogRegCost_bestcurve(2:end)-LogRegCost_bestcurve(1:end-1))/dloglambda_path;
            LogDataCurvature_bestcurve=(LogDataCost_bestcurve(1:end-2)+LogDataCost_bestcurve(3:end)-2*LogDataCost_bestcurve(2:end-1))/(dloglambda_path^2);
            LogRegCurvature_bestcurve=(LogRegCost_bestcurve(1:end-2)+LogRegCost_bestcurve(3:end)-2*LogRegCost_bestcurve(2:end-1))/(dloglambda_path^2);
            TotalCurvature_bestcurve=LogDataCurvature_bestcurve+LogRegCurvature_bestcurve;
            
            % Indexing note:  curvature has size numlambdasubsamples-2.
            % ie, it is missing the first and last elements.  The "find"
            % commands below have been modified so that the resulting index
            % points to the correct position in the original list.
            
            % Identify best lambda:
            maxcurve=max(TotalCurvature_bestcurve(ceil(.5*bestwavelength/dloglambda_path):floor(numlambdasubsamples-.5*bestwavelength/dloglambda_path)));
            bestlambdaind=1+find(TotalCurvature_bestcurve==maxcurve,1,'first');
            bestlambda=exp(LogLambda_path(bestlambdaind));
            bestmodelind=find(abs(LogLambda_points-log(bestlambda))==min(abs(LogLambda_points-log(bestlambda))),1,'first');
            
            % Compute max/min allowable lambdas:
            bestlambdaind_minallowable=find(TotalCurvature_bestcurve(1:bestlambdaind)<.5*maxcurve,1,'last');
            if isempty(bestlambdaind_minallowable)
                bestlambdaind_minallowable=1;
            end
            bestlambdaind_maxallowable=bestlambdaind-1+find(TotalCurvature_bestcurve(bestlambdaind-1:end)<.5*maxcurve,1,'first');
            if isempty(bestlambdaind_maxallowable)
                bestlambdaind_maxallowable=numlambdasubsamples;
            end
            bestlambda_minallowable=exp(LogLambda_path(bestlambdaind_minallowable));
            bestlambda_maxallowable=exp(LogLambda_path(bestlambdaind_maxallowable));
            
            % Compute uncertainty:
            bestlambda_loguncertainty=.5*bestwavelength;
            
            % Display results:
            disp(['Min acceptable lambda=',num2str(bestlambda_minallowable)])
            disp(['Best lambda=',num2str(bestlambda)])
            disp(['Max acceptable lambda=',num2str(bestlambda_maxallowable)])
            disp(['Lambda uncertainty ratio=',num2str(exp(bestlambda_loguncertainty))])
            disp(['Drag coeff scale=',num2str(dragcoeffstdev)])
            disp(['Velocity scale=',num2str(velrms)])
            disp(['Best model Jobs=',num2str(DataCost(bestmodelind))])
            disp(['Best model Jreg=',num2str(RegularizationCost(bestmodelind))])
            
            % Make a structure with L-curve fitting parameters:
            LCurveFittingParameters=struct('numlambdasubsamples',numlambdasubsamples,'numlambdawavelengths',numlambdawavelengths,'maxlambdawavelength',maxlambdawavelength,...
                'outlierthreshold',outlierthreshold,'excludenonmonojobs',excludenonmonojobs,'excludenonmonojreg',excludenonmonojreg,'num2excludefromstart',num2excludefromstart,...
                'lambdainterpmethod',lambdainterpmethod,'curvenormstyle',curvenormstyle,'curvevarstyle',curvevarstyle);
            
            % Save best values:
            save(outputfile,'*best*','*_curves','LambdaWavelengths','IsOutlier','LCurveFittingParameters','-append')
            
            %toc
            %return
            
            %% Make L-Curve Figure:
            
            % Compute boxes for L-Curve figure:
            Boxes=cell(5,1);
            Boxes{1}=[horzaxesbuffer_lcurve,vertaxesbuffer_lcurve,lcurvefrac-horzaxesbuffer_lcurve-horznotextbuffer_lcurve,1-titlebuffer_lcurve-vertaxesbuffer_lcurve];
            availablehorzspace=1-lcurvefrac-2*(horzaxesbuffer_lcurve+horznotextbuffer_lcurve);
            availablevertspace=1-2*(vertaxesbuffer_lcurve+titlebuffer_lcurve);
            for ii=1:4
                d1=ceil(ii/2);
                d2=ii-(d1-1)*2;
                Boxes{ii+1}=[lcurvefrac+d2*horzaxesbuffer_lcurve+(d2-1)*(horznotextbuffer_lcurve+availablehorzspace/2),vertaxesbuffer_lcurve+(2-d1)*(vertaxesbuffer_lcurve+titlebuffer_lcurve+availablevertspace/2),availablehorzspace/2,availablevertspace/2];
            end
            
            % Make figure:
            figure
            
            % Plot cost vs cost, log scale:
            % Make subplot:
            subplot('Position',Boxes{1})
            % Plot model costs:
            hgoodpoints=plot(RegularizationCost(~IsOutlier),DataCost(~IsOutlier),'Marker','o','MarkerFaceColor','k','Color','k','LineStyle','none');
            hold on
            if sum(IsOutlier)>0
                houtliers=plot(RegularizationCost(IsOutlier),DataCost(IsOutlier),'Marker','o','MarkerFaceColor','m','Color','m','LineStyle','none');
            end
            % Plot best curve:
            hcurve=plot(exp(LogRegCost_bestcurve),exp(LogDataCost_bestcurve),'k','LineWidth',1);
            set(gca,'XScale','log')
            set(gca,'YScale','log')
            % Add lambda labels:
            for ii=1:4:numlambda
                thisind=find(LogLambda_path>=LogLambda_points(ii),1,'first');
                if thisind==numlambdasubsamples
                    thisind=thisind-1;
                end
                thisangle=atan2(LogDataCost_bestcurve(thisind+1)-LogDataCost_bestcurve(thisind),LogRegCost_bestcurve(thisind+1)-LogRegCost_bestcurve(thisind));
                plot(exp(log(RegularizationCost(ii))+[0,.2*cos(thisangle+pi/2)]),exp(log(DataCost(ii))+[0,.2*sin(thisangle+pi/2)]),'k','LineWidth',1)
                text(exp(log(RegularizationCost(ii))+.2*cos(thisangle+pi/2)),exp(log(DataCost(ii))+.2*sin(thisangle+pi/2)),['\lambda=',num2str(Lambda(ii))],...
                    'Color','k','HorizontalAlignment','right','VerticalAlignment','top','FontSize',12)
            end
            % Plot and label best value:
            angle1=atan2(LogDataCost_bestcurve(bestlambdaind+1)-LogDataCost_bestcurve(bestlambdaind),...
                LogRegCost_bestcurve(bestlambdaind+1)-LogRegCost_bestcurve(bestlambdaind));
            plot(exp(LogRegCost_bestcurve(bestlambdaind)),exp(LogDataCost_bestcurve(bestlambdaind)),'Color','r','MarkerFaceColor','r','Marker','d')
            lastpow10=10^floor(log10(bestlambda));
            plot(exp(LogRegCost_bestcurve(bestlambdaind)-[0,.2*cos(angle1+pi/2)]),exp(LogDataCost_bestcurve(bestlambdaind)-[0,.2*sin(angle1+pi/2)]),'r','LineWidth',1)
            text(exp(LogRegCost_bestcurve(bestlambdaind)+.2),exp(LogDataCost_bestcurve(bestlambdaind)+.2),{'Best value:';['\lambda=',num2str(.1*lastpow10*round(10*bestlambda/lastpow10))]},...
                'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k','FontSize',12)
            % Plot max/min allowable:
            hcorner=plot(exp(LogRegCost_bestcurve(bestlambdaind_minallowable:bestlambdaind_maxallowable)),...
                exp(LogDataCost_bestcurve(bestlambdaind_minallowable:bestlambdaind_maxallowable)),'Color','r','LineWidth',2);
            plot(exp(LogRegCost_bestcurve(bestlambdaind_minallowable)),exp(LogDataCost_bestcurve(bestlambdaind_minallowable)),'Color','r','MarkerFaceColor','r','Marker','d')
            plot(exp(LogRegCost_bestcurve(bestlambdaind_maxallowable)),exp(LogDataCost_bestcurve(bestlambdaind_maxallowable)),'Color','r','MarkerFaceColor','r','Marker','d')
            % Label max/min allowable:
            angle1=atan2(LogDataCost_bestcurve(bestlambdaind_minallowable+1)-LogDataCost_bestcurve(bestlambdaind_minallowable),...
                LogRegCost_bestcurve(bestlambdaind_minallowable+1)-LogRegCost_bestcurve(bestlambdaind_minallowable));
            angle2=atan2(LogDataCost_bestcurve(bestlambdaind_maxallowable)-LogDataCost_bestcurve(bestlambdaind_maxallowable-1),...
                LogRegCost_bestcurve(bestlambdaind_maxallowable)-LogRegCost_bestcurve(bestlambdaind_maxallowable-1));
            lastpow10_1=10^floor(log10(bestlambda_minallowable));
            lastpow10_2=10^floor(log10(bestlambda_maxallowable));
            plot(exp(LogRegCost_bestcurve(bestlambdaind_minallowable)-[0,.2*cos(angle1+pi/2)]),exp(LogDataCost_bestcurve(bestlambdaind_minallowable)-[0,.2*sin(angle1+pi/2)]),'r','LineWidth',1)
            text(exp(LogRegCost_bestcurve(bestlambdaind_minallowable)-.2*cos(angle1+pi/2)),exp(LogDataCost_bestcurve(bestlambdaind_minallowable)-.2*sin(angle1+pi/2)),...
                {'Min value:';['\lambda=',num2str(.1*lastpow10_1*round(10*bestlambda_minallowable/lastpow10_1))]},...
                'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k','FontSize',12)
            plot(exp(LogRegCost_bestcurve(bestlambdaind_maxallowable)-[0,.2*cos(angle2+pi/2)]),exp(LogDataCost_bestcurve(bestlambdaind_maxallowable)-[0,.2*sin(angle2+pi/2)]),'r','LineWidth',1)
            text(exp(LogRegCost_bestcurve(bestlambdaind_maxallowable)-.2*cos(angle2+pi/2)),exp(LogDataCost_bestcurve(bestlambdaind_maxallowable)-.2*sin(angle2+pi/2)),...
                {'Max value:';['\lambda=',num2str(.1*lastpow10_2*round(10*bestlambda_maxallowable/lastpow10_2))]},...
                'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k','FontSize',12)
            % Compute axes limits:
            thesexlims=10.^[floor(log10(min(RegularizationCost))),ceil(log10(max(RegularizationCost)))];
            theseylims=10.^[floor(log10(min([DataCost(~IsOutlier);exp(LogDataCost_bestcurve(bestlambdaind)-.5)]))),ceil(log10(max(DataCost(~IsOutlier))))];
            % Set up axes etc:
            xlim(thesexlims)
            ylim(theseylims)
            set(gca,'PlotBoxAspectRatio',[(log(thesexlims(2))-log(thesexlims(1)))/(log(theseylims(2))-log(theseylims(1))),1,1])
            xlabel('Regularization cost J_{reg} (unitless)','FontSize',13)
            ylabel('Data cost J_{obs} (unitless)','FontSize',13)
            set(gca,'FontSize',13)
            if sum(IsOutlier)>0
                legend([hgoodpoints,houtliers,hcurve,hcorner],'Individual model results (good points)','Individual model results (outliers)','Smoothed tradeoff curve','Corner region','location','NorthEast','FontSize',13)
            else
                legend([hgoodpoints,hcurve,hcorner],'Individual model results','Smoothed tradeoff curve','Corner region','location','NorthEast','FontSize',13)
            end
            title('a) L-Curve Plot (cost vs cost)','FontSize',14)
            
            % Plot both components against lambda values:
            % Make subplot:
            subplot('Position',Boxes{2})
            % Plot good model points:
            h1=plot(Lambda(~IsOutlier),DataCost(~IsOutlier),'Marker','o','MarkerFaceColor','b','Color','b','LineStyle','none');
            hold on
            h2=plot(Lambda(~IsOutlier),RegularizationCost(~IsOutlier),'Marker','o','MarkerFaceColor','r','Color','r','LineStyle','none');
            % Plot bad model points:
            if sum(IsOutlier)>0
                plot(Lambda(IsOutlier),DataCost(IsOutlier),'Marker','o','MarkerFaceColor','none','Color','b','LineStyle','none')
                plot(Lambda(IsOutlier),RegularizationCost(IsOutlier),'Marker','o','MarkerFaceColor','none','Color','r','LineStyle','none')
            end
            % Plot smooth curve:
            h3=plot(exp(LogLambda_path),exp(LogDataCost_bestcurve),'b','LineWidth',1);
            h4=plot(exp(LogLambda_path),exp(LogRegCost_bestcurve),'r','LineWidth',1);
            % Compute axes limits:
            thesexlims=10.^[floor(log10(lambdamin)),ceil(log10(lambdamax))];
            theseylims=10.^[floor(log10(min([DataCost(~IsOutlier);RegularizationCost(~IsOutlier)]))),ceil(log10(max([DataCost(~IsOutlier);RegularizationCost(~IsOutlier)])))];
            % Set up axes etc:
            set(gca,'XScale','log')
            set(gca,'YScale','log')
            xlim(thesexlims)
            ylim(theseylims)
            set(gca,'XTick',10.^(log10(thesexlims(1)):log10(thesexlims(2))))
            xlabel('\lambda value (unitless)','FontSize',13)
            ylabel('Cost Components J (unitless)','FontSize',13)
            set(gca,'FontSize',13)
            legend([h1;h2;h3;h4],'J_{obs} (models)','J_{reg} (models)','J_{obs} (curve)','J_{reg} (curve)','location','SouthWest','FontSize',13)
            %             if sum(IsOutlier)>0
            %                 legend('J_{obs} (individual models)','J_{reg} (individual models)','J_{obs} (outlier models)','J_{reg} (outlier models)','J_{obs} (smooth curve)','J_{reg} (smooth curve)','location','NorthWest')
            %             else
            %                 legend('J_{obs} (individual models)','J_{reg} (individual models)','J_{obs} (smooth curve)','J_{reg} (smooth curve)','location','NorthWest')
            %             end
            title('b) Cost components','FontSize',14)
            
            % Plot gradients against lambda values:
            % Make subplot:
            subplot('Position',Boxes{3})
            % Plot smooth curves:
            plot(exp(.5*(LogLambda_path(1:end-1)+LogLambda_path(2:end))),LogDataGradient_bestcurve,'b','LineWidth',1)
            hold on
            plot(exp(.5*(LogLambda_path(1:end-1)+LogLambda_path(2:end))),LogRegGradient_bestcurve,'r','LineWidth',1)
            % Plot zero line:
            plot(thesexlims,[0,0],'--k')
            % Compute axes limits:
            thesexlims=10.^[floor(log10(lambdamin)),ceil(log10(lambdamax))];
            theseylims=[min([LogDataGradient_bestcurve;LogRegGradient_bestcurve]),max([LogDataGradient_bestcurve;LogRegGradient_bestcurve])];
            theseylims=theseylims+[-1,1]*.125*(theseylims(2)-theseylims(1));
            % Set up axes etc:
            set(gca,'XScale','log')
            xlim(thesexlims)
            ylim(theseylims)
            set(gca,'XTick',10.^(log10(thesexlims(1)):log10(thesexlims(2))))
            xlabel('\lambda value (unitless)','FontSize',13)
            ylabel('Cost gradient, $\frac{d(ln(J))}{d(ln(\lambda))}$','interpreter','latex','FontSize',13)
            set(gca,'FontSize',13)
            legend('J_{obs}','J_{reg}','location','NorthWest','FontSize',13)
            title('c) Cost gradients','FontSize',14)
            
            % Plot cost curvature against lambda values:
            % Make subplot:
            subplot('Position',Boxes{4})
            % Plot best curve curvatures:
            plot(exp(LogLambda_path(2:end-1)),LogDataCurvature_bestcurve,'b','LineWidth',1)
            hold on
            plot(exp(LogLambda_path(2:end-1)),LogRegCurvature_bestcurve,'r','LineWidth',1)
            plot(exp(LogLambda_path(2:end-1)),TotalCurvature_bestcurve,'k','LineWidth',1)
            % Show best value and uncertainty:
            plot(bestlambda,maxcurve,'Color','k','MarkerFaceColor','k','Marker','o')
            text(bestlambda,maxcurve*1.025,{'Best';'value'},'HorizontalAlignment','center','VerticalAlignment','bottom','FontSize',12)
            plot(bestlambda*[1/exp(bestlambda_loguncertainty),exp(bestlambda_loguncertainty)],maxcurve*[1,1],'k')
            plot((bestlambda/exp(bestlambda_loguncertainty))*[1,1],maxcurve*[.975,1.025],'k')
            plot(bestlambda*exp(bestlambda_loguncertainty)*[1,1],maxcurve*[.975,1.025],'k')
            % Show max/min lambda and uncertainties:
            %plot(bestlambda*[1,1],maxcurve*[.5,1],'--k')
            %plot([bestlambda_minallowable,bestlambda_maxallowable],.5*maxcurve*[1,1],'--k')
            plot(bestlambda_minallowable,.5*maxcurve,'Color','k','MarkerFaceColor','k','Marker','o')
            plot(bestlambda_minallowable*[1/exp(bestlambda_loguncertainty),exp(bestlambda_loguncertainty)],.5*maxcurve*[1,1],'k')
            plot((bestlambda_minallowable/exp(bestlambda_loguncertainty))*[1,1],.5*maxcurve*[.975,1.025],'k')
            plot(bestlambda_minallowable*exp(bestlambda_loguncertainty)*[1,1],.5*maxcurve*[.975,1.025],'k')
            plot(bestlambda_minallowable*[1,.75],.5*maxcurve+[0,.01],'k')
            text(bestlambda_minallowable*.75,.5*maxcurve+.01,{'Min';'value'},'HorizontalAlignment','right','VerticalAlignment','bottom','Color','k','FontSize',12)
            plot(bestlambda_maxallowable,.5*maxcurve,'Color','k','MarkerFaceColor','k','Marker','o')
            plot(bestlambda_maxallowable*[1/exp(bestlambda_loguncertainty),exp(bestlambda_loguncertainty)],.5*maxcurve*[1,1],'k')
            plot((bestlambda_maxallowable/exp(bestlambda_loguncertainty))*[1,1],.5*maxcurve*[.975,1.025],'k')
            plot(bestlambda_maxallowable*exp(bestlambda_loguncertainty)*[1,1],.5*maxcurve*[.975,1.025],'k')
             plot(bestlambda_maxallowable*[1,2],.5*maxcurve+[0,.01],'k')
            text(bestlambda_maxallowable*2,.5*maxcurve+.01,{'Max';'value'},'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k','FontSize',12)
            % Determine axes limits:
            thesexlims=10.^[floor(log10(lambdamin)),ceil(log10(lambdamax))];
            theseylims=TotalCurvature_bestcurve(bestlambdaind-1)*[-.25,1.25];
            % Plot zero line:
            plot(thesexlims,[0,0],'--k')
            % Set up axes etc:
            xlim(thesexlims)
            ylim(theseylims)
            set(gca,'XScale','log')
            set(gca,'XTick',10.^(log10(thesexlims(1)):log10(thesexlims(2))))
            xlabel('\lambda value (unitless)','FontSize',13)
            ylabel('Cost curvature, $\frac{d^{2}(ln(J))}{d(ln(\lambda))^{2}}$','interpreter','latex','FontSize',13)
            set(gca,'FontSize',13)
            title({'d) Cost curvature';'(selection of optimal \lambda values)'},'FontSize',14)
            if log(bestlambda)>mean(log([lambdamin,lambdamax]))
                locationstring='NorthWest';
            else
                locationstring='NorthEast';
            end
            legend('J_{obs}','J_{reg}','Total','location',locationstring,'FontSize',13)
            
            % Plot variance against wavelength:
            % Make subplot:
            subplot('Position',Boxes{5})
            % Plot RMS values:
            plot(LambdaWavelengths,ScatterVariance_curves,'b','LineWidth',1)
            hold on
            plot(LambdaWavelengths,CurvatureVariance_curves,'r','LineWidth',1)
            plot(LambdaWavelengths,TotalVariance_curves,'k','LineWidth',1)
            % Plot and label best value:
            plot(LambdaWavelengths(bestwavelengthind),TotalVariance_curves(bestwavelengthind),'Color','k','MarkerFaceColor','k','Marker','o')
            lastpow10=min(1,10^floor(log10(bestwavelength)));
            text(LambdaWavelengths(bestwavelengthind)+.5,TotalVariance_curves(bestwavelengthind)+.35,{['Best wavelength=',num2str(.01*lastpow10*round(100*bestwavelength/lastpow10))];...
                ['Equiv. uncertainty ratio=',num2str(.01*round(100*exp(bestlambda_loguncertainty)))]},'HorizontalAlignment','center','VerticalAlignment','bottom','Color','k','FontSize',12)
            plot(LambdaWavelengths(bestwavelengthind)+[0,0.5],TotalVariance_curves(bestwavelengthind)+[0,.35],'k')
            % Determine axes limits:
            thesexlims=[0,maxlambdawavelength];
            theseylims=max(TotalVariance_curves)*[0,1.05];
            % Set up axes etc:
            xlim(thesexlims)
            ylim(theseylims)
            xlabel('Wavelength (distance in ln(\lambda) space)','FontSize',13)
            ylabel('Variance (unitless)','FontSize',13)
            set(gca,'FontSize',13)
            title({'e) Variance of tradeoff curve';'(selection of optimal smoothing)'},'FontSize',14)
            legend('Scatter variance','Curvature variance','Total variance','location','NorthEast','FontSize',13)
            
            
            % Save figure:
            thisfigname=[outputfile(1:end-4),'_LCurvePlots6a.png'];
            set(gcf,'PaperSize',lcurvepagesize)
            set(gcf,'PaperPosition',[0,0,lcurvepagesize])
            print('-dpng',thisfigname,['-r',num2str(resolution)])
            
            % Close figure:
            close
            
            %toc
            %return
            
            
        end
        
        %% Make results comparison figure:
        
        % Check whether to make L-curve comparison figs:
        if makelcurvemapcomparisonfigs
            
            % Make figure:
            figure(1)
            
            % Compute number of subplots:
            numtotalplots=numlambda;
            numvertplots=ceil(numtotalplots/4);
            numhorzplots=4;
            % Compute subplot dimensions:
            % Compute available space:
            availablehorzspace=1-numhorzplots*(horznotextbuffer_map+horzaxesbuffer_map+horzaxeslabelbuffer_map);
            availablevertspace=1-numvertplots*(titlebuffer_map+vertaxesbuffer_map+vertaxeslabelbuffer_map);
            % Compute boxes:
            Boxes=cell(numvertplots,numhorzplots);
            for ii=1:numvertplots
                for jj=1:numhorzplots
                    Boxes{numvertplots-ii+1,jj}=[jj*(horzaxesbuffer_map+horzaxeslabelbuffer_map)+(jj-1)*(horznotextbuffer_map+availablehorzspace/numhorzplots),...
                        ii*(vertaxesbuffer_map+vertaxeslabelbuffer_map)+(ii-1)*(titlebuffer_map+availablevertspace/numvertplots),...
                        availablehorzspace/numhorzplots,availablevertspace/numvertplots];
                end
            end
            % Make boxes transpose:
            Boxes=Boxes';
            
            % Make the faded jet colormap:
            FadedJet=colormap('jet');
            FadedJet(3*round(size(FadedJet,1)/8)+1:round(size(FadedJet,1)/2),1)=linspace(0,1,round(size(FadedJet,1)/8))';
            FadedJet(3*round(size(FadedJet,1)/8)+1:round(size(FadedJet,1)/2),3)=1;
            FadedJet(round(size(FadedJet,1)/2)+1:5*round(size(FadedJet,1)/8),1)=1;
            FadedJet(round(size(FadedJet,1)/2)+1:5*round(size(FadedJet,1)/8),3)=linspace(1,0,round(size(FadedJet,1)/8))';
            % Make the clipped jet colormap:
            ClippedJet=colormap('jet');
            ClippedJet=ClippedJet(round(size(ClippedJet,1)/8):end-round(size(ClippedJet,1)/8),:);
            % Get regular jet:
            RegularJet=colormap('jet');
            
            % Check whether to make drag comparisons:
            if includedrag
                
                % Identify elements to plot drag for:
                TheseElements=IceAreaFraction_elements>=.5&OceanAreaFraction_elements<=.5;
                
                % Drag Coefficient Figure:
                % Loop through lambda values:
                for ii=1:numlambda
                    % Construct this file name:
                    slashnum=strfind(outputfile,'/');
                    thisbase=10^floor(log10(Lambda(ii)));
                    thisname=[outputfile(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))];
                    % Load this model:
                    load([outputfile(1:slashnum(end)),thisname,'.mat'],'md')
                    % Plot drag coefficient:
                    % Make subplot:
                    subplot('Position',Boxes{ii})
                    % Plot patch:
                    patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d)))
                    % Set up shit:
                    caxis(log10([mindragcoeff,maxdragcoeff]))
                    xlim(xlims/1000)
                    ylim(ylims/1000)
                    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
                    xlabel('Easting (km)')
                    ylabel('Northing (km)')
                    hc=colorbar;
                    if ~strcmp(nsource,'none')
                        set(get(hc,'Ylabel'),'String',['log_{10}(Pa^{1/2}(m/s)^{-1/',num2str(2*m),'})'])
                    else
                        if r==1
                            set(get(hc,'Ylabel'),'String',['log_{10}((m/s)^{-1/',num2str(2*m),'})'])
                        else
                            set(get(hc,'Ylabel'),'String',['log_{10}(Pa^{',num2str((1-r)/2),'}(m/s)^{-1/',num2str(2*m),'})'])
                        end
                    end
                    title([char(96+ii),') Drag Coefficient for Lambda=',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))])
                end
                % Save figure:
                thispagesize=[numhorzplots*widthpermap,numvertplots*heightpermap];
                thisfigname=[outputfile(1:end-4),'_DragCoeffComparison.png'];
                set(gcf,'PaperSize',thispagesize)
                set(gcf,'PaperPosition',[0,0,thispagesize])
                print('-dpng',thisfigname,['-r',num2str(resolution)])
                
                % Drag Figure:
                % Loop through lambda values:
                for ii=1:numlambda
                    % Construct this file name:
                    slashnum=strfind(outputfile,'/');
                    thisbase=10^floor(log10(Lambda(ii)));
                    thisname=[outputfile(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))];
                    % Load this model:
                    load([outputfile(1:slashnum(end)),thisname,'.mat'],'md')
                    % Compute drag:
                    if exist('CoeffConversionParameters','var')
                        ThisDragCoeff_vertices=md.friction.coefficient(1:numvertices2d);
                    else
                        ThisDragCoeff_vertices=md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d);
                    end
                    ThisDrag=(ThisDragCoeff_vertices.^2).*((md.results.StressbalanceSolution.Vel(1:numvertices2d)/secondsperyear).^(1/m));
                    if ~strcmp(nsource,'none')
                        ThisDrag=ThisDrag.*(md.friction.effective_pressure.^r);
                    end
                    % Set drag to zero under floating shelf and in ice-free areas:
                    ThisDrag(OceanLevelSet_vertices<0|IceLevelSet_vertices>0)=0;
                    % Plot Drag:
                    % Make subplot:
                    subplot('Position',Boxes{ii})
                    % Make plot:
                    patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',ThisDrag/1000)
                    % Set up shit:
                    caxis(draglims)
                    xlim(xlims/1000)
                    ylim(ylims/1000)
                    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
                    xlabel('Easting (km)')
                    ylabel('Northing (km)')
                    hc=colorbar;
                    set(get(hc,'Ylabel'),'String','kPa')
                    title([char(96+ii),') Basal Drag for Lambda=',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))])
                end
                % Save figure:
                thispagesize=[numhorzplots*widthpermap,numvertplots*heightpermap];
                thisfigname=[outputfile(1:end-4),'_DragComparison.png'];
                set(gcf,'PaperSize',thispagesize)
                set(gcf,'PaperPosition',[0,0,thispagesize])
                print('-dpng',thisfigname,['-r',num2str(resolution)])
                
            end
            
            % Check whether to make rheology comparisons:
            if includerheology
                
                % Determine which elements to plot rheology for:
                ActiveMask_bbar_elements=mean(ActiveMask_bbar(elements2d),2);
                TheseElements=ActiveMask_bbar_elements>=.5;
                
                % Bbar Figure:
                % Loop through lambda values:
                for ii=1:numlambda
                    % Construct this file name:
                    slashnum=strfind(outputfile,'/');
                    thisbase=10^floor(log10(Lambda(ii)));
                    thisname=[outputfile(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))];
                    % Load this model:
                    load([outputfile(1:slashnum(end)),thisname,'.mat'],'md')
                    % Get this Bbar field:
                    if numlayers>1 && ~strcmp(flowequation,'SSA')
                        ThisBbar=DepthAverage(md,md.results.StressbalanceSolution.MaterialsRheologyBbar);
                    else
                        ThisBbar=md.results.StressbalanceSolution.MaterialsRheologyBbar;
                    end
                    % Plot Bbar coefficient:
                    % Make subplot:
                    subplot('Position',Boxes{ii})
                    % Make Plot:
                    patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(ThisBbar))
                    % Set up shit:
                    caxis(log10([minbbar,maxbbar]))
                    xlim(xlims/1000)
                    ylim(ylims/1000)
                    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
                    xlabel('Easting (km)')
                    ylabel('Northing (km)')
                    hc=colorbar;
                    %set(get(hc,'Ylabel'),'String',['log_{10}(Pa s^{1/',num2str(n),'})'])
                    set(get(hc,'Ylabel'),'String',['log_{10}(Pa s^{1/',num2str(n),'})'])
                    title([char(96+ii),') Vertically Averaged Rheology (Bbar) for Lambda=',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))])
                end
                % Save figure:
                thispagesize=[numhorzplots*widthpermap,numvertplots*heightpermap];
                thisfigname=[outputfile(1:end-4),'_BbarComparison.png'];
                set(gcf,'PaperSize',thispagesize)
                set(gcf,'PaperPosition',[0,0,thispagesize])
                print('-dpng',thisfigname,['-r',num2str(resolution)])
                
                % Viscosity Figure:
                % Loop through lambda values:
                for ii=1:numlambda
                    % Construct this file name:
                    slashnum=strfind(outputfile,'/');
                    thisbase=10^floor(log10(Lambda(ii)));
                    thisname=[outputfile(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))];
                    % Load this model:
                    load([outputfile(1:slashnum(end)),thisname,'.mat'],'md')
                    % Compute viscosity:
                    md.materials.rheology_B=md.results.StressbalanceSolution.MaterialsRheologyBbar;
                    ThisViscosity_elements=ComputeViscosity(md);
                    % Compute vertically averaged viscosity:
                    if md.mesh.numberoflayers>1
                        DZhat=(md.mesh.z(1+numvertices2d*(1:numlayers-1)')-md.mesh.z(1+numvertices2d*(0:numlayers-2)'))/md.geometry.thickness(1);
                        ThisViscosityBar_elements2d=sum(reshape(ThisViscosity_elements,[numelements2d,numlayers-1]).*repmat(DZhat',[numelements2d,1]),2);
                    else
                        ThisViscosityBar_elements2d=ThisViscosity_elements;
                    end
                    % Plot viscosity:
                    % Make subplot:
                    subplot('Position',Boxes{ii})
                    % Make plot:
                    patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d,'EdgeColor','none','FaceColor','flat','FaceVertexCData',log10(ThisViscosityBar_elements2d))
                    % Set up shit:
                    caxis(viscositylims)
                    xlim(xlims/1000)
                    ylim(ylims/1000)
                    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
                    xlabel('Easting (km)')
                    ylabel('Northing (km)')
                    hc=colorbar;
                    set(get(hc,'Ylabel'),'String','log_{10}(Pa*s)')
                    title([char(96+ii),') Vertically Averaged Viscosity for Lambda=',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))])
                end
                % Save figure:
                thispagesize=[numhorzplots*widthpermap,numvertplots*heightpermap];
                thisfigname=[outputfile(1:end-4),'_ViscosityComparison.png'];
                set(gcf,'PaperSize',thispagesize)
                set(gcf,'PaperPosition',[0,0,thispagesize])
                print('-dpng',thisfigname,['-r',num2str(resolution)])
                
            end
            
            % Close figure:
            close
            
        end
        
        %% Make FFT comparison figure:
        
        % Check whether to make L-curve comparison figs:
        if makelcurvefftcomparisonfig
            
            % Hardcoded things:
            lambdaskip=4;
            dx_grid=500;
            
            % Generate grid:
            X_grid=dx_grid*[floor(xlims(1)/dx_grid):1:ceil(xlims(2)/dx_grid)];
            Y_grid=dx_grid*[floor(ylims(1)/dx_grid):1:ceil(ylims(2)/dx_grid)]';
            ysize=length(Y_grid);
            xsize=length(X_grid);
            
            % Generate positive wavenumber grid:
            ks=1/dx_grid;
            Kx=ks*(0:(xsize/2))/xsize;
            Ky=ks*(0:(ysize/2))/ysize;
            
            % Generate full wavenumber grid:
            if rem(xsize,2)==1
                Kx=[Kx,fliplr(-Kx(2:end))];
            else
                Kx=[Kx,fliplr(-Kx(2:end-1))];
            end
            if rem(ysize,2)==1
                Ky=[Ky,fliplr(-Ky(2:end))]';
            else
                Ky=[Ky,fliplr(-Ky(2:end-1))]';
            end
            [Kx_mesh,Ky_mesh]=meshgrid(Kx,Ky);
            Kmag=sqrt(Kx_mesh.^2+Ky_mesh.^2);
            
            % Generate 1D spectral sampling:
            Kmag_1d=ks*(0:(min(xsize,ysize)/2))/min(xsize,ysize);
            num1dk=length(Kmag_1d);
            dk=Kmag_1d(2);
            
            % Pre-allocate:
            Power_mask_1d=zeros(num1dk,1);
            Power_result_1d=zeros(num1dk,numlambda);
            
            % Interpolate active mask to the grid:
            if includedrag
                ActiveMask_grid=InterpFromMeshToGrid(elements2d,x2d,y2d,double(ActiveMask_drag),X_grid,Y_grid,0);
                ActiveMask_elements=mean(ActiveMask_drag(elements2d),2);
            else
                ActiveMask_grid=InterpFromMeshToGrid(elements2d,x2d,y2d,double(ActiveMask_bbar),X_grid,Y_grid,0);
                ActiveMask_elements=mean(ActiveMask_bbar(elements2d),2);
            end
            
            % Loop through lambda values:
            maskmultiplier=Inf;
            for ii=1:lambdaskip:numlambda
                % Communicate:
                disp(['ii=',num2str(ii),'/',num2str(numlambda)])
                % Construct this file name:
                slashnum=strfind(outputfile,'/');
                thisbase=10^floor(log10(Lambda(ii)));
                thisname=[outputfile(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))];
                % Load this model:
                load([outputfile(1:slashnum(end)),thisname,'.mat'],'md')
                % Check what we're inverting:  (ASSUMES ONLY INVERTING ONE)
                if includedrag
                    % Get relevant drag coefficient:
                    if exist('CoeffConversionParameters','var')
                        ThisDragCoeff_vertices=md.friction.coefficient(1:numvertices2d);
                    else
                        ThisDragCoeff_vertices=md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d);
                    end
                    % Compute this mean:
                    %meanresult=sum(Area.*ActiveMask_elements.*mean(md.results.StressbalanceSolution.FrictionCoefficient(elements2d),2))/sum(Area.*ActiveMask_elements);
                    meanresult=sum(Area.*ActiveMask_elements.*log(max(1,mean(ThisDragCoeff_vertices(elements2d),2)/mindragcoeff)))/sum(Area.*ActiveMask_elements);
                    % Interpolate to regular grid:
                    %Result_grid=InterpFromMeshToGrid(elements2d,x2d,y2d,md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d).*FastFlowWeight_vertices,X_grid,Y_grid,0);
                    %Result_grid=InterpFromMeshToGrid(elements2d,x2d,y2d,md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d).*ActiveMask_drag,X_grid,Y_grid,0);
                    Result_grid=InterpFromMeshToGrid(elements2d,x2d,y2d,log(max(1,ThisDragCoeff_vertices/mindragcoeff)).*ActiveMask_drag,X_grid,Y_grid,0);
                else
                    % Compute this mean:
                    meanresult=sum(Area.*ActiveMask_elements.*log(max(1,mean(md.results.StressbalanceSolution.MaterialsRheologyBbar(elements2d),2)/minbbar)))/sum(Area.*ActiveMask_elements);
                    % Interpolate to regular grid:
                    Result_grid=InterpFromMeshToGrid(elements2d,x2d,y2d,log(max(1,md.results.StressbalanceSolution.MaterialsRheologyBbar(1:numvertices2d)/minbbar)).*ActiveMask_bbar,X_grid,Y_grid,0);
                end
                % Update overall min mean result:
                maskmultiplier=min(maskmultiplier,meanresult);
                % Ensure exactly zero outside of the active mask:
                Result_grid(ActiveMask_grid<=.01)=0;
                % Take the fft:
                Result_fft=fft2(Result_grid);
                % Loop through 1D spectral sampling:
                for jj=1:length(Kmag_1d)
                    TheseK=abs(Kmag-Kmag_1d(jj))<4*dk;
                    TheseWeights=exp(-.5*(((Kmag(TheseK)-Kmag_1d(jj))/dk).^2));
                    Power_result_1d(jj,ii)=sum(TheseWeights.*(abs(Result_fft(TheseK)).^2))/sum(TheseWeights);
                end
            end
            
            % Take the spectrum of the active mask:
            % Take the fft:
            ActiveMask_fft=fft2(ActiveMask_grid*maskmultiplier);
            % Loop through 1D spectral sampling:
            for jj=1:length(Kmag_1d)
                TheseK=abs(Kmag-Kmag_1d(jj))<4*dk;
                TheseWeights=exp(-.5*(((Kmag(TheseK)-Kmag_1d(jj))/dk).^2));
                Power_mask_1d(jj)=sum(TheseWeights.*(abs(ActiveMask_fft(TheseK)).^2))/sum(TheseWeights);
            end
            
            % Make a figure:
            figure
            % Get regular jet:
            RegularJet=colormap('jet');
            % Call first subplot:
            subplot(2,1,1)
            % Plot mask spectrum:
            plot(Kmag_1d*1000,Power_mask_1d,'k')
            hold on
            % Loop through lambda values:
            for ii=1:lambdaskip:numlambda
                % Get this color:
                thiscolor=interp1(linspace(0,1,size(RegularJet,1))',RegularJet,max(0,min(1,(log(Lambda(ii))-log(lambdamin))/(log(lambdamax)-log(lambdamin)))));
                % Plot:
                plot(Kmag_1d*1000,Power_result_1d(:,ii),'Color',thiscolor)
            end
            % Set up shit:
            caxis(log10([lambdamin,lambdamax]))
            set(gca,'XScale','log')
            set(gca,'YScale','log')
            xlabel('Wavenumber (1/km)')
            ylabel('Power')
            hc=colorbar;
            set(get(hc,'Ylabel'),'String','log_{10}(Lambda)')
            thesexlims=get(gca,'Xlim');
            theseylims=get(gca,'Ylim');
            text(thesexlims(1),theseylims(1),'Black=Spectrum of Active Mask','HorizontalAlignment','left','VerticalAlignment','bottom','Color','k')
            title('a) Spectral Analysis of Inversion Results')
            % Call second subplot:
            subplot(2,1,2)
            % Loop through lambda values:
            for ii=1:lambdaskip:numlambda
                % Get this color:
                thiscolor=interp1(linspace(0,1,size(RegularJet,1))',RegularJet,max(0,min(1,(log(Lambda(ii))-log(lambdamin))/(log(lambdamax)-log(lambdamin)))));
                % Plot:
                plot(Kmag_1d*1000,Power_result_1d(:,ii)./Power_mask_1d,'Color',thiscolor)
                hold on
            end
            % Set up shit:
            caxis(log10([lambdamin,lambdamax]))
            set(gca,'XScale','log')
            set(gca,'YScale','log')
            xlabel('Wavenumber (1/km)')
            ylabel('Power/Power(mask)')
            hc=colorbar;
            set(get(hc,'Ylabel'),'String','log_{10}(Lambda)')
            title('b) Spectra Normalized by Mask Spectrum')
            % Save figure:
            thispagesize=[10,15];
            thisfigname=[outputfile(1:end-4),'_SpectralAnalysis.png'];
            set(gcf,'PaperSize',thispagesize)
            set(gcf,'PaperPosition',[0,0,thispagesize])
            print('-dpng',thisfigname,['-r',num2str(resolution)])
            % Close figure:
            close
            
        end
        
        %% Make Flux Divergence Comparison Figure:
        
        if makelcurvedivergencecomparisonfig
            
            % Get active mask on the elements:
            if includedrag && includerheology
                ActiveMask_elements=mean(ActiveMask_bbar(elements2d)|ActiveMask_drag(elements2d),2);
            elseif includedrag
                ActiveMask_elements=mean(ActiveMask_drag(elements2d),2);
            else
                ActiveMask_elements=mean(ActiveMask_bbar(elements2d),2);
            end
            
            % Pre-allocate:
            FluxDivergenceRMS=zeros(numlambda,1);
            
            % Loop through lambda values:
            for ii=1:numlambda
                % Construct this file name:
                slashnum=strfind(outputfile,'/');
                thisbase=10^floor(log10(Lambda(ii)));
                thisname=[outputfile(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(ii)/thisbase))];
                % Load this model:
                load([outputfile(1:slashnum(end)),thisname,'.mat'],'md')
                % Compute flux divergence:
                % Compute depth-average velocity:
                if numlayers>1 && ~strcmp(flowequation,'SSA')
                    md.geometry.thickness=repmat(md.geometry.thickness,[numlayers,1]);
                    Ubar_vertices=DepthAverage(md,md.results.StressbalanceSolution.Vx);
                    Vbar_vertices=DepthAverage(md,md.results.StressbalanceSolution.Vy);
                    md.geometry.thickness=md.geometry.thickness(1:numvertices2d);
                else
                    Ubar_vertices=md.results.StressbalanceSolution.Vx;
                    Vbar_vertices=md.results.StressbalanceSolution.Vy;
                end
                % Compute ice flux:
                FluxX_vertices=Ubar_vertices.*md.geometry.thickness(1:numvertices2d);
                FluxY_vertices=Vbar_vertices.*md.geometry.thickness(1:numvertices2d);
                % Compute flux divergence:
                FluxDivergence_elements=sum(BasisGradX_elements.*FluxX_vertices(elements2d),2)+sum(BasisGradY_elements.*FluxY_vertices(elements2d),2);
                % Compute RMS:
                FluxDivergenceRMS(ii)=sqrt(sum((FluxDivergence_elements.^2).*Area.*ActiveMask_elements)/sum(Area.*ActiveMask_elements));
            end
            
            % Make a figure:
            figure(1)
            % Plot:
            plot(Lambda,FluxDivergenceRMS,'k')
            % Set up shit:
            set(gca,'XScale','log')
            xlabel('Lambda Value (unitless)')
            ylabel('RMS of Flux Divergence (m/yr)')
            title('RMS variability of Flux Divergence in Inversion Active Area')
            % Save figure:
            thispagesize=[10,8];
            thisfigname=[outputfile(1:end-4),'_FluxDivergenceRMS.png'];
            set(gcf,'PaperSize',thispagesize)
            set(gcf,'PaperPosition',[0,0,thispagesize])
            print('-dpng',thisfigname,['-r',num2str(resolution)])
            % Close figure:
            close
            
        end
        
        
        %% Load Best Model:
        
        % Check if we need to load pre-conversion file:
        if exist('CoeffConversionParameters','var')
            % Load old L-curve data:
            load(outputfile,'OldLCurveData')
            % Get model ind:
            thismodelind=find(DataCost==min(DataCost),1,'first');
            %             if isempty(CoeffConversionParameters.newflowequation)
            %                 thismodelind=OldLCurveData.bestmodelind;
            %             else
            %                 thismodelind=find(abs(log(Lambda)-log(OldLCurveData.bestlambda_minallowable))==min(abs(log(Lambda)-log(OldLCurveData.bestlambda_minallowable))),1,'first');
            %             end
            % Construct file name:
            slashnum=strfind(CoeffConversionParameters.inputfile,'/');
            thisbase=10^floor(log10(Lambda(thismodelind)));
            thisname=[CoeffConversionParameters.inputfile(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(thismodelind)/thisbase))];
            % Load original model:
            OriginalModel=load([CoeffConversionParameters.inputfile(1:slashnum(end)),thisname,'.mat'],'md');
        else
            thismodelind=bestmodelind;
        end
        
        % Construct file name for best value:
        slashnum=strfind(outputfile,'/');
        thisbase=10^floor(log10(Lambda(thismodelind)));
        thisname=[outputfile(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(thismodelind)/thisbase))];
        
        % Load model for best model:
        load([outputfile(1:slashnum(end)),thisname,'.mat'],'md')
        
        
    end
    %toc
    %return
    
    %% Compute misfit components:
    
    % Isolate 2D variables from 3D grid:
    ModelVx=md.results.StressbalanceSolution.Vx(end-numvertices2d+1:end);
    ModelVy=md.results.StressbalanceSolution.Vy(end-numvertices2d+1:end);
    ModelVel=md.results.StressbalanceSolution.Vel(end-numvertices2d+1:end);
    ObsVx=md.inversion.vx_obs(end-numvertices2d+1:end);
    ObsVy=md.inversion.vy_obs(end-numvertices2d+1:end);
    ObsVel=md.inversion.vel_obs(end-numvertices2d+1:end);
    IceLevelSet=md.mask.ice_levelset(1:numvertices2d);
    OceanLevelSet=md.mask.ocean_levelset(1:numvertices2d);
    if includedrag
        if exist('CoeffConversionParameters','var')
            ModelDragCoeff=md.friction.coefficient(1:numvertices2d);
            OriginalDragCoeff=OriginalModel.md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d);
        else
            ModelDragCoeff=md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d);
            OriginalDragCoeff=md.friction.coefficient(1:numvertices2d);
        end
        if ~strcmp(nsource,'none')
            EffectivePressure_vertices=md.friction.effective_pressure(1:numvertices2d);
        end
    end
    if includerheology
        ModelBbar=md.results.StressbalanceSolution.MaterialsRheologyBbar(1:numvertices2d);
        OriginalBbar=md.materials.rheology_B(1:numvertices2d);
    end
    CostFuncCoeffs=md.inversion.cost_functions_coefficients(end-numvertices2d+1:end,:);
    
    %toc
    %return
    
    % Assign area term:
    if includedrag && includerheology
        thisarea=totalicearea;
    elseif includedrag
        thisarea=totalgroundedarea;
    else
        ActiveMask_bbar_elements=mean(ActiveMask_bbar(elements2d),2);
        thisarea=sum(Area.*ActiveMask_bbar_elements);
    end
    % Pre-allocate:
    RawMisfits=zeros(numvertices2d,length(costfunctions));
    NormMisfits=zeros(numvertices2d,length(costfunctions));
    % Loop through each component:
    for ii=1:length(costfunctions)
        % Check component:
        if costfunctions(ii)==101 % vel absolute misfit
            RawMisfits(:,ii)=.5*((ModelVx-ObsVx).^2+(ModelVy-ObsVy).^2);
            NormMisfits(:,ii)=thisarea*(CostFuncCoeffs(:,ii)/intrinsicweights(ii)).*RawMisfits(:,ii)/(secondsperyear^2);
        elseif costfunctions(ii)==102  % vel relative misfit
            % I HAVENT DONE THIS PART
        elseif costfunctions(ii)==102.1  % homebrew vel relative misfit
            RawMisfits(:,ii)=.5*((ModelVx-ObsVx).^2+(ModelVy-ObsVy).^2)./(ObsVel.^2+velstabilizer^2);
            NormMisfits(:,ii)=thisarea*((ObsVel/secondsperyear).^2+(velstabilizer/secondsperyear)^2).*(CostFuncCoeffs(:,ii)/intrinsicweights(ii)).*RawMisfits(:,ii);
        elseif costfunctions(ii)==103   % vel log misfit
            RawMisfits(:,ii)=log(ModelVel./ObsVel).^2;
            NormMisfits(:,ii)=thisarea*4*((issmmeanvel/secondsperyear)^2)*(CostFuncCoeffs(:,ii)/intrinsicweights(ii)).*RawMisfits(:,ii);
        elseif costfunctions(ii)==501   % drag coeff gradient
            % Compute gradients:
            DragCoeffGradX_elements=sum(BasisGradX_elements.*ModelDragCoeff(elements2d),2);
            DragCoeffGradY_elements=sum(BasisGradY_elements.*ModelDragCoeff(elements2d),2);
            % Compute (squared) magnitude:
            DragCoeffGradMag_elements=DragCoeffGradX_elements.^2+DragCoeffGradY_elements.^2;
            % Average from elements onto vertices:
            % pre-allocate:
            GradMag_vertices=zeros(numvertices2d,1);
            Denom_vertices=zeros(numvertices2d,1);
            % loop through elements:
            for jj=1:numelements2d
                % Add to the running sums:
                GradMag_vertices(elements2d(jj,:))=GradMag_vertices(elements2d(jj,:))+DragCoeffGradMag_elements(jj);
                Denom_vertices(elements2d(jj,:))=Denom_vertices(elements2d(jj,:))+1;
            end
            % finish computing means:
            GradMag_vertices=GradMag_vertices./Denom_vertices;
            % Compute misfit:
            RawMisfits(:,ii)=GradMag_vertices;
            NormMisfits(:,ii)=thisarea*.5*(CostFuncCoeffs(:,ii)/intrinsicweights(ii)).*RawMisfits(:,ii);
        elseif costfunctions(ii)==502  % bbar gradient
            % Compute gradients:
            BbarGradX_elements=sum(BasisGradX_elements.*ModelBbar(elements2d),2);
            BbarGradY_elements=sum(BasisGradY_elements.*ModelBbar(elements2d),2);
            % Compute (squared) magnitude:
            BbarGradMag_elements=BbarGradX_elements.^2+BbarGradY_elements.^2;
            % Average from elements onto vertices:
            % pre-allocate:
            GradMag_vertices=zeros(numvertices2d,1);
            Denom_vertices=zeros(numvertices2d,1);
            % loop through elements:
            for jj=1:numelements2d
                % Add to the running sums:
                GradMag_vertices(elements2d(jj,:))=GradMag_vertices(elements2d(jj,:))+BbarGradMag_elements(jj);
                Denom_vertices(elements2d(jj,:))=Denom_vertices(elements2d(jj,:))+1;
            end
            % finish computing means:
            GradMag_vertices=GradMag_vertices./Denom_vertices;
            % Compute misfit:
            RawMisfits(:,ii)=GradMag_vertices;
            NormMisfits(:,ii)=thisarea*.5*(CostFuncCoeffs(:,ii)/intrinsicweights(ii)).*RawMisfits(:,ii);
        end
    end
    % Compute additional misfits:
    VelMagMisfit=ModelVel-ObsVel;  % m/yr
    VelMagNormMisfit=(ModelVel-ObsVel)./ObsVel;
    
    %% Compute histograms:
    
    % Note: I want area-weighted probability density, so I cannot use
    % Matlab's built-in histcounts function.  The output for each should
    % have units of area fraction per d(log(metric)).  The only exception
    % is for the two magnitude misfits (these aren't included in the cost
    % function, just my own metrics), for which the denominator is just
    % d(metric), ie, linear scale instead of log scale.
    
    % Average misfits from vertices to elements:
    RawMisfits_elements=zeros(numelements2d,length(costfunctions));
    NormMisfits_elements=zeros(numelements2d,length(costfunctions));
    for ii=1:length(costfunctions)
        RawMisfits_elements(:,ii)=mean(reshape(RawMisfits(elements2d,ii),[],3),2);
        NormMisfits_elements(:,ii)=mean(reshape(NormMisfits(elements2d,ii),[],3),2);
    end
    RawMagMisfit_elements=mean(VelMagMisfit(elements2d),2);
    NormMagMisfit_elements=mean(VelMagNormMisfit(elements2d),2);
    
    % Compute element weight based on the active masks:
    if includedrag && includerheology
        ElementWeight=mean(ActiveMask_bbar(elements2d)|ActiveMask_drag(elements2d),2);
    elseif includedrag
        ElementWeight=mean(ActiveMask_drag(elements2d),2);
    else
        ElementWeight=mean(ActiveMask_bbar(elements2d),2);
    end
    
    % Make bin edges:
    % pre-allocate:
    Edges_rawmisfits=zeros(numbins+1,length(costfunctions));
    Edges_normmisfits=zeros(numbins+1,length(costfunctions));
    % Generate mask for limits computation:
    if includedrag && includerheology
        MaskForLimits=repmat(sum(isnan(NormMisfits)|isinf(NormMisfits),2)==0&...
            isnan(VelMagNormMisfit)==0&isinf(VelMagNormMisfit)==0&ObsVel>minvelobs4hist&...
            GoodData_obs_vertices>.5&(ActiveMask_drag|ActiveMask_bbar),[1,length(costfunctions)]);
        for ii=1:length(costfunctions)
            MaskForLimits(:,ii)=MaskForLimits(:,ii)&NormMisfits(:,ii)~=0;
        end
    elseif includedrag
        MaskForLimits=repmat(sum(isnan(NormMisfits)|isinf(NormMisfits),2)==0&...
            isnan(VelMagNormMisfit)==0&isinf(VelMagNormMisfit)==0&ObsVel>minvelobs4hist&...
            GoodData_obs_vertices>.5&ActiveMask_drag,[1,length(costfunctions)]);
        for ii=1:length(costfunctions)
            MaskForLimits(:,ii)=MaskForLimits(:,ii)&NormMisfits(:,ii)~=0;
        end
    else
        MaskForLimits=repmat(sum(isnan(NormMisfits)|isinf(NormMisfits),2)==0&...
            isnan(VelMagNormMisfit)==0&isinf(VelMagNormMisfit)==0&ObsVel>minvelobs4hist&...
            GoodData_obs_vertices>.5&ActiveMask_bbar,[1,length(costfunctions)]);
        for ii=1:length(costfunctions)
            MaskForLimits(:,ii)=MaskForLimits(:,ii)&NormMisfits(:,ii)~=0;
        end
    end
    % Loop through misfits:
    for ii=1:length(costfunctions)
        % Compute medians, use for min/max:
        thisrawmedian=median(RawMisfits(MaskForLimits(:,ii),ii));
        thisnormmedian=median(NormMisfits(MaskForLimits(:,ii),ii));
        thisrawmin=10.^(log10(thisrawmedian)-4);
        thisrawmax=10.^(log10(thisrawmedian)+4);
        thisnormmin=10.^(log10(thisnormmedian)-4);
        thisnormmax=10.^(log10(thisnormmedian)+4);
        % Identify min and max:
        %thisrawmin=min(RawMisfits(MaskForLimits,ii));
        %thisrawmax=max(RawMisfits(MaskForLimits,ii));
        %thisnormmin=min(NormMisfits(MaskForLimits,ii));
        %thisnormmax=max(NormMisfits(MaskForLimits,ii));
        % Compute bin edges on log scale:
        Edges_rawmisfits(:,ii)=exp(linspace(log(thisrawmin),log(thisrawmax),numbins+1)');
        Edges_normmisfits(:,ii)=exp(linspace(log(thisnormmin),log(thisnormmax),numbins+1)');
    end
    % Compute bin edges for mag misfits:
    %thisrawmin=min(VelMagMisfit(MaskForLimits));
    %thisrawmax=max(VelMagMisfit(MaskForLimits));
    %thisnormmin=min(VelMagNormMisfit(MaskForLimits));
    %thisnormmax=max(VelMagNormMisfit(MaskForLimits));
    %thisrawmin=rawmaglims(1);
    %thisrawmax=rawmaglims(2);
    thisnormmin=normmaglims(1);
    thisnormmax=normmaglims(2);
    rawmagrms=sqrt(sum((RawMagMisfit_elements.^2).*Area.*ElementWeight)/sum(Area.*ElementWeight));
    Edges_rawmagmisfit=rawmagrms*linspace(-3,3,numbins+1)';
    Edges_normmagmisfit=linspace(thisnormmin,thisnormmax,numbins+1)';
    
    % Compute probability density:
    % Pre-allocate:
    ProbDensity_rawmisfits=zeros(numbins,length(costfunctions));
    ProbDensity_normmisfits=zeros(numbins,length(costfunctions));
    ProbDensity_rawmagmisfit=zeros(numbins,1);
    ProbDensity_normmagmisfit=zeros(numbins,1);
    % Loop through bins:
    for ii=1:numbins
        % Add up area in this bin:
        for jj=1:length(costfunctions)
            ProbDensity_rawmisfits(ii,jj)=sum(Area(RawMisfits_elements(:,jj)>=Edges_rawmisfits(ii,jj)&RawMisfits_elements(:,jj)<Edges_rawmisfits(ii+1,jj)).*...
                ElementWeight(RawMisfits_elements(:,jj)>=Edges_rawmisfits(ii,jj)&RawMisfits_elements(:,jj)<Edges_rawmisfits(ii+1,jj)));
            ProbDensity_normmisfits(ii,jj)=sum(Area(NormMisfits_elements(:,jj)>=Edges_normmisfits(ii,jj)&NormMisfits_elements(:,jj)<Edges_normmisfits(ii+1,jj)).*...
                ElementWeight(NormMisfits_elements(:,jj)>=Edges_normmisfits(ii,jj)&NormMisfits_elements(:,jj)<Edges_normmisfits(ii+1,jj)));
        end
        ProbDensity_rawmagmisfit(ii)=sum(Area(RawMagMisfit_elements>=Edges_rawmagmisfit(ii)&RawMagMisfit_elements<Edges_rawmagmisfit(ii+1)).*...
            ElementWeight(RawMagMisfit_elements>=Edges_rawmagmisfit(ii)&RawMagMisfit_elements<Edges_rawmagmisfit(ii+1)));
        ProbDensity_normmagmisfit(ii)=sum(Area(NormMagMisfit_elements>=Edges_normmagmisfit(ii)&NormMagMisfit_elements<Edges_normmagmisfit(ii+1)).*...
            ElementWeight(NormMagMisfit_elements>=Edges_normmagmisfit(ii)&NormMagMisfit_elements<Edges_normmagmisfit(ii+1)));
    end
    % Normalize area counts:
    ProbDensity_rawmisfits=ProbDensity_rawmisfits./repmat(thisarea*log(Edges_rawmisfits(2,:)./Edges_rawmisfits(1,:)),[numbins,1]);
    ProbDensity_normmisfits=ProbDensity_normmisfits./repmat(thisarea*log(Edges_normmisfits(2,:)./Edges_normmisfits(1,:)),[numbins,1]);
    ProbDensity_rawmagmisfit=ProbDensity_rawmagmisfit/(thisarea*(Edges_rawmagmisfit(2)-Edges_rawmagmisfit(1)));
    ProbDensity_normmagmisfit=ProbDensity_normmagmisfit/(thisarea*(Edges_normmagmisfit(2)-Edges_normmagmisfit(1)));
    
    %% Make histogram figure:
    
    % Compute number of subplots:
    numtotalplots=length(costfunctions)+1+2;
    numvertplots=ceil(numtotalplots/2);
    numhorzplots=2;
    
    % Compute subplot dimensions:
    % Compute available space:
    availablehorzspace=1-numhorzplots*(horznotextbuffer_hist+horzaxesbuffer_hist+horzaxeslabelbuffer_hist);
    availablevertspace=1-numvertplots*(titlebuffer_hist+vertaxesbuffer_hist+vertaxeslabelbuffer_hist);
    % Compute boxes:
    Boxes=cell(numvertplots,numhorzplots);
    for ii=1:numvertplots
        for jj=1:numhorzplots
            Boxes{numvertplots-ii+1,jj}=[jj*(horzaxesbuffer_hist+horzaxeslabelbuffer_hist)+(jj-1)*(horznotextbuffer_hist+availablehorzspace/numhorzplots),...
                ii*(vertaxesbuffer_hist+vertaxeslabelbuffer_hist)+(ii-1)*(titlebuffer_hist+availablevertspace/numvertplots),...
                availablehorzspace/numhorzplots,availablevertspace/numvertplots];
        end
    end
    
    % Make figure:
    figure
    
    % All normalized misfits:
    % Make subplot:
    subplot('Position',Boxes{1,1})
    % Loop through misfits:
    for ii=1:length(costfunctions)
        % Make this histogram:
        histogram('BinEdges',Edges_normmisfits(:,ii)','BinCounts',ProbDensity_normmisfits(:,ii),'DisplayStyle','stairs')
        hold on
    end
    % Set up axes:
    set(gca,'XScale','log')
    xmedian=round(log10(median(sum(NormMisfits(max(MaskForLimits,[],2),:),2))));
    thesexlims=10.^(xmedian+[-4,4]);
    xlim(thesexlims)
    set(gca,'XTick',10.^[ceil(log10(thesexlims(1))):floor(log10(thesexlims(2)))])
    xlabel('Normalized Misfit')
    ylabel('Area-weighted probability density')
    % Make legend:
    legendtext=cell(length(costfunctions),1);
    for ii=1:length(costfunctions)
        if costfunctions(ii)==101
            legendtext{ii}='Velocity absolute misfit';
        elseif costfunctions(ii)==102 || costfunctions(ii)==102.1
            legendtext{ii}='Velocity relative misfit';
        elseif costfunctions(ii)==103
            legendtext{ii}='Velocity logarithmic misfit';
        elseif costfunctions(ii)==501
            legendtext{ii}='Drag regularization';
        elseif costfunctions(ii)==502
            legendtext{ii}='Rheology regularization';
        end
    end
    legend(legendtext,'location','NorthEast')
    % Make a title:
    title({'a) All Normalized Misfits';'(included in cost function)'})
    
    % Loop through absolute misfits:
    for ii=1:length(costfunctions)
        % Make subplot:
        subplot('Position',Boxes{ceil((ii+1)/2),rem(ii,2)+1})
        % Make histogram:
        histogram('BinEdges',Edges_rawmisfits(:,ii)','BinCounts',ProbDensity_rawmisfits(:,ii),'DisplayStyle','stairs','EdgeColor','k')
        % Set up axes:
        set(gca,'XScale','log')
        xmedian=round(log10(median(RawMisfits(MaskForLimits(:,ii),ii))));
        thesexlims=10.^(xmedian+[-4,4]);
        xlim(thesexlims)
        set(gca,'XTick',10.^[ceil(log10(thesexlims(1))):floor(log10(thesexlims(2)))])
        ylabel('Area-weighted probability density')
        % Check what misfit this is, set x-axis label:
        if costfunctions(ii)==101
            xlabel('Velocity absolute misfit, (m/yr)^2');
        elseif costfunctions(ii)==102 || costfunctions(ii)==102.1
            xlabel('Velocity relative misfit, unitless');
        elseif costfunctions(ii)==103
            xlabel('Velocity logarithmic misfit, ln(v/v_{obs})^2');
        elseif costfunctions(ii)==501
            if m==1
                if ~strcmp(nsource,'none')
                    if r==1
                        xlabel('Drag coefficient squared gradient, (m/s)^{-1}m^{-2}');
                    else
                        xlabel('Drag coefficient squared gradient, Pa^{',num2str(1-r),'}(m/s)^{-1}m^{-2}');
                    end
                else
                    xlabel('Drag coefficient squared gradient, Pa(m/s)^{-1}m^{-2}');
                end
            else
                if ~strcmp(nsource,'none')
                    if r==1
                        xlabel(['Drag coefficient squared gradient, (m/s)^{-1/',num2str(m),'}m^{-2}']);
                    else
                        xlabel(['Drag coefficient squared gradient, Pa^{',num2str(1-r),'}(m/s)^{-1/',num2str(m),'}m^{-2}']);
                    end
                else
                    xlabel(['Drag coefficient squared gradient, Pa(m/s)^{-1/',num2str(m),'}m^{-2}']);
                end
            end
        elseif costfunctions(ii)==502
            xlabel(['Bbar squared gradient, Pa^{2}s^{2/',num2str(n),'}m^{-2}']);
        end
        % Make title:
        title([char(97+ii),') Cost function component ',num2str(costfunctions(ii))])
    end
    
    % Velocity absolute magnitude misfit:
    % Make subplot:
    subplot('Position',Boxes{ceil((length(costfunctions)+2)/2),rem(length(costfunctions)+1,2)+1})
    % Make histogram:
    histogram('BinEdges',Edges_rawmagmisfit','BinCounts',ProbDensity_rawmagmisfit,'DisplayStyle','stairs','EdgeColor','k')
    % Set up axes:
    theseylims=get(gca,'Ylim');
    hold on
    plot([0,0],theseylims,'--k')
    %xlim(rawmaglims)
    xlabel({'Velocity magnitude difference,';'|v|-|v_{obs}| (m/yr)'})
    ylabel('Area-weighted probability density')
    % Make title:
    title({[char(97+length(costfunctions)+1),') Additional fit metric: absolute magnitude difference'];'(not included in cost function)'})
    
    % Velocity relative magnitude misfit:
    % Make subplot:
    subplot('Position',Boxes{ceil((length(costfunctions)+3)/2),rem(length(costfunctions)+2,2)+1})
    % Make histogram:
    histogram('BinEdges',Edges_normmagmisfit','BinCounts',ProbDensity_normmagmisfit,'DisplayStyle','stairs','EdgeColor','k')
    % Set up axes:
    theseylims=get(gca,'Ylim');
    hold on
    plot([0,0],theseylims,'--k')
    xlim(normmaglims)
    xlabel({'Velocity normalized magnitude difference,';'(|v|-|v_{obs}|)/|v_{obs}| (unitless)'})
    ylabel('Area-weighted probability density')
    % Make title:
    title({[char(97+length(costfunctions)+2),') Additional fit metric: relative magnitude difference'];'(not included in cost function)'})
    
    % Save figure:
    thispagesize=[numhorzplots*widthperhist,numvertplots*heightperhist];
    if dolcurve && exist('Lambda','var')
        thisfigname=[outputfile(1:slashnum(end)),thisname,'_histograms.png'];
    else
        thisfigname=[outputfile(1:end-4),'_histograms.png'];
    end
    set(gcf,'PaperSize',thispagesize)
    set(gcf,'PaperPosition',[0,0,thispagesize])
    print('-dpng',thisfigname,['-r',num2str(resolution)])
    
    % Close figure:
    close
    
    %% Make map figure:
    
    % Make figure:
    figure
    
    % Compute number of subplots:
    numtotalplots=4+4*includedrag+(includedrag&~strcmp(nsource,'none'))+4*includerheology+length(costfunctions)+1;
    numvertplots=ceil(numtotalplots/4);
    numhorzplots=4;
    
    % Compute subplot dimensions:
    % Compute available space:
    availablehorzspace=1-numhorzplots*(horznotextbuffer_map+horzaxesbuffer_map+horzaxeslabelbuffer_map);
    availablevertspace=1-numvertplots*(titlebuffer_map+vertaxesbuffer_map+vertaxeslabelbuffer_map);
    % Compute boxes:
    Boxes=cell(numvertplots,numhorzplots);
    for ii=1:numvertplots
        for jj=1:numhorzplots
            Boxes{numvertplots-ii+1,jj}=[jj*(horzaxesbuffer_map+horzaxeslabelbuffer_map)+(jj-1)*(horznotextbuffer_map+availablehorzspace/numhorzplots),...
                ii*(vertaxesbuffer_map+vertaxeslabelbuffer_map)+(ii-1)*(titlebuffer_map+availablevertspace/numvertplots),...
                availablehorzspace/numhorzplots,availablevertspace/numvertplots];
        end
    end
    % Make boxes transpose:
    Boxes=Boxes';
    
    % Make the faded jet colormap:
    FadedJet=colormap('jet');
    FadedJet(3*round(size(FadedJet,1)/8)+1:round(size(FadedJet,1)/2),1)=linspace(0,1,round(size(FadedJet,1)/8))';
    FadedJet(3*round(size(FadedJet,1)/8)+1:round(size(FadedJet,1)/2),3)=1;
    FadedJet(round(size(FadedJet,1)/2)+1:5*round(size(FadedJet,1)/8),1)=1;
    FadedJet(round(size(FadedJet,1)/2)+1:5*round(size(FadedJet,1)/8),3)=linspace(1,0,round(size(FadedJet,1)/8))';
    
    % Make the clipped jet colormap:
    ClippedJet=colormap('jet');
    ClippedJet=ClippedJet(round(size(ClippedJet,1)/8):end-round(size(ClippedJet,1)/8),:);
    
    % Get regular jet:
    RegularJet=colormap('jet');
    
    % Initialize plot counter:
    thisplot=1;
    
    % Observed velocity:
    % Make subplot:
    subplot('Position',Boxes{thisplot})
    % Plot velocity:
    patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d,'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(ObsVel))
    caxis([-0.5,3.5])
    xlim(xlims/1000)
    ylim(ylims/1000)
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','log_{10}(m/yr)')
    title('a) Observed Velocity Magnitude')
    % Advance plot counter:
    thisplot=thisplot+1;
    
    % Model velocity:
    % Make subplot:
    subplot('Position',Boxes{thisplot})
    % Plot velocity:
    patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d,'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(ModelVel))
    caxis([-0.5,3.5])
    xlim(xlims/1000)
    ylim(ylims/1000)
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','log_{10}(m/yr)')
    title('b) Modelled Velocity Magnitude')
    % Advance plot counter:
    thisplot=thisplot+1;
    
    % Check whether to plot drag:
    if includedrag
        
        % Identify elements to plot drag for:
        TheseElements=IceAreaFraction_elements>=.5&OceanAreaFraction_elements<=.5;
        
        % Drag coefficient:
        % Make subplot:
        subplot('Position',Boxes{thisplot})
        % Plot velocity:
        patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(ModelDragCoeff))
        caxis(log10(prctile(ModelDragCoeff(ActiveMask_drag),[.1,99.9])))
        xlim(xlims/1000)
        ylim(ylims/1000)
        set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
        xlabel('Easting (km)')
        ylabel('Northing (km)')
        hc=colorbar;
        if ~strcmp(nsource,'none')
            if r==1
                set(get(hc,'Ylabel'),'String',['log_{10}((m/s)^{-1/',num2str(2*m),'})'])
            else
                set(get(hc,'Ylabel'),'String',['log_{10}(Pa^{',num2str((1-r)/2),'}(m/s)^{-1/',num2str(2*m),'})'])
            end
        else
            set(get(hc,'Ylabel'),'String',['log_{10}(Pa^{1/2}(m/s)^{-1/',num2str(2*m),'})'])
        end
        if exist('CoeffConversionParameters','var')
            title([char(96+thisplot),') Final Drag Coefficient (after conversion)'])
        else
            title([char(96+thisplot),') Final Drag Coefficient (after inversion)'])
        end
        % Advance plot counter:
        thisplot=thisplot+1;
        
        % Compute basal drag:
        if numlayers>1 && ~strcmp(flowequation,'SSA')
            Drag=(ModelDragCoeff.^2).*((md.results.StressbalanceSolution.Vel(1:numvertices2d)/secondsperyear).^(1/m));
        else
            Drag=(ModelDragCoeff.^2).*((md.results.StressbalanceSolution.Vel/secondsperyear).^(1/m));
        end
        if ~strcmp(nsource,'none')
            Drag=Drag.*(EffectivePressure_vertices.^r);
        end
        % Set drag to zero under floating shelf and in ice-free areas:
        Drag(OceanLevelSet_vertices<0|IceLevelSet_vertices>0)=0;
        
        % Drag:
        % Make subplot:
        subplot('Position',Boxes{thisplot})
        % Plot velocity:
        patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Drag/1000)
        caxis(draglims)
        xlim(xlims/1000)
        ylim(ylims/1000)
        set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
        xlabel('Easting (km)')
        ylabel('Northing (km)')
        hc=colorbar;
        set(get(hc,'Ylabel'),'String','kPa')
        title([char(96+thisplot),') Basal Drag'])
        % Advance plot counter:
        thisplot=thisplot+1;
        
        % Original Drag coefficient:
        % Make subplot:
        subplot('Position',Boxes{thisplot})
        % Plot velocity:
        patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(OriginalDragCoeff))
        if exist('CoeffConversionParameters','var')
            caxis(log10(prctile(OriginalDragCoeff(ActiveMask_drag),[.1,99.9])))
        else
            caxis(log10(prctile(ModelDragCoeff(ActiveMask_drag),[.1,99.9])))
        end
        xlim(xlims/1000)
        ylim(ylims/1000)
        set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
        xlabel('Easting (km)')
        ylabel('Northing (km)')
        hc=colorbar;
        if exist('CoeffConversionParameters','var')
            thism=CoeffConversionParameters.oldm;
            thisr=CoeffConversionParameters.oldr;
            thisnsource=CoeffConversionParameters.oldnsource;
        else
            thism=m;
            thisr=r;
            thisnsource=nsource;
        end
        if ~strcmp(thisnsource,'none')
            if thisr==1
                set(get(hc,'Ylabel'),'String',['log_{10}((m/s)^{-1/',num2str(2*thism),'})'])
            else
                set(get(hc,'Ylabel'),'String',['log_{10}(Pa^{',num2str((1-thisr)/2),'}(m/s)^{-1/',num2str(2*thism),'})'])
            end
        else
            set(get(hc,'Ylabel'),'String',['log_{10}(Pa^{1/2}(m/s)^{-1/',num2str(2*thism),'})'])
        end
        if exist('CoeffConversionParameters','var')
            title([char(96+thisplot),') Initial Drag Coefficient (before conversion)'])
        else
            title([char(96+thisplot),') Initial Drag Coefficient (before inversion)'])
        end
        % Advance plot counter:
        thisplot=thisplot+1;
        
        % Change in drag coefficient:
        % Make subplot:
        subplot('Position',Boxes{thisplot})
        % Plot variable
        patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(ModelDragCoeff./OriginalDragCoeff))
        if exist('CoeffConversionParameters','var')
            originalmin=min(OriginalDragCoeff(IceLevelSet<0&OceanLevelSet>0));
            dragcoeffratiolims=dragcoeffratiolims+...
                sum(Area.*mean(ActiveMask_drag(elements2d),2).*mean(log10(ModelDragCoeff(elements2d)),2))/totalgroundedarea-...
                sum(Area.*mean(ActiveMask_drag(elements2d),2).*mean(log10(max(originalmin,OriginalDragCoeff(elements2d))),2))/totalgroundedarea;
        end
        caxis(dragcoeffratiolims)
        clength=size(FadedJet,1);
        if exist('CoeffConversionParameters','var')
            ThisFadedJet=FadedJet;
        else
            ThisFadedJet=interp1(linspace(-1,1,clength)',FadedJet,[linspace(-1,0,round(clength*(-dragcoeffratiolims(1)/(dragcoeffratiolims(2)-dragcoeffratiolims(1)))))';...
                linspace(0,1,clength-round(clength*(-dragcoeffratiolims(1)/(dragcoeffratiolims(2)-dragcoeffratiolims(1)))))'],'linear');
        end
        colormap(gca,ThisFadedJet)
        xlim(xlims/1000)
        ylim(ylims/1000)
        set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
        xlabel('Easting (km)')
        ylabel('Northing (km)')
        hc=colorbar;
        set(get(hc,'Ylabel'),'String',['log_{10}(C_{final}/C_{initial})'])
        title([char(96+thisplot),') Change in Drag Coefficient'])
        % Advance plot counter:
        thisplot=thisplot+1;
        
        % Effective Pressure:
        if ~strcmp(nsource,'none')
            % Make subplot:
            subplot('Position',Boxes{thisplot})
            % Plot velocity:
            patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',EffectivePressure_vertices/1000)
            caxis(prctile(EffectivePressure_vertices(ActiveMask_drag),[.1,99.9])/1e3)
            xlim(xlims/1000)
            ylim(ylims/1000)
            set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
            xlabel('Easting (km)')
            ylabel('Northing (km)')
            hc=colorbar;
            set(get(hc,'Ylabel'),'String','kPa')
            title([char(96+thisplot),') Effective Pressure'])
            % Advance plot counter:
            thisplot=thisplot+1;
        end
        
    end
    
    % Check whether to plot rheology:
    if includerheology
        
        % Determine which elements to plot rheology for:
        ActiveMask_bbar_elements=mean(ActiveMask_bbar(elements2d),2);
        TheseElements=ActiveMask_bbar_elements>=.5;
        
        % Bbar coefficient:
        % Make subplot:
        subplot('Position',Boxes{thisplot})
        % Plot Bbar:
        patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(ModelBbar))
        caxis(log10(prctile(ModelBbar(ActiveMask_bbar),[.1,99.9])))
        xlim(xlims/1000)
        ylim(ylims/1000)
        set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
        xlabel('Easting (km)')
        ylabel('Northing (km)')
        hc=colorbar;
        set(get(hc,'Ylabel'),'String',['log_{10}(Pa s^{1/',num2str(n),'})'])
        %set(get(hc,'Ylabel'),'String',['Pa s^{1/',num2str(n),'}'])
        title({[char(96+thisplot),') Final Vertically Averaged Rheology'];'(Bbar after inversion)'})
        % Advance plot counter:
        thisplot=thisplot+1;
        
        % Assign viscosity result to materials section:
        md.materials.rheology_B=md.results.StressbalanceSolution.MaterialsRheologyBbar;
        % Compute viscosity:
        Viscosity_elements=ComputeViscosity(md);
        % Compute vertically averaged viscosity:
        if numlayers>1 && ~strcmp(flowequation,'SSA')
            DZhat=(md.mesh.z(1+numvertices2d*(1:numlayers-1)')-md.mesh.z(1+numvertices2d*(0:numlayers-2)'))/md.geometry.thickness(1);
            ViscosityBar_elements2d=sum(reshape(Viscosity_elements,[numelements2d,numlayers-1]).*repmat(DZhat',[numelements2d,1]),2);
        else
            ViscosityBar_elements2d=Viscosity_elements;
        end
        
        % Viscosity:
        % Make subplot:
        subplot('Position',Boxes{thisplot})
        % Plot viscosity:
        patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d,'EdgeColor','none','FaceColor','flat','FaceVertexCData',log10(ViscosityBar_elements2d))
        caxis(viscositylims)
        xlim(xlims/1000)
        ylim(ylims/1000)
        set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
        xlabel('Easting (km)')
        ylabel('Northing (km)')
        hc=colorbar;
        set(get(hc,'Ylabel'),'String','log_{10}(Pa*s)')
        title([char(96+thisplot),') Vertically Averaged Viscosity'])
        % Advance plot counter:
        thisplot=thisplot+1;
        
        % Original Bbar coefficient:
        % Make subplot:
        subplot('Position',Boxes{thisplot})
        % Plot Bbar:
        patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(OriginalBbar))
        caxis(log10(prctile(ModelBbar(ActiveMask_bbar),[.1,99.9])))
        xlim(xlims/1000)
        ylim(ylims/1000)
        set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
        xlabel('Easting (km)')
        ylabel('Northing (km)')
        hc=colorbar;
        set(get(hc,'Ylabel'),'String',['log_{10}(Pa s^{1/',num2str(n),'})'])
        %set(get(hc,'Ylabel'),'String',['Pa s^{1/',num2str(n),'}'])
        title({[char(96+thisplot),') Initial Vertically Averaged Rheology'];'(Bbar before inversion)'})
        % Advance plot counter:
        thisplot=thisplot+1;
        
        % Change in Bbar:
        % Make subplot:
        subplot('Position',Boxes{thisplot})
        % Plot variable
        patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(ModelBbar./OriginalBbar))
        caxis(bbarratiolims)
        clength=size(FadedJet,1);
        ThisFadedJet=interp1(linspace(-1,1,clength)',FadedJet,[linspace(-1,0,round(clength*(-bbarratiolims(1)/(bbarratiolims(2)-bbarratiolims(1)))))';...
            linspace(0,1,clength-round(clength*(-bbarratiolims(1)/(bbarratiolims(2)-bbarratiolims(1)))))'],'linear');
        colormap(gca,ThisFadedJet)
        xlim(xlims/1000)
        ylim(ylims/1000)
        set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
        xlabel('Easting (km)')
        ylabel('Northing (km)')
        hc=colorbar;
        set(get(hc,'Ylabel'),'String','log_{10}(B_{final}/B_{initial})')
        title([char(96+thisplot),') Change in Bbar Coefficient'])
        % Advance plot counter:
        thisplot=thisplot+1;
        
    end
    
    % Identify elements to plot misfit for:
    if includedrag && includerheology==0
        TheseElements=IceAreaFraction_elements>=.5&OceanAreaFraction_elements<=.5;
    else
        TheseElements=IceAreaFraction_elements>=.5;
    end
    
    % Loop through cost function components:
    for ii=1:length(costfunctions)
        % Make subplot:
        %subplot('Position',Boxes{2+ceil((ii)/2),rem(ii+1,2)+1})
        subplot('Position',Boxes{thisplot})
        % Plot normalized misfit:
        patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(NormMisfits(:,ii)))
        cmedian=round(log10(median(NormMisfits(MaskForLimits(:,ii),ii))));
        theseclims=cmedian+[-3,3];
        caxis(theseclims)
        xlim(xlims/1000)
        ylim(ylims/1000)
        set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
        xlabel('Easting (km)')
        ylabel('Northing (km)')
        hc=colorbar;
        set(get(hc,'Ylabel'),'String','log_{10}(normalized misfit)')
        if costfunctions(ii)==101
            title({[char(96+thisplot),') Cost function component:'];'Absolute velocity misfit'})
        elseif costfunctions(ii)==102 || costfunctions(ii)==102.1
            title({[char(96+thisplot),') Cost function component:'];'Relative velocity misfit'})
        elseif costfunctions(ii)==103
            title({[char(96+thisplot),') Cost function component:'];'Logarithmic velocity misfit'})
        elseif costfunctions(ii)==501
            title({[char(96+thisplot),') Cost function component:'];'Drag regularizer'})
        elseif costfunctions(ii)==502
            title({[char(96+thisplot),') Cost function component:'];'Rheology regularizer'})
        end
        % Advance plot counter:
        thisplot=thisplot+1;
    end
    
    % Plot total cost function:
    % Make subplot:
    %subplot('Position',Boxes{2+ceil((length(costfunctions)+1)/2),rem(length(costfunctions)+2,2)+1})
    subplot('Position',Boxes{thisplot})
    % Plot normalized misfit:
    patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(TheseElements,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(sum(NormMisfits.*repmat(intrinsicweights,[numvertices2d,1]),2)))
    cmedian=round(log10(median(sum(NormMisfits(max(MaskForLimits,[],2),:),2))));
    theseclims=cmedian+[-3,3];
    caxis(theseclims)
    xlim(xlims/1000)
    ylim(ylims/1000)
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','log_{10}(normalized misfit)')
    title([char(96+thisplot),') Total cost function'])
    % Advance plot counter:
    thisplot=thisplot+1;
    
    % Plot velocity magnitude misfit:
    % Make subplot:
    %subplot('Position',Boxes{2+ceil((length(costfunctions)+2)/2),rem(length(costfunctions)+3,2)+1})
    subplot('Position',Boxes{thisplot})
    % Plot normalized misfit:
    patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d,'EdgeColor','none','FaceColor','interp','FaceVertexCData',VelMagMisfit)
    caxis(rawmaglims)
    colormap(gca,FadedJet)
    xlim(xlims/1000)
    ylim(ylims/1000)
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','|v|-|v_{obs}| (m/yr)')
    title({[char(96+thisplot),') Velocity magnitude misfit'];'(not in cost function)'})
    % Advance plot counter:
    thisplot=thisplot+1;
    
    % Plot normalized velocity magnitude misfit:
    % Make subplot:
    %subplot('Position',Boxes{2+ceil((length(costfunctions)+3)/2),rem(length(costfunctions)+4,2)+1})
    subplot('Position',Boxes{thisplot})
    % Plot normalized misfit:
    patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d,'EdgeColor','none','FaceColor','interp','FaceVertexCData',VelMagNormMisfit)
    caxis(normmaglims)
    colormap(gca,FadedJet)
    xlim(xlims/1000)
    ylim(ylims/1000)
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)')
    ylabel('Northing (km)')
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','(|v|-|v_{obs}|)/|v_{obs}| (unitless)')
    title({[char(96+thisplot),') Velocity normalized magnitude misfit'];'(not in cost function)'})
    % Advance plot counter:
    thisplot=thisplot+1;
    
    
    % Save figure:
    thispagesize=[numhorzplots*widthpermap,numvertplots*heightpermap];
    if dolcurve && exist('Lambda','var')
        thisfigname=[outputfile(1:slashnum(end)),thisname,'_maps.png'];
    else
        thisfigname=[outputfile(1:end-4),'_maps.png'];
    end
    set(gcf,'PaperSize',thispagesize)
    set(gcf,'PaperPosition',[0,0,thispagesize])
    print('-dpng',thisfigname,['-r',num2str(resolution)])
    
    % Close figure:
    close
    
    %% Misfit progression figure:
    
    if ~exist('CoeffConversionParameters','var')
        
        % Make figure:
        figure
        % Plot variables:
        plot(md.results.StressbalanceSolution.J(:,end),'k')
        hold on
        for ii=1:length(costfunctions)
            plot(md.results.StressbalanceSolution.J(:,ii))
        end
        % Fix up axes etc:
        ylim(10.^[floor(log10(min(md.results.StressbalanceSolution.J(:)))),ceil(log10(max(md.results.StressbalanceSolution.J(:))))])
        set(gca,'YScale','log')
        xlim([0,size(md.results.StressbalanceSolution.J,1)])
        xlabel('Iteration')
        ylabel('Cost Function')
        title('Misfit Progression During Inversion')
        % Make legend:
        legendtext=cell(length(costfunctions)+1,1);
        legendtext{1}='Total';
        for ii=1:length(costfunctions)
            legendtext{ii+1}=num2str(costfunctions(ii));
        end
        legend(legendtext,'location','NorthEast')
        
        % Save figure:
        thispagesize=[6,5];
        if dolcurve && exist('Lambda','var')
            thisfigname=[outputfile(1:slashnum(end)),thisname,'_misfitprogression.png'];
        else
            thisfigname=[outputfile(1:end-4),'_misfitprogression.png'];
        end
        set(gcf,'PaperSize',thispagesize)
        set(gcf,'PaperPosition',[0,0,thispagesize])
        print('-dpng',thisfigname,['-r',num2str(resolution)])
        
        % Close figure:
        close
        
    end
    
end
 %end
% Final display:
disp('Done!')
toc

%% Cutting Room Floor:


% %% Lambda Selection Via 2D Curve-Fitting:
% % Generate sorted list of non-outlier data points:
% LogDataCost_points=log(DataCost(~IsOutlier));
% LogRegCost_points=log(RegularizationCost(~IsOutlier));
% LogLambda_points=log(Lambda(~IsOutlier));
% [LogRegCost_points,SortedInd]=sort(LogRegCost_points,'descend');
% LogDataCost_points=LogDataCost_points(SortedInd);
% LogLambda_points=LogLambda_points(SortedInd);
%
% % Compute distance in log(cost)/log(cost) space:
% Dist_points=[0;cumsum(sqrt(diff(LogDataCost_points).^2+diff(LogRegCost_points).^2))];
%
% % Compute distance subsampling:
% Dist_path=linspace(0,max(Dist_points),numlambdasubsamples)';
%
% % Interpolate log(variables) onto the new distance subsampling:
% LogDataCost_path=interp1(Dist_points,LogDataCost_points,Dist_path,lambdainterpmethod);
% LogRegCost_path=interp1(Dist_points,LogRegCost_points,Dist_path,lambdainterpmethod);
% LogLambda_path=interp1(Dist_points,LogLambda_points,Dist_path,lambdainterpmethod);
%
% % Reinterpolate to approximately even spacing if necessary:
% if ~strcmp(lambdainterpmethod,'linear')
%     NewDist_path=linspace(0,max(Dist_path),numlambdasubsamples)';
%     LogDataCost_path=interp1(Dist_path,LogDataCost_path,NewDist_path,'linear');
%     LogRegCost_path=interp1(Dist_path,LogRegCost_path,NewDist_path,'linear');
%     LogLambda_path=interp1(Dist_path,LogLambda_path,NewDist_path,'linear');
%     Dist_path=[0;cumsum(sqrt(diff(LogDataCost_path).^2+diff(LogRegCost_path).^2))];
% end
%
% % Compute ddist:
% ddist=Dist_path(2)-Dist_path(1);
%
% % Compute minimum distance wavelength:
% mindistwavelength=.5*(mean(Dist_points(2:end)-Dist_points(1:end-1)));
%
% % Generate wavelength sampling:
% DistWavelengths=linspace(mindistwavelength,maxlambdawavelength,numlambdawavelengths);
%
% % Compute characteristic scales:
% overalldist=sqrt((LogDataCost_points(end)-LogDataCost_points(1)).^2+(LogRegCost_points(end)-LogRegCost_points(1)).^2);
% scattervarscale=(.5*overalldist)^2;
%
% % Pre-allocate:
% LogDataCost_curves=zeros(numlambdasubsamples,numlambdawavelengths);
% LogRegCost_curves=zeros(numlambdasubsamples,numlambdawavelengths);
% LogLambda_curves=zeros(numlambdasubsamples,numlambdawavelengths);
% Dist_curves=zeros(numlambdasubsamples,numlambdawavelengths);
% CurvatureRMS_curves=zeros(1,numlambdawavelengths);
% ScatterRMS_curves=zeros(1,numlambdawavelengths);
%
% % Loop through potential smoothings:
% for ii=1:numlambdawavelengths
%     % Check wavelength:
%     if DistWavelengths(ii)>0
%         % Generate expanded vectors:
%         numexpandedsamples=round(2*DistWavelengths(ii)/ddist);
%         DistExpanded=[ddist*linspace(-numexpandedsamples,-1,numexpandedsamples)';Dist_path;Dist_path(end)+ddist*linspace(1,numexpandedsamples,numexpandedsamples)'];
%         DataCostExpanded=[LogDataCost_path(1)*ones(numexpandedsamples,1);LogDataCost_path;LogDataCost_path(end)*ones(numexpandedsamples,1)];
%         RegCostExpanded=[LogRegCost_path(1)*ones(numexpandedsamples,1);LogRegCost_path;LogRegCost_path(end)*ones(numexpandedsamples,1)];
%         LambdaExpanded=[LogLambda_path(1)*ones(numexpandedsamples,1);LogLambda_path;LogLambda_path(end)*ones(numexpandedsamples,1)];
%         % Generate relative distance matrix:
%         RelDistMatrix=repmat(DistExpanded,[1,numlambdasubsamples+2*numexpandedsamples])-repmat(DistExpanded',[numlambdasubsamples+2*numexpandedsamples,1]);
%         % Generate smoothing matrix:
%         SmoothingMatrix=exp(-.5*((RelDistMatrix*sqrt(2*log(2))/(.5*DistWavelengths(ii))).^2));
%         SmoothingMatrix=SmoothingMatrix./repmat(sum(SmoothingMatrix,2),[1,numlambdasubsamples+2*numexpandedsamples]);
%         % Smooth the variables:
%         DataCostExpanded=SmoothingMatrix*DataCostExpanded;
%         RegCostExpanded=SmoothingMatrix*RegCostExpanded;
%         LambdaExpanded=SmoothingMatrix*LambdaExpanded;
%         % Identify the portions that get close enough to the
%         % endpoints:
%         ind1=find(abs(DataCostExpanded-LogDataCost_path(1))<=endpointdisttolerance*overalldist&abs(RegCostExpanded-LogRegCost_path(1))<=endpointdisttolerance*overalldist,1,'last');
%         ind2=find(abs(DataCostExpanded-LogDataCost_path(end))<=endpointdisttolerance*overalldist&abs(RegCostExpanded-LogRegCost_path(end))<=endpointdisttolerance*overalldist,1,'first');
%         % Recompute distance along curve:
%         ThisDist=[0;cumsum(sqrt(diff(DataCostExpanded).^2+diff(RegCostExpanded).^2))];
%         % Reference distance to start point:
%         ThisDist=ThisDist-ThisDist(ind1);
%         % Compute evenly spaced distance vector:
%         NewDist=linspace(0,ThisDist(ind2),numlambdasubsamples)';
%         % Reinterpolate onto evenly spaced vector:
%         LogDataCost_curves(:,ii)=interp1(ThisDist,DataCostExpanded,NewDist,lambdainterpmethod);
%         LogRegCost_curves(:,ii)=interp1(ThisDist,RegCostExpanded,NewDist,lambdainterpmethod);
%         LogLambda_curves(:,ii)=interp1(ThisDist,LambdaExpanded,NewDist,lambdainterpmethod);
%         % Recompute distance one last time:
%         Dist_curves(:,ii)=[0;cumsum(sqrt(diff(LogDataCost_curves(:,ii)).^2+diff(LogRegCost_curves(:,ii)).^2))];
%         %                 % Extract the portions we need:
%         %                 LogDataCost_curves(:,ii)=DataCostExpanded(numexpandedsamples+1:end-numexpandedsamples);
%         %                 LogRegCost_curves(:,ii)=RegCostExpanded(numexpandedsamples+1:end-numexpandedsamples);
%         %                 LogLambda_curves(:,ii)=LambdaExpanded(numexpandedsamples+1:end-numexpandedsamples);
%         %                 % Recompute distance along curve:
%         %                 Dist_curves(:,ii)=[0;cumsum(sqrt(diff(LogDataCost_curves(:,ii)).^2+diff(LogRegCost_curves(:,ii)).^2))];
%         %                 % Reinterpolate to roughly even distance spacing along curve:
%         %                 NewDist=linspace(0,max(Dist_curves(:,ii)),numdistsamples)';
%         %                 LogDataCost_curves(:,ii)=interp1(Dist_curves(:,ii),LogDataCost_curves(:,ii),NewDist,distinterpmethod);
%         %                 LogRegCost_curves(:,ii)=interp1(Dist_curves(:,ii),LogRegCost_curves(:,ii),NewDist,distinterpmethod);
%         %                 LogLambda_curves(:,ii)=interp1(Dist_curves(:,ii),LogLambda_curves(:,ii),NewDist,distinterpmethod);
%         %                 Dist_curves(:,ii)=[0;cumsum(sqrt(diff(LogDataCost_curves(:,ii)).^2+diff(LogRegCost_curves(:,ii)).^2))];
%     else
%         % No smoothing required:
%         LogDataCost_curves(:,ii)=LogDataCost_path;
%         LogRegCost_curves(:,ii)=LogRegCost_path;
%         LogLambda_curves(:,ii)=LogLambda_path;
%         Dist_curves(:,ii)=Dist_path;
%     end
%     % Compute curvature variance:
%     ThisDDist=.5*(Dist_curves(3:end,ii)-Dist_curves(1:end-2,ii));
%     ThisLogDataCurvature=(LogDataCost_curves(1:end-2,ii)+LogDataCost_curves(3:end,ii)-2*LogDataCost_curves(2:end-1,ii))./(ThisDDist.^2);
%     ThisLogRegCurvature=(LogRegCost_curves(1:end-2,ii)+LogRegCost_curves(3:end,ii)-2*LogRegCost_curves(2:end-1,ii))./(ThisDDist.^2);
%     CurvatureRMS_curves(ii)=sqrt(mean(ThisLogDataCurvature.^2+ThisLogRegCurvature.^2));
%     % Compute scatter variance:
%     TheseDist=sqrt((repmat(LogDataCost_points,[1,numlambdasubsamples])-repmat(LogDataCost_curves(:,ii)',[numlambda-sum(IsOutlier),1])).^2+...
%         (repmat(LogRegCost_points,[1,numlambdasubsamples])-repmat(LogRegCost_curves(:,ii)',[numlambda-sum(IsOutlier),1])).^2);
%     TheseMinDist=min(TheseDist,[],2);
%     ScatterRMS_curves(ii)=sqrt(mean(TheseMinDist.^2));
% end
%
% % Determine best wavelength:
% TotalRMS_curves=ScatterRMS_curves+CurvatureRMS_curves;
% bestwavelengthind=find(TotalRMS_curves==min(TotalRMS_curves),1,'first');
% bestwavelength=DistWavelengths(bestwavelengthind);
%
% % Isolate best curve:
% LogDataCost_bestcurve=LogDataCost_curves(:,bestwavelengthind);
% LogRegCost_bestcurve=LogRegCost_curves(:,bestwavelengthind);
% LogLambda_bestcurve=LogLambda_curves(:,bestwavelengthind);
% Dist_bestcurve=Dist_curves(:,bestwavelengthind);
%
% % Compute curvature of best curve:
% ddist_bestcurve=Dist_bestcurve(2)-Dist_bestcurve(1);
% LogDataCurvature_bestcurve=(LogDataCost_bestcurve(1:end-2)+LogDataCost_bestcurve(3:end)-2*LogDataCost_bestcurve(2:end-1))/(ddist_bestcurve^2);
% LogRegCurvature_bestcurve=(LogRegCost_bestcurve(1:end-2)+LogRegCost_bestcurve(3:end)-2*LogRegCost_bestcurve(2:end-1))/(ddist_bestcurve^2);
% TotalCurvature_bestcurve=LogDataCurvature_bestcurve+LogRegCurvature_bestcurve;
%
% % Guess lambda uncertainty:
% DLogLambdaDDist_bestcurve=diff(LogLambda_bestcurve)./diff(Dist_bestcurve);
% lambdauncertaintyguess=.5*bestwavelength*abs(mean(DLogLambdaDDist_bestcurve));
%
% % Identify best lambda:
% meandloglambda=mean(diff(LogLambda_bestcurve));
% bestlambdaind=1+find(TotalCurvature_bestcurve==max(TotalCurvature_bestcurve(ceil(.5*lambdauncertaintyguess/meandloglambda):floor(numlambdasubsamples-1-.5*lambdauncertaintyguess/meandloglambda))),1,'first');
% bestdist=Dist_bestcurve(bestlambdaind);
% bestlambda=exp(LogLambda_bestcurve(bestlambdaind));
% bestmodelind=find(abs(log(Lambda)-log(bestlambda))==min(abs(log(Lambda)-log(bestlambda))),1,'first');
%
% % Compute uncertainty in best lambda:
% bestlambda_loguncertainty=.5*bestwavelength*abs(interp1(.5*(Dist_bestcurve(1:end-1)+Dist_bestcurve(2:end)),DLogLambdaDDist_bestcurve,bestdist,lambdainterpmethod));
%
% % Compute max/min allowable lambdas:
% bestlambda_minallowable=exp(log(bestlambda)-bestlambda_loguncertainty);
% bestlambda_maxallowable=exp(log(bestlambda)+bestlambda_loguncertainty);
% bestlambdaind_minallowable=find(LogLambda_bestcurve>=log(bestlambda_minallowable),1,'first');
% bestlambdaind_maxallowable=find(LogLambda_bestcurve<=log(bestlambda_maxallowable),1,'last');


% %% Make L-Curve Figure for Lambda Selection via 2D Curve-Fitting:
%
% % Compute boxes for L-Curve figure:
% if strcmp(lcurve_orientation,'horz')
%     Boxes={[horzaxesbuffer_lcurve,vertaxesbuffer_lcurve,lcurvefrac-horzaxesbuffer_lcurve-horznotextbuffer_lcurve,1-titlebuffer_lcurve-vertaxesbuffer_lcurve],...
%         [lcurvefrac+horzaxesbuffer_lcurve,.5+vertaxesbuffer_lcurve,1-lcurvefrac-horzaxesbuffer_lcurve-horznotextbuffer_lcurve,.5-titlebuffer_lcurve-vertaxesbuffer_lcurve],...
%         [lcurvefrac+horzaxesbuffer_lcurve,vertaxesbuffer_lcurve,1-lcurvefrac-horzaxesbuffer_lcurve-horznotextbuffer_lcurve,.5-titlebuffer_lcurve-vertaxesbuffer_lcurve]};
% else
%     Boxes={[horzaxesbuffer_lcurve,1-lcurvefrac+vertaxesbuffer_lcurve,1-horzaxesbuffer_lcurve-horznotextbuffer_lcurve,lcurvefrac-titlebuffer_lcurve-vertaxesbuffer_lcurve],...
%         [horzaxesbuffer_lcurve,vertaxesbuffer_lcurve,.5-horzaxesbuffer_lcurve-horznotextbuffer_lcurve,1-lcurvefrac-titlebuffer_lcurve-vertaxesbuffer_lcurve],...
%         [.5+horzaxesbuffer_lcurve,vertaxesbuffer_lcurve,.5-horzaxesbuffer_lcurve-horznotextbuffer_lcurve,1-lcurvefrac-titlebuffer_lcurve-vertaxesbuffer_lcurve]};
% end
%
% % Make figure:
% figure
%
% % Plot cost vs cost, log scale:
% % Make subplot:
% subplot('Position',Boxes{1})
% % Plot model costs:
% plot(RegularizationCost(~IsOutlier),DataCost(~IsOutlier),'Marker','o','MarkerFaceColor','k','Color','k','LineStyle','none')
% hold on
% if sum(IsOutlier)>0
%     plot(RegularizationCost(IsOutlier),DataCost(IsOutlier),'Marker','o','MarkerFaceColor','m','Color','m','LineStyle','none')
% end
% % Plot best curve:
% plot(exp(LogRegCost_bestcurve),exp(LogDataCost_bestcurve),'k','LineWidth',1)
% set(gca,'XScale','log')
% set(gca,'YScale','log')
% % Add lambda labels:
% %         for ii=1:3:numlambda
% %             text(RegularizationCost(ii),DataCost(ii),num2str(Lambda(ii)),'Color','k','HorizontalAlignment','right','VerticalAlignment','top')
% %         end
% % Highlight best value:
% plot(exp(LogRegCost_bestcurve(bestlambdaind)),exp(LogDataCost_bestcurve(bestlambdaind)),'Color','r','MarkerFaceColor','r','Marker','d')
% lastpow10=10^floor(log10(bestlambda));
% text(exp(LogRegCost_bestcurve(bestlambdaind)+.25),exp(LogDataCost_bestcurve(bestlambdaind)+.25),{'Best Lambda=';num2str(.1*lastpow10*round(10*bestlambda/lastpow10))},'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k')
% % Highlight allowable range:
% plot(exp(LogRegCost_bestcurve(bestlambdaind_minallowable:bestlambdaind_maxallowable)),...
%     exp(LogDataCost_bestcurve(bestlambdaind_minallowable:bestlambdaind_maxallowable)),'Color','r','LineWidth',2)
% angle1=atan2(LogDataCost_bestcurve(bestlambdaind_minallowable+1)-LogDataCost_bestcurve(bestlambdaind_minallowable),...
%     LogRegCost_bestcurve(bestlambdaind_minallowable+1)-LogRegCost_bestcurve(bestlambdaind_minallowable));
% angle2=atan2(LogDataCost_bestcurve(bestlambdaind_maxallowable)-LogDataCost_bestcurve(bestlambdaind_maxallowable-1),...
%     LogRegCost_bestcurve(bestlambdaind_maxallowable)-LogRegCost_bestcurve(bestlambdaind_maxallowable-1));
% plot(exp(LogRegCost_bestcurve(bestlambdaind_minallowable)+.2*cos(angle1+pi/2)*[-1,1]),...
%     exp(LogDataCost_bestcurve(bestlambdaind_minallowable)+.2*sin(angle1+pi/2)*[-1,1]),'Color','r','LineWidth',2)
% plot(exp(LogRegCost_bestcurve(bestlambdaind_maxallowable)+.2*cos(angle2+pi/2)*[-1,1]),...
%     exp(LogDataCost_bestcurve(bestlambdaind_maxallowable)+.2*sin(angle2+pi/2)*[-1,1]),'Color','r','LineWidth',2)
% lastpow10_1=10^floor(log10(bestlambda_minallowable));
% lastpow10_2=10^floor(log10(bestlambda_maxallowable));
% if bestlambdaind_maxallowable<numlambdasubsamples
%     text(exp(LogRegCost_bestcurve(bestlambdaind)-.25),exp(LogDataCost_bestcurve(bestlambdaind)-.25),...
%         {'Acceptable Range=';[num2str(.1*lastpow10_1*round(10*bestlambda_minallowable/lastpow10_1)),'-',num2str(.1*lastpow10_2*round(10*bestlambda_maxallowable/lastpow10_2))]},...
%         'HorizontalAlignment','right','VerticalAlignment','top','Color','k')
% else
%     text(exp(LogRegCost_bestcurve(bestlambdaind_minallowable)-.25),exp(LogDataCost_bestcurve(bestlambdaind_minallowable)-.25),...
%         {'Acceptable Range=';[num2str(.1*lastpow10_1*round(10*bestlambda_minallowable/lastpow10_1)),'-',num2str(.1*lastpow10_2*round(10*bestlambda_maxallowable/lastpow10_2))]},...
%         'HorizontalAlignment','right','VerticalAlignment','top','Color','k')
% end
% % Compute axes limits:
% thesexlims=10.^[floor(log10(min(RegularizationCost(~IsOutlier)))),ceil(log10(max(RegularizationCost(~IsOutlier))))];
% theseylims=10.^[floor(log10(min([DataCost(~IsOutlier);exp(LogDataCost_bestcurve(bestlambdaind)-.5)]))),ceil(log10(max(DataCost(~IsOutlier))))];
% % Indicate smoothing wavelength:
% thiscentx=exp(log(thesexlims(2))-.5*bestwavelength-.25);
% thiscenty=exp(log(theseylims(1))+.55*[log(theseylims(2))-log(theseylims(1))]);
% plot(thiscentx*[exp(-.5*bestwavelength),exp(.5*bestwavelength)],thiscenty*[1,1],'k')
% plot(thiscentx*exp(-.5*bestwavelength)*[1,1],thiscenty*[.9,1.1],'k')
% plot(thiscentx*exp(.5*bestwavelength)*[1,1],thiscenty*[.9,1.1],'k')
% text(thiscentx,thiscenty*1.1,{'Smoothing';'wavelength'},'HorizontalAlignment','center','VerticalAlignment','bottom')
% % Set up axes etc:
% xlim(thesexlims)
% ylim(theseylims)
% set(gca,'PlotBoxAspectRatio',[(log(thesexlims(2))-log(thesexlims(1)))/(log(theseylims(2))-log(theseylims(1))),1,1])
% xlabel('Regularization Cost (unitless)')
% ylabel('Data Cost (unitless)')
% if sum(IsOutlier)>0
%     legend('Raw model results (good points)','Raw model results (outliers)','Smoothed tradeoff curve','location','NorthEast')
% else
%     legend('Raw model results','Smoothed tradeoff curve','location','NorthEast')
% end
% title('a) L-Curve Plot (cost vs cost)')
%
% % Plot cost curvature against lambda values:
% % Make subplot:
% subplot('Position',Boxes{2})
% % Plot best curve curvatures:
% plot(exp(LogLambda_bestcurve(2:end-1)),LogDataCurvature_bestcurve,'b')
% hold on
% plot(exp(LogLambda_bestcurve(2:end-1)),LogRegCurvature_bestcurve,'r')
% plot(exp(LogLambda_bestcurve(2:end-1)),TotalCurvature_bestcurve,'k')
% % Show best value and uncertainty:
% plot(exp(LogLambda_bestcurve(bestlambdaind)),TotalCurvature_bestcurve(bestlambdaind-1),'Color','k','MarkerFaceColor','k','Marker','o')
% plot([bestlambda_minallowable,bestlambda_maxallowable],TotalCurvature_bestcurve(bestlambdaind-1)*[1,1],'k')
% plot(bestlambda_minallowable*[1,1],TotalCurvature_bestcurve(bestlambdaind-1)*[.975,1.025],'k')
% plot(bestlambda_maxallowable*[1,1],TotalCurvature_bestcurve(bestlambdaind-1)*[.975,1.025],'k')
% text(bestlambda,TotalCurvature_bestcurve(bestlambdaind-1)*1.025,{'Best lambda';'and uncertainty'},'HorizontalAlignment','center','VerticalAlignment','bottom')
% % Determine axes limits:
% thesexlims=10.^[floor(log10(min(Lambda))),ceil(log10(max(Lambda)))];
% theseylims=TotalCurvature_bestcurve(bestlambdaind-1)*[-.25,1.25];
% % Plot zero line:
% plot(thesexlims,[0,0],'--k')
% % Set up axes etc:
% xlim(thesexlims)
% ylim(theseylims)
% set(gca,'XScale','log')
% set(gca,'XTick',10.^(log10(thesexlims(1)):log10(thesexlims(2))))
% xlabel('Lambda Value (unitless)')
% ylabel('Curvature (d^2(log(cost))/ds^2)')
% title({'b) Cost Curvature vs Lambda';'(selection of optimal tradeoff)'})
% if log(bestlambda)>mean(log([lambdamin,lambdamax]))
%     locationstring='NorthWest';
% else
%     locationstring='NorthEast';
% end
% legend('Data','Regularization','Total curvature','location',locationstring)
%
% % Plot variance against wavelength:
% % Make subplot:
% subplot('Position',Boxes{3})
% % Plot RMS values:
% plot(DistWavelengths/log(10),ScatterRMS_curves,'b')
% hold on
% plot(DistWavelengths/log(10),CurvatureRMS_curves,'r')
% plot(DistWavelengths/log(10),TotalRMS_curves,'k')
% % Plot and label best value:
% plot(DistWavelengths(bestwavelengthind)/log(10),TotalRMS_curves(bestwavelengthind),'Color','k','MarkerFaceColor','k','Marker','o')
% lastpow10=min(1,10^floor(log10(bestwavelength/log(10))));
% text(DistWavelengths(bestwavelengthind)/log(10),TotalRMS_curves(bestwavelengthind),{'Best Wavelength=';num2str(.01*lastpow10*round(100*bestwavelength/(log(10)*lastpow10)))},'HorizontalAlignment','center','VerticalAlignment','bottom','Color','k')
% % Determine axes limits:
% thesexlims=[0,maxlambdawavelength/log(10)];
% theseylims=max(TotalRMS_curves)*[0,1.25];
% % Set up axes etc:
% xlim(thesexlims)
% ylim(theseylims)
% xlabel('Wavelength (distance in log/log space)')
% ylabel('RMS variability (unitless)')
% title({'c) Curve Variance vs Wavelength';'(selection of optimal smoothing)'})
% legend('Scatter','Curvature','Total variance','location','NorthEast')
%
% % Save figure:
% thisfigname=[outputfile(1:end-4),'_LCurvePlots.png'];
% set(gcf,'PaperSize',lcurvepagesize)
% set(gcf,'PaperPosition',[0,0,lcurvepagesize])
% print('-dpng',thisfigname,['-r',num2str(resolution)])



% if isfield(md.mesh,'numberofelements2d')==0
%     patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements,'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(md.results.StressbalanceSolution.Vel))
% else
%     patch('Vertices',[md.mesh.x2d/1000,md.mesh.y2d/1000],'Faces',md.mesh.elements2d,'EdgeColor','none','FaceColor','interp','FaceVertexCData',log10(md.results.StressbalanceSolution.Vel(md.mesh.vertexonsurface==1)))
% end
% caxis([0,3.5])
% xlims=get(gca,'Xlim');
% ylims=get(gca,'Ylim');
% set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
% xlabel('Easting (km)')
% ylabel('Northing (km)')
% hc=colorbar;
% set(get(hc,'Ylabel'),'String','log_{10}(m/yr)')
% title(['Velocity solution for test forward model with flow equation=',flowequation,' and sliding m=',num2str(m)])
% thisfigname=[outputfile(1:end-4),'_velmap.png'];
% set(gcf,'PaperSize',pagesize1)
% set(gcf,'PaperPosition',[0,0,pagesize1])
% print('-dpng',thisfigname,['-r',num2str(resolution)])

% % Commands for column average velocity and sliding fraction:
% Zhat=(md.mesh.z(1:numvertices2d:end)-md.mesh.z(1))/md.geometry.thickness(1);
% DZhat=[.5*Zhat(2);.5*(Zhat(3:end)-Zhat(1:end-2));.5*(1-Zhat(end-1))];
% Ubar=zeros(numvertices2d,1);
% for ii=1:numvertices2d
%     Ubar(ii)=sum(DZhat.*md.results.StressbalanceSolution.Vel(ii:numvertices2d:end));
% end
% SlidingVelocity=md.results.StressbalanceSolution.Vel(1:numvertices2d);
% SurfaceVelocity=md.results.StressbalanceSolution.Vel((numlayers-1)*numvertices2d+1:end);


% % Plot costs:
% if includedrag && includerheology
%     error('Sorry, I dont know how to plot this yet.')
% elseif includedrag
%     % Plot:
%     plot(DragWavelengths/1000,DataCost,'Marker','d','MarkerFaceColor','b','Color','b')
%     hold on
%     plot(DragWavelengths/1000,RegularizationCost,'Marker','d','MarkerFaceColor','r','Color','r')
%     plot(DragWavelengths/1000,DataCost.*RegularizationCost,'Marker','d','MarkerFaceColor','k','Color','k')
%     % Determine x limits:
%     thesexlims=10.^[floor(log10(min(DragWavelengths/1000))),ceil(log10(max(DragWavelengths/1000)))];
% else
%     % Plot:
%     plot(BbarWavelengths/1000,DataCost,'Marker','d','MarkerFaceColor','b','Color','b')
%     hold on
%     plot(BbarWavelengths/1000,RegularizationCost,'Marker','d','MarkerFaceColor','r','Color','r')
%     plot(BbarWavelengths/1000,DataCost.*RegularizationCost,'Marker','d','MarkerFaceColor','k','Color','k')
%     % Determine x limits:
%     thesexlims=10.^[floor(log10(min(BbarWavelengths/1000))),ceil(log10(max(BbarWavelengths/1000)))];
% end

% % Plot cost gradients:
% if includedrag && includerheology
%     error('Sorry, I dont know how to plot this yet.')
% elseif includedrag
%     % Compute cost gradient (on log-log scale):
%     DataCostGradient=diff(log(DataCost))./diff(log(DragWavelengths));
%     RegCostGradient=diff(log(RegularizationCost))./diff(log(DragWavelengths));
%     WavelengthMidpoints=sqrt(DragWavelengths(1:end-1).*DragWavelengths(2:end));
%     CombinedCostGradient=sqrt(DataCostGradient.^2+RegCostGradient.^2);
%     % Plot:
%     %             plot(WavelengthMidpoints/1000,DataCostGradient,'Marker','d','MarkerFaceColor','b','Color','b')
%     %             hold on
%     %             plot(WavelengthMidpoints/1000,-RegCostGradient,'Marker','d','MarkerFaceColor','r','Color','r')
%     %             plot(WavelengthMidpoints/1000,CombinedCostGradient,'Marker','d','MarkerFaceColor','k','Color','k')
%     %             histogram('BinEdges',DragWavelengths'/1000,'BinCounts',max(0,DataCostGradient),'DisplayStyle','stairs','EdgeColor','b')
%     %             hold on
%     %             histogram('BinEdges',DragWavelengths'/1000,'BinCounts',max(0,-RegCostGradient),'DisplayStyle','stairs','EdgeColor','r')
%     %             histogram('BinEdges',DragWavelengths'/1000,'BinCounts',max(0,CombinedCostGradient),'DisplayStyle','stairs','EdgeColor','k')
%     plot([DragWavelengths(1);reshape(repmat(DragWavelengths(2:end-1)',[2,1]),[],1);DragWavelengths(end)]/1000,reshape(repmat(DataCostGradient',[2,1]),[],1),'Color','b')
%     hold on
%     plot([DragWavelengths(1);reshape(repmat(DragWavelengths(2:end-1)',[2,1]),[],1);DragWavelengths(end)]/1000,reshape(repmat(-RegCostGradient',[2,1]),[],1),'Color','r')
%     plot([DragWavelengths(1);reshape(repmat(DragWavelengths(2:end-1)',[2,1]),[],1);DragWavelengths(end)]/1000,reshape(repmat(CombinedCostGradient',[2,1]),[],1),'Color','k')
%     % Determine x limits:
%     thesexlims=10.^[floor(log10(min(DragWavelengths/1000))),ceil(log10(max(DragWavelengths/1000)))];
%     % Identify best midpoint:
%     bestmidind=find(CombinedCostGradient==max(CombinedCostGradient),1,'first');
%     bestmidwavelength=WavelengthMidpoints(bestmidind);
%     % Identify best wavelength and lambda value:
%     if CombinedCostGradient(bestmidind+1)>CombinedCostGradient(bestmidind-1)
%         % Larger value is better:
%         bestind=bestmidind+1;
%     else
%         % Smaller value is better:
%         bestind=bestmidind;
%     end
%     % Record best lambda and wavelength:
%     bestwavelength=DragWavelengths(bestind);
%     bestlambda=Lambda(bestind);
%     % Label best wavelength:
%     plot(DragWavelengths(bestind)/1000,CombinedCostGradient(bestmidind),'Marker','d','MarkerFaceColor','k','COlor','k')
%     text(DragWavelengths(bestind)/1000,CombinedCostGradient(bestmidind),...
%         {'Best Wavelength=';[num2str(.01*round(100*bestwavelength/1000)),' km']},...
%         'HorizontalAlignment','center','VerticalAlignment','bottom')
% else
%     % Compute cost gradient (on log-log scale):
%     DataCostGradient=diff(log(DataCost))./diff(log(BbarWavelengths));
%     RegCostGradient=diff(log(RegularizationCost))./diff(log(BbarWavelengths));
%     WavelengthMidpoints=sqrt(BbarWavelengths(1:end-1).*BbarWavelengths(2:end));
%     CombinedCostGradient=sqrt(DataCostGradient.^2+RegCostGradient.^2);
%     % Plot:
%     %             plot(WavelengthMidpoints/1000,DataCostGradient,'Marker','d','MarkerFaceColor','b','Color','b')
%     %             hold on
%     %             plot(WavelengthMidpoints/1000,-RegCostGradient,'Marker','d','MarkerFaceColor','r','Color','r')
%     %             plot(WavelengthMidpoints/1000,CombinedCostGradient,'Marker','d','MarkerFaceColor','k','Color','k')
%     %             histogram('BinEdges',BbarWavelengths'/1000,'BinCounts',max(0,DataCostGradient),'DisplayStyle','stairs','EdgeColor','b')
%     %             hold on
%     %             histogram('BinEdges',BbarWavelengths'/1000,'BinCounts',max(0,-RegCostGradient),'DisplayStyle','stairs','EdgeColor','r')
%     %             histogram('BinEdges',BbarWavelengths'/1000,'BinCounts',max(0,CombinedCostGradient),'DisplayStyle','stairs','EdgeColor','k')
%     plot([BbarWavelengths(1);reshape(repmat(BbarWavelengths(2:end-1)',[2,1]),[],1);BbarWavelengths(end)]/1000,reshape(repmat(DataCostGradient',[2,1]),[],1),'Color','b')
%     hold on
%     plot([BbarWavelengths(1);reshape(repmat(BbarWavelengths(2:end-1)',[2,1]),[],1);BbarWavelengths(end)]/1000,reshape(repmat(-RegCostGradient',[2,1]),[],1),'Color','r')
%     plot([BbarWavelengths(1);reshape(repmat(BbarWavelengths(2:end-1)',[2,1]),[],1);BbarWavelengths(end)]/1000,reshape(repmat(CombinedCostGradient',[2,1]),[],1),'Color','k')
%     % Determine x limits:
%     thesexlims=10.^[floor(log10(min(BbarWavelengths/1000))),ceil(log10(max(BbarWavelengths/1000)))];
%     % Identify best midpoint:
%     bestmidind=find(CombinedCostGradient==max(CombinedCostGradient),1,'first');
%     bestmidwavelength=WavelengthMidpoints(bestmidind);
%     % Identify best wavelength and lambda value:
%     if CombinedCostGradient(bestmidind+1)>CombinedCostGradient(bestmidind-1)
%         % Larger value is better:
%         bestind=bestmidind+1;
%     else
%         % Smaller value is better:
%         bestind=bestmidind;
%     end
%     % Record best lambda and wavelength:
%     bestwavelength=BbarWavelengths(bestind);
%     bestlambda=Lambda(bestind);
%     % Label best wavelength:
%     plot(BbarWavelengths(bestind)/1000,CombinedCostGradient(bestmidind),'Marker','d','MarkerFaceColor','k','COlor','k')
%     text(BbarWavelengths(bestind)/1000,CombinedCostGradient(bestmidind),...
%         {'Best Wavelength=';[num2str(.01*round(100*bestwavelength/1000)),' km']},...
%         'HorizontalAlignment','center','VerticalAlignment','bottom')
% end

% % Plot cost gradients:
% plot([Lambda(1);reshape(repmat(Lambda(2:end-1)',[2,1]),[],1);Lambda(end)],reshape(repmat(DataCostGradient',[2,1]),[],1),'Color','b')
% hold on
% plot([Lambda(1);reshape(repmat(Lambda(2:end-1)',[2,1]),[],1);Lambda(end)],reshape(repmat(-RegCostGradient',[2,1]),[],1),'Color','r')
% plot([Lambda(1);reshape(repmat(Lambda(2:end-1)',[2,1]),[],1);Lambda(end)],reshape(repmat(CombinedCostGradient',[2,1]),[],1),'Color','k')


% % Pre-allocate gradients:
% SurfGradX_elements=zeros(numelements2d,1);
% SurfGradY_elements=zeros(numelements2d,1);
% % Loop through elements to compute surface gradient:
% for jj=1:numelements2d
%     % Construct matrix:
%     ThisMatrix=[x2d(elements2d(jj,:)),y2d(elements2d(jj,:)),ones(3,1)];
%     % Solve matrix:
%     thissolution=ThisMatrix\[md.geometry.surface(elements2d(jj,:))];
%     % Assign gradient values:
%     SurfGradX_elements(jj)=thissolution(1);
%     SurfGradY_elements(jj)=thissolution(2);
% end


% % Compute scatter estimates:
% % These are RMS scatter as a function of lambda
% % Compute squared deviations:
% DataDeviation=(log(DataCost)-log(DataCostSmooth)).^2;
% RegDeviation=(log(RegularizationCost)-log(RegCostSmooth)).^2;
% % COmpute weighting matrix:
% Weight=exp(-.5*(((repmat(log(Lambda),[1,numlambda])-repmat(log(Lambda'),[numlambda,1]))/log(lambdasmoothing)).^2));
% Weight=Weight./repmat(sum(Weight,2),[1,numlambda]);
% % Compute RMS scatter:
% DataScatter=sqrt(Weight*DataDeviation);
% RegScatter=sqrt(Weight*RegDeviation);
% TotalScatter=sqrt(DataScatter.^2+RegScatter.^2);


% numlambdasubsamples=10*numlambda;
%
% % Identify best lambda smoothing value:
% % Generate wavelength sampling:
% LambdaSmoothings=exp(linspace(log(minlambdasmoothing),log(maxlambdasmoothing),numlambdasmoothing));
% % Generate lambda subsampling:
% LambdaSubSample=exp(linspace(log(lambdamin),log(lambdamax),numlambdasubsamples)');
% % Interpolate costs to the subsampling:
% DataCostSubSample=exp(interp1(log(Lambda),log(DataCost),log(LambdaSubSample),'pchip'));
% RegCostSubSample=exp(interp1(log(Lambda),log(RegularizationCost),log(LambdaSubSample),'pchip'));
% % Pre-allocate:
% LogDataCost_curves=zeros(numlambdasubsamples,numlambdasmoothing);
% LogRegCost_curves=zeros(numlambdasubsamples,numlambdasmoothing);
% DataCostGradient=zeros(numlambdasubsamples-1,numlambdasmoothing);
% RegCostGradient=zeros(numlambdasubsamples-1,numlambdasmoothing);
% LocalCurvature=zeros(numlambdasubsamples-2,numlambdasmoothing);
% % Compute dloglambda:
% dloglambda=log(LambdaSubSample(2))-log(LambdaSubSample(1));  % Assumes lambda evenly spaced on log scale
% % Loop through potential smoothings:
% for smoothind=1:numlambdasmoothing
%     % Smooth cost vectors on log scale:
%     if LambdaSmoothings(smoothind)~=1
%         LogDataCost_curves(:,smoothind)=exp(intuitive_lowpass(log(DataCostSubSample),log(LambdaSmoothings(smoothind))/(sqrt(2)*dloglambda),lambdasmoothbcstyle));
%         LogRegCost_curves(:,smoothind)=exp(intuitive_lowpass(log(RegCostSubSample),log(LambdaSmoothings(smoothind))/(sqrt(2)*dloglambda),lambdasmoothbcstyle));
%     else
%         LogDataCost_curves(:,smoothind)=DataCostSubSample;
%         LogRegCost_curves(:,smoothind)=RegCostSubSample;
%     end
%     % Compute cost gradients (on log-log scale):
%     DataCostGradient(:,smoothind)=diff(log(LogDataCost_curves(:,smoothind)))./dloglambda;
%     RegCostGradient(:,smoothind)=diff(log(LogRegCost_curves(:,smoothind)))./dloglambda;
%     % Compute local slope and trig terms:
%     DeltaData=log(LogDataCost_curves(3:end,smoothind))-log(LogDataCost_curves(1:end-2,smoothind));
%     DeltaReg=log(LogRegCost_curves(3:end,smoothind))-log(LogRegCost_curves(1:end-2,smoothind));
%     DeltaS=sqrt(DeltaData.^2+DeltaReg.^2);
%     SinTheta=DeltaData./DeltaS;  % data=Y
%     CosTheta=DeltaReg./DeltaS;   % reg=X
%     % Compute local curvature:
%     % loop through lambda samples:
%     for ii=2:numlambdasubsamples-1
%         % Rotate into local coordinate system:
%         TheseXPrime=CosTheta(ii-1)*log(LogRegCost_curves(ii-1:ii+1,smoothind))+SinTheta(ii-1)*log(LogDataCost_curves(ii-1:ii+1,smoothind));
%         TheseYPrime=CosTheta(ii-1)*log(LogDataCost_curves(ii-1:ii+1,smoothind))-SinTheta(ii-1)*log(LogRegCost_curves(ii-1:ii+1,smoothind));
%         % Compute local gradient:
%         DYPrimeDXPrime=(TheseYPrime(2:end)-TheseYPrime(1:end-1))./(TheseXPrime(2:end)-TheseXPrime(1:end-1));
%         % Compute xprime midpoints:
%         TheseMidXPrime=.5*(TheseXPrime(1:end-1)+TheseXPrime(2:end));
%         % Compute local curvature:
%         LocalCurvature(ii-1,smoothind)=-(DYPrimeDXPrime(2)-DYPrimeDXPrime(1))/(TheseMidXPrime(2)-TheseMidXPrime(1));
%     end
% end
% % Compute curvature of curvature:
% CurveOfCurvature=(LocalCurvature(1:end-2,:)+LocalCurvature(3:end,:)-2*LocalCurvature(2:end-1,:))/(dloglambda^2);
% % Compute RMS of curvature of curvature:
% RMSCurveOfCurvature=sqrt(mean(CurveOfCurvature.^2,1));
% % Compute curvature of RMS:
% dloglambdasmoothing=log(LambdaSmoothings(2))-log(LambdaSmoothings(1));
% CurveOfRMSCurveOfCurvature=(RMSCurveOfCurvature(1:end-2)+RMSCurveOfCurvature(3:end)-2*RMSCurveOfCurvature(2:end-1))/(dloglambdasmoothing^2);
% % Identify best lambda smoothing:
% bestsmoothind=1+find(CurveOfRMSCurveOfCurvature==max(CurveOfRMSCurveOfCurvature),1,'first');
% bestlambdasmoothing=LambdaSmoothings(bestsmoothind);
% % Identify best lambda:
% bestlambdaind=1+find(LocalCurvature(:,bestsmoothind)==max(LocalCurvature(:,bestsmoothind)),1,'first');
% bestlambda=LambdaSubSample(bestlambdaind);
% error('end here')
% % Identify best lambda smoothing value:
% % Generate lambda smoothing sampling:
% LambdaSmoothings=exp(linspace(log(minlambdasmoothing),log(maxlambdasmoothing),numlambdasmoothing));
% % Pre-allocate:
% LogDataCost_curves=zeros(numlambda,numlambdasmoothing);
% LogRegCost_curves=zeros(numlambda,numlambdasmoothing);
% DataCostGradient=zeros(numlambda-1,numlambdasmoothing);
% RegCostGradient=zeros(numlambda-1,numlambdasmoothing);
% LocalCurvature=zeros(numlambda-2,numlambdasmoothing);
% % Compute dloglambda:
% dloglambda=log(Lambda(2))-log(Lambda(1));  % Assumes lambda evenly spaced on log scale
% % Loop through potential smoothings:
% for smoothind=1:numlambdasmoothing
%     % Smooth cost vectors on log scale:
%     if LambdaSmoothings(smoothind)~=1
%         LogDataCost_curves(:,smoothind)=exp(intuitive_lowpass(log(DataCost),log(LambdaSmoothings(smoothind))/(sqrt(2)*dloglambda),lambdasmoothbcstyle));
%         LogRegCost_curves(:,smoothind)=exp(intuitive_lowpass(log(RegularizationCost),log(LambdaSmoothings(smoothind))/(sqrt(2)*dloglambda),lambdasmoothbcstyle));
%     else
%         LogDataCost_curves(:,smoothind)=DataCost;
%         LogRegCost_curves(:,smoothind)=RegularizationCost;
%     end
%     % Compute cost gradients (on log-log scale):
%     DataCostGradient(:,smoothind)=diff(log(LogDataCost_curves(:,smoothind)))./dloglambda;
%     RegCostGradient(:,smoothind)=diff(log(LogRegCost_curves(:,smoothind)))./dloglambda;
%     % Compute local slope and trig terms:
%     DeltaData=log(LogDataCost_curves(3:end,smoothind))-log(LogDataCost_curves(1:end-2,smoothind));
%     DeltaReg=log(LogRegCost_curves(3:end,smoothind))-log(LogRegCost_curves(1:end-2,smoothind));
%     DeltaS=sqrt(DeltaData.^2+DeltaReg.^2);
%     SinTheta=DeltaData./DeltaS;  % data=Y
%     CosTheta=DeltaReg./DeltaS;   % reg=X
%     % Compute local curvature:
%     % loop through lambda samples:
%     for ii=2:numlambda-1
%         % Rotate into local coordinate system:
%         TheseXPrime=CosTheta(ii-1)*log(LogRegCost_curves(ii-1:ii+1,smoothind))+SinTheta(ii-1)*log(LogDataCost_curves(ii-1:ii+1,smoothind));
%         TheseYPrime=CosTheta(ii-1)*log(LogDataCost_curves(ii-1:ii+1,smoothind))-SinTheta(ii-1)*log(LogRegCost_curves(ii-1:ii+1,smoothind));
%         % Compute local gradient:
%         DYPrimeDXPrime=(TheseYPrime(2:end)-TheseYPrime(1:end-1))./(TheseXPrime(2:end)-TheseXPrime(1:end-1));
%         % Compute xprime midpoints:
%         TheseMidXPrime=.5*(TheseXPrime(1:end-1)+TheseXPrime(2:end));
%         % Compute local curvature:
%         LocalCurvature(ii-1,smoothind)=-(DYPrimeDXPrime(2)-DYPrimeDXPrime(1))/(TheseMidXPrime(2)-TheseMidXPrime(1));
%     end
% end
% % Compute curvature of curvature:
% CurveOfCurvature=(LocalCurvature(1:end-2,:)+LocalCurvature(3:end,:)-2*LocalCurvature(2:end-1,:))/(dloglambda^2);
% % Compute RMS of curvature of curvature:
% RMSCurveOfCurvature=sqrt(mean(CurveOfCurvature.^2,1));
% % Compute curvature of RMS:
% dloglambdasmoothing=log(LambdaSmoothings(2))-log(LambdaSmoothings(1));
% CurveOfRMSCurveOfCurvature=(RMSCurveOfCurvature(1:end-2)+RMSCurveOfCurvature(3:end)-2*RMSCurveOfCurvature(2:end-1))/(dloglambdasmoothing^2);
% % Identify best lambda smoothing:
% bestsmoothind=1+find(CurveOfRMSCurveOfCurvature==max(CurveOfRMSCurveOfCurvature),1,'first');
% bestlambdasmoothing=LambdaSmoothings(bestsmoothind);
% % Identify best lambda:
% bestlambdaind=1+find(LocalCurvature(:,bestsmoothind)==max(LocalCurvature(:,bestsmoothind)),1,'first');
% bestlambda=Lambda(bestlambdaind);

% %% Make L-Curve Figure:
% figure
%
% % Plot cost vs cost, log scale:
% % Make subplot:
% subplot(3,2,1)
% % Plot costs:
% plot(RegularizationCost,DataCost,'Marker','o','MarkerFaceColor','k','Color','k','LineStyle','none')
% hold on
% plot(LogRegCost_curves(:,bestsmoothind),LogDataCost_curves(:,bestsmoothind),'--k')
% set(gca,'XScale','log')
% set(gca,'YScale','log')
% % Add lambda labels:
% for ii=1:3:numlambda
%     text(RegularizationCost(ii),DataCost(ii),num2str(Lambda(ii)),'Color','k','HorizontalAlignment','right','VerticalAlignment','top')
% end
% % Highlight best value:
% hold on
% plot(RegularizationCost(bestlambdaind),DataCost(bestlambdaind),'Color','r','MarkerFaceColor','r','Marker','o')
% lastpow10=10^floor(log10(bestlambda));
% text(RegularizationCost(bestlambdaind),DataCost(bestlambdaind),{'Best Lambda=';num2str(.01*lastpow10*round(100*bestlambda/lastpow10))},'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k')
% % Set up axes etc:
% thesexlims=10.^[floor(log10(min(RegularizationCost))),ceil(log10(max(RegularizationCost)))];
% theseylims=10.^[floor(log10(min(DataCost))),ceil(log10(max(DataCost)))];
% xlim(thesexlims)
% ylim(theseylims)
% set(gca,'PlotBoxAspectRatio',[(log(thesexlims(2))-log(thesexlims(1)))/(log(theseylims(2))-log(theseylims(1))),1,1])
% xlabel('Regularization Cost (unitless)')
% ylabel('Data Cost (unitless)')
% legend('Raw Results','Best Smoothed Curve','location','NorthEast')
% title('a) L-Curve Plot (cost vs cost)')
%
% % Plot costs against lambda values:
% % Make subplot:
% subplot(3,2,2)
% % Plot costs:
% plot(Lambda,DataCost,'Marker','d','MarkerFaceColor','b','Color','b','LineStyle','none')
% hold on
% plot(Lambda,RegularizationCost,'Marker','s','MarkerFaceColor','r','Color','r','LineStyle','none')
% plot(Lambda,LogDataCost_curves(:,bestsmoothind),'--b')
% plot(Lambda,LogRegCost_curves(:,bestsmoothind),'--r')
% % Highlight best value:
% set(gca,'YScale','log')
% theseylims=10.^[floor(log10(min([DataCost;RegularizationCost]))),ceil(log10(max([DataCost;RegularizationCost])))];
% plot(bestlambda*[1,1],theseylims,'--k')
% % Determine x limits:
% thesexlims=10.^[floor(log10(min(Lambda))),ceil(log10(max(Lambda)))];
% set(gca,'XScale','log')
% % Set up axes etc:
% xlim(thesexlims)
% set(gca,'XTick',10.^(log10(thesexlims(1)):log10(thesexlims(2))))
% ylim(theseylims)
% set(gca,'YTick',10.^(log10(theseylims(1)):log10(theseylims(2))))
% xlabel('Lambda Value (unitless)')
% ylabel('Cost Function (unitless)')
% title('b) Cost vs Lambda')
% %text(thesexlims(1),theseylims(1),{'Note: gradient and curvature';'computed from smoothed curves'},'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k')
% legend('Data Cost (raw)','Reg Cost (raw)','Data Cost (smooth)','Reg Cost (smooth)','location','SouthWest')
%
% % Plot cost gradient against lambda values:
% % Make subplot:
% subplot(3,2,3)
% % Get lambda midpoints:
% LambdaMid=sqrt(Lambda(1:end-1).*Lambda(2:end));
% % Plot all cost gradients:
% h31=plot(LambdaMid,DataCostGradient,'Color',[.5,.5,1],'LineWidth',.5);
% hold on
% h32=plot(LambdaMid,-RegCostGradient,'Color',[1,.5,.5],'LineWidth',.5);
% % Plot best cost gradients:
% h33=plot(LambdaMid,DataCostGradient(:,bestsmoothind),'b','LineWidth',1);
% h34=plot(LambdaMid,-RegCostGradient(:,bestsmoothind),'r','LineWidth',1);
% % Determine limits:
% thesexlims=10.^[floor(log10(min(Lambda))),ceil(log10(max(Lambda)))];
% theseylims=max([DataCostGradient(:,bestsmoothind);-RegCostGradient(:,bestsmoothind)])*[-.25,1.25];
% % Add zero line:
% plot(thesexlims,[0,0],'--k')
% % Highlight best value:
% plot(bestlambda*[1,1],theseylims,'--k')
% % Set up axes etc:
% xlim(thesexlims)
% ylim(theseylims)
% set(gca,'XScale','log')
% set(gca,'XTick',10.^(log10(thesexlims(1)):log10(thesexlims(2))))
% xlabel('Lambda Value (unitless)')
% ylabel({'Cost Function Gradient';'d(log(cost))/d(log(lambda))'})
% title('c) Cost Gradient vs Lambda')
% legend([h31(1);h32(1);h33;h34],'Data (all wavelengths)','-Reg (all wavelengths)','Data (best wavelength)','-Reg (best wavelength)','location','NorthWest')
%
% % Plot cost curvature against lambda values:
% % Make subplot:
% subplot(3,2,4)
% % Plot all cost curvatures:
% h41=plot(Lambda(2:end-1),LocalCurvature,'Color',[.5,.5,.5],'LineWidth',.5);
% hold on
% % Plot best cost curvature:
% h42=plot(Lambda(2:end-1),LocalCurvature(:,bestsmoothind),'k','LineWidth',1);
% % Label best value:
% plot(Lambda(bestlambdaind),LocalCurvature(bestlambdaind-1,bestsmoothind),'Color','r','MarkerFaceColor','r','Marker','o')
% lastpow10=10^floor(log10(bestlambda));
% text(Lambda(bestlambdaind),LocalCurvature(bestlambdaind-1,bestsmoothind),{'Best Lambda=';num2str(.01*lastpow10*round(100*bestlambda/lastpow10))},'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k')
% % Determine axes limits:
% thesexlims=10.^[floor(log10(min(Lambda))),ceil(log10(max(Lambda)))];
% theseylims=LocalCurvature(bestlambdaind-1,bestsmoothind)*[-.25,1.25];
% % Show smoothing wavelength:
% plot(bestlambda*[1/sqrt(bestlambdasmoothing),sqrt(bestlambdasmoothing)],.1*theseylims(2)*[1,1],'k')
% plot(bestlambda*(1/sqrt(bestlambdasmoothing))*[1,1],theseylims(2)*[.05,.15],'k')
% plot(bestlambda*sqrt(bestlambdasmoothing)*[1,1],theseylims(2)*[.05,.15],'k')
% text(bestlambda,theseylims(2)*.15,{'Best wavelength';'(lambda uncertainty)'},'HorizontalAlignment','center','VerticalAlignment','bottom')
% % Plot zero line:
% plot(thesexlims,[0,0],'--k')
% % Set up axes etc:
% xlim(thesexlims)
% ylim(theseylims)
% set(gca,'XScale','log')
% set(gca,'XTick',10.^(log10(thesexlims(1)):log10(thesexlims(2))))
% xlabel('Lambda Value (unitless)')
% ylabel('Local Curvature (unitless)')
% title({'d) Cost Curvature vs Lambda';'(selection of optimal lambda)'})
% legend([h41(1);h42],'All wavelengths','Best Wavelength','location','NorthEast')
%
% % Plot curvature variability against lambda smoothing:
% % Make subplot:
% subplot(3,2,5)
% % Plot RMS of curvature of curvature:
% plot(log10(LambdaSmoothings),RMSCurveOfCurvature,'k')
% % Label best value:
% hold on
% plot(log10(LambdaSmoothings(bestsmoothind)),RMSCurveOfCurvature(bestsmoothind),'Color','r','MarkerFaceColor','r','Marker','o')
% lastpow10=10^floor(log10(bestlambdasmoothing));
% text(log10(LambdaSmoothings(bestsmoothind)),RMSCurveOfCurvature(bestsmoothind),{['Best Wavelength=',num2str(.01*round(100*log10(bestlambdasmoothing))),','];...
%     ['or factor of ',num2str(.01*lastpow10*round(100*bestlambdasmoothing/lastpow10))]},'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k')
% % Determine axes limits:
% thesexlims=[floor(log10(min(LambdaSmoothings))),ceil(log10(max(LambdaSmoothings)))];
% theseylims=max(RMSCurveOfCurvature)*[0,1.25];
% % Set up axes etc:
% xlim(thesexlims)
% ylim(theseylims)
% xlabel('Smoothing Wavelength (log_{10} scale)')
% ylabel({'Curvature Variability';'(RMS of d^{2}(curve)/d\lambda^{2})'})
% title('e) Cost Curvature Variability vs Smoothing Wavelength')
%
% % Plot curve of curvature variability against lambda smoothing:
% % Make subplot:
% subplot(3,2,6)
% % Plot curvature of RMS of curvature of curvature:
% plot(log10(LambdaSmoothings(2:end-1)),CurveOfRMSCurveOfCurvature,'k')
% % Label best value:
% hold on
% plot(log10(LambdaSmoothings(bestsmoothind)),CurveOfRMSCurveOfCurvature(bestsmoothind-1),'Color','r','MarkerFaceColor','r','Marker','o')
% lastpow10=10^floor(log10(bestlambdasmoothing));
% text(log10(LambdaSmoothings(bestsmoothind)),CurveOfRMSCurveOfCurvature(bestsmoothind-1),{['Best Wavelength=',num2str(.01*round(100*log10(bestlambdasmoothing))),','];...
%     ['or factor of ',num2str(.01*lastpow10*round(100*bestlambdasmoothing/lastpow10))]},'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k')
% % Determine axes limits:
% thesexlims=[floor(log10(min(LambdaSmoothings))),ceil(log10(max(LambdaSmoothings)))];
% theseylims=CurveOfRMSCurveOfCurvature(bestsmoothind-1)*[-.25,1.25];
% % Plot zero line:
% plot(thesexlims,[0,0],'--k')
% % Set up axes etc:
% xlim(thesexlims)
% ylim(theseylims)
% xlabel('Smoothing Wavelength (log_{10} scale)')
% ylabel({'Curve of Curvature Variability','(d^{2}(RMS)/d(wavelength)^{2})'})
% title({'f) Curvature of Curvature Variability vs Smoothing Wavelength';'(selection of optimal wavelength)'})
%
% % Save figure:
% thispagesize=[14,13];
% thisfigname=[outputfile(1:end-4),'_LCurvePlots.png'];
% set(gcf,'PaperSize',thispagesize)
% set(gcf,'PaperPosition',[0,0,thispagesize])
% print('-dpng',thisfigname,['-r',num2str(resolution)])

%% Function for computing strain rates and viscosity:

function Viscosity_elements=ComputeViscosity(md)

% Hardcoded secondsperyear:
secondsperyear=60*60*24*365.25;

% Get rheology exponent:
n=md.materials.rheology_n(1);  % assumes constant n

% Extract mesh properites:
try
    x2d=md.mesh.x2d;
    y2d=md.mesh.y2d;
    elements2d=md.mesh.elements2d;
    numlayers=md.mesh.numberoflayers;
catch
    x2d=md.mesh.x;
    y2d=md.mesh.y;
    elements2d=md.mesh.elements;
    numlayers=1;
end
numelements2d=size(elements2d,1);
numvertices2d=length(x2d);

% pre-allocate horizontal gradients:
DUDX_ud=zeros(numelements2d*numlayers,1);  % these all live on top/bottom element faces
DUDY_ud=zeros(numelements2d*numlayers,1);
DVDX_ud=zeros(numelements2d*numlayers,1);
DVDY_ud=zeros(numelements2d*numlayers,1);
if numlayers>1
    DWDX_ud=zeros(numelements2d*numlayers,1);
    DWDY_ud=zeros(numelements2d*numlayers,1);
end

% pre-allocate vertical gradients:
if numlayers>1
    DUDZ_vertedges=zeros(numvertices2d*(numlayers-1),1);   % these all live on vertical element edges
    DVDZ_vertedges=zeros(numvertices2d*(numlayers-1),1);
    DWDZ_vertedges=zeros(numvertices2d*(numlayers-1),1);
end

% Loop through 2D elements to compute horizontal gradients:
for ii=1:numelements2d
    % Construct matrix for this element:
    ThisMatrix=[x2d(elements2d(ii,:)),y2d(elements2d(ii,:)),ones(3,1)];
    % Loop through vertical layers:
    for jj=1:numlayers
        % Solve matrix for u-plane:
        thissolution=ThisMatrix\[md.results.StressbalanceSolution.Vx((jj-1)*numvertices2d+elements2d(ii,:))];
        % Assign gradient values:
        DUDX_ud((jj-1)*numelements2d+ii)=thissolution(1);
        DUDY_ud((jj-1)*numelements2d+ii)=thissolution(2);
        % Solve matrix for v-plane:
        thissolution=ThisMatrix\[md.results.StressbalanceSolution.Vy((jj-1)*numvertices2d+elements2d(ii,:))];
        % Assign gradient values:
        DVDX_ud((jj-1)*numelements2d+ii)=thissolution(1);
        DVDY_ud((jj-1)*numelements2d+ii)=thissolution(2);
        if numlayers>1
            % Solve matrix for w-plane:
            thissolution=ThisMatrix\[md.results.StressbalanceSolution.Vz((jj-1)*numvertices2d+elements2d(ii,:))];
            % Assign gradient values:
            DWDX_ud((jj-1)*numelements2d+ii)=thissolution(1);
            DWDY_ud((jj-1)*numelements2d+ii)=thissolution(2);
        end
    end
end

% Loop through 2D vertices to compute vertical gradients:
if numlayers>1
    for ii=1:numvertices2d
        % Get this DZ vector:
        ThisDZ=md.mesh.z(ii+numvertices2d*(1:numlayers-1)')-md.mesh.z(ii+numvertices2d*(0:numlayers-2)');
        % Compute all three velocity gradients:
        DUDZ_vertedges(ii+numvertices2d*(0:numlayers-2)')=(md.results.StressbalanceSolution.Vx(ii+numvertices2d*(1:numlayers-1)')-md.results.StressbalanceSolution.Vx(ii+numvertices2d*(0:numlayers-2)'))./ThisDZ;
        DVDZ_vertedges(ii+numvertices2d*(0:numlayers-2)')=(md.results.StressbalanceSolution.Vy(ii+numvertices2d*(1:numlayers-1)')-md.results.StressbalanceSolution.Vy(ii+numvertices2d*(0:numlayers-2)'))./ThisDZ;
        DWDZ_vertedges(ii+numvertices2d*(0:numlayers-2)')=(md.results.StressbalanceSolution.Vz(ii+numvertices2d*(1:numlayers-1)')-md.results.StressbalanceSolution.Vz(ii+numvertices2d*(0:numlayers-2)'))./ThisDZ;
    end
end

% Average horizontal gradients from top/bottom faces onto element centers:
if numlayers>1
    DUDX_elements=.5*(DUDX_ud(1:numelements2d*(numlayers-1))+DUDX_ud(numelements2d+1:numelements2d*numlayers));
    DUDY_elements=.5*(DUDY_ud(1:numelements2d*(numlayers-1))+DUDY_ud(numelements2d+1:numelements2d*numlayers));
    DVDX_elements=.5*(DVDX_ud(1:numelements2d*(numlayers-1))+DVDX_ud(numelements2d+1:numelements2d*numlayers));
    DVDY_elements=.5*(DVDY_ud(1:numelements2d*(numlayers-1))+DVDY_ud(numelements2d+1:numelements2d*numlayers));
    DWDX_elements=.5*(DWDX_ud(1:numelements2d*(numlayers-1))+DWDX_ud(numelements2d+1:numelements2d*numlayers));
    DWDY_elements=.5*(DWDY_ud(1:numelements2d*(numlayers-1))+DWDY_ud(numelements2d+1:numelements2d*numlayers));
else
    DUDX_elements=DUDX_ud;
    DUDY_elements=DUDY_ud;
    DVDX_elements=DVDX_ud;
    DVDY_elements=DVDY_ud;
end

% Average vertical gradients from vertical edges to element centers:
if numlayers>1
    % pre-allocate:
    DUDZ_elements=zeros(numelements2d*(numlayers-1),1);
    DVDZ_elements=zeros(numelements2d*(numlayers-1),1);
    DWDZ_elements=zeros(numelements2d*(numlayers-1),1);
    % Loop through layers:
    for ii=1:numlayers-1
        % Perform averaging:
        DUDZ_elements((ii-1)*numelements2d+1:ii*numelements2d)=mean(DUDZ_vertedges((ii-1)*numvertices2d+elements2d),2);
        DVDZ_elements((ii-1)*numelements2d+1:ii*numelements2d)=mean(DVDZ_vertedges((ii-1)*numvertices2d+elements2d),2);
        DWDZ_elements((ii-1)*numelements2d+1:ii*numelements2d)=mean(DWDZ_vertedges((ii-1)*numvertices2d+elements2d),2);
    end
end

% Compute strain rate components from velocity gradients:
if numlayers>1
    StrainRateXZ_elements=.5*(DUDZ_elements+DWDX_elements);
    StrainRateYZ_elements=.5*(DVDZ_elements+DWDY_elements);
end
StrainRateXY_elements=.5*(DUDY_elements+DVDX_elements);
StrainRateXX_elements=DUDX_elements;
StrainRateYY_elements=DVDY_elements;
if numlayers>1
    StrainRateZZ_elements=DWDZ_elements;
end

% Compute effective strain rate:
if numlayers>1
    StrainRateMag_elements=sqrt(StrainRateXY_elements.^2+StrainRateXZ_elements.^2+StrainRateYZ_elements.^2+.5*(StrainRateXX_elements.^2+StrainRateYY_elements.^2+StrainRateZZ_elements.^2));
else
    StrainRateMag_elements=sqrt(StrainRateXY_elements.^2+.5*(StrainRateXX_elements.^2+StrainRateYY_elements.^2));
end

% Average rheology onto the elements:
B_elements=mean(md.materials.rheology_B(md.mesh.elements),2);

% Compute viscosity:
Viscosity_elements=B_elements.*((StrainRateMag_elements/secondsperyear).^(-(n-1)/n));

end
%
% % Compute vertically averaged viscosity:
% DZhat=(md.mesh.z(1+numvertices2d*(1:numlayers-1)')-md.mesh.z(1+numvertices2d*(0:numlayers-2)'))/md.geometry.thickness(1);
% ViscosityBar_elements2d=sum(reshape(Viscosity_elements,[numelements2d,numlayers-1]).*repmat(DZhat',[numelements2d,1]),2);
