% BestDragFigure_v1

% Lea-Sophie Hoeyns, 01/21/2024

% This script makes figures for the best basal drag map and the basal drag
% coefficient map including three zoom-in boxes for Pine Island and 
% Thwaites Glacier, as well as for the Siple Coast region. 

% Colormaps citation: 
% [1] Zhaoxu Liu / slandarer (2024). 200 colormap (https://www.mathworks.com/matlabcentral/fileexchange/120088-200-colormap), 
% MATLAB Central File Exchange. Retrieved February 13, 2024.
% [2] Crameri, Fabio. (2021). Scientific colour maps (7.0.1). Zenodo. https://doi.org/10.5281/zenodo.5501399

% Best drag figure for the experiment: Ncuas, m=3, 1D Temp run evaluated at 
% lambda_best=0.5
% Figure1. Best basal drag map with zoom-in boxes 
% Figure2. Best basal drag coefficient map with zoom-in boxes 
% Figure3. Sliding ration v_b/v_s 

clearvars
close all
tic

%% DATA PATHS 
issmscriptpath='/albedo/work/projects/landice/ISSM/ISSM_installations/ISSM_v4p22/trunk/';
functionpath='~/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Functions';
slancm='~/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Functions/slanCM';
addpath(genpath(issmscriptpath),functionpath,slancm)
% File names:
% use the lambda_best model result! 
lcurvefile='/albedo/scratch/user/lhoeyns/1_ISSM_execution/Zwischenspeicher/OldLcurveRange/Run_1-1/WAIS.dragInv_HO_m3_Ncuas_1DTemp_LambdaBest-0.5_22-Jan-2024.mat'; % cell array of strings, include trailing ".mat", first one is reference
figname='/albedo/home/lhoeyns/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Plotscripts/Figurescripts-Manuscript/Zenodo/BestDragFigure_v1_final.png';
figname_coeff='/albedo/home/lhoeyns/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Plotscripts/Figurescripts-Manuscript/Zenodo/BestDragCoeffFigure_v1_final.png';
figname_ratio='/albedo/home/lhoeyns/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Plotscripts/Figurescripts-Manuscript/Zenodo/BestDrag_SlidingRatioFigure_v1_final.png';

%% PARAMETERS 
% Map limits: 
%xlims=[-2000,600];                    % km
%ylims=[-1400,600];                    % km
% Pine Island Glacier zoom-in
xlims_pi=[-1760,-1400];                % km
ylims_pi=[-310,-60];                   % km
% Thwaites Glacier zoom-in
xlims_thw=[-1590,-1230];               % km
ylims_thw=[-600,-320];                 % km
% Siple Coast zoom-in
xlims_sip=[-990,-200];                 % km
ylims_sip=[-950,-370];                 % km
% Axes limits:
xtick=500;                             % km
ytick=500;                             % km
% Caxis limits:
dragcoeffclims=[0,0.6];                % [1x2] 
draglims=[0,200];                      % [1x2] kPa
dragtick=50;                           
clims=[0,100];                         % [1x2] %


% Aesthetic settings:
% Colormaps:
dragcoeffcmap=crameri('-RomaO');       % valid colormap string
flipdragcoeffcmap=0;                   % logical
dragcmap=slanCM('twilight_s');         % valid colormap string
flipdragcmap=0;                        % logical
cmap=crameri('-batlowK');              % valid colormap string

% Fontsizes:
axesfontsize=16;                       % points
axesfontsize2=22;                      % points
titlefontsize=16;                      % points
colorbarfontsize=16;                   % points
colorbarfontsize2=22;                  % points
titlefontsize2=16;                     % points
titlefontsize3=22;                     % points
% Linewidths:
axeslinewidth=2;                       % points
axeslinewidth2=1;                      % points

% Material parameters:
secondsperyear=60*60*24*365.25;

% Output size and resolution:
TileSpacing='compact';%'tight';
Padding='compact';
resolution=300;                        % dpi

%% PREPARATION  

% Load basic model:
load(lcurvefile,'md','MyInversionParameters')

% Unpack inversion parameters:
unpack(MyInversionParameters)

% 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)));
% Compute ocean and ice area fractions:
OceanAreaFraction_elements=(1-mean(OceanLevelSet_vertices(elements2d),2))/2;
IceAreaFraction_elements=(1-mean(IceLevelSet_vertices(elements2d),2))/2;

% Generate the active mask:
ActiveMask_drag=IceLevelSet_vertices<=0&OceanLevelSet_vertices>=0&md.mesh.vertexonboundary(1:numvertices2d)==0;
ActiveMask_elements=mean(ActiveMask_drag(elements2d),2);

% Extract variables we need:
if numlayers>1 && ~strcmp(flowequation,'SSA')
    ModelVel=md.results.StressbalanceSolution.Vel(md.mesh.vertexonsurface==1);
    ModelBasalVel=md.results.StressbalanceSolution.Vel(md.mesh.vertexonbase==1);
    %ObsVel=md.inversion.vel_obs(md.mesh.vertexonsurface==1);
    %OriginalDragCoeff=OriginalModel.md.results.StressbalanceSolution.FrictionCoefficient(1:numvertices2d);
    ModelDragCoeff=md.results.StressbalanceSolution.FrictionCoefficient(md.mesh.vertexonbase==1);
    if ~strcmp(nsource,'none')
        EffectivePressure_vertices=md.friction.effective_pressure(md.mesh.vertexonbase==1);
    end
else
    ModelVel=md.results.StressbalanceSolution.Vel;
    ModelBasalVel=ModelVel;
    ObsVel=md.inversion.vel_obs;
    ModelDragCoeff=md.results.StressbalanceSolution.FrictionCoefficient;
    if ~strcmp(nsource,'none')
        EffectivePressure_vertices=md.friction.effective_pressure;
    end
end
 
% Compute drag:
Drag=(ModelDragCoeff.^2).*((ModelBasalVel/secondsperyear).^(1/m));
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;

% Compute slip ratio: 
slipratio = (ModelBasalVel./ModelVel).*100; % in percent (%)

%% FIGURE Best Drag Map  
FigH = figure('Position', get(0, 'Screensize'));
tiledlayout(3,3,'TileSpacing',TileSpacing,'Padding',Padding);

% Plot zoom-in Pine Island 
nexttile(7,[1,1])
% Draw main figure:
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Drag/1e3)

% Set up color scale and colorbar:
colormap(gca,dragcmap)
caxis(draglims)

% Set up axes:
xlim(xlims_pi)
ylim(ylims_pi)
set(gca,'LineWidth',axeslinewidth2)
set(gca,'XTick',[])
set(gca,'YTick',[])
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])
box on 
% Make title:
title('\rm (a) Zoom-in Pine Island Glacier','FontSize',titlefontsize2)  

% Plot zoom-in Thwaites
nexttile(8,[1,1])
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Drag/1e3)
% Set up color scale and colorbar:
colormap(gca,dragcmap)
caxis(draglims)
% Set up axes:
xlim(xlims_thw)
ylim(ylims_thw)
set(gca,'LineWidth',axeslinewidth2)
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])
set(gca,'XTick',[])
set(gca,'YTick',[])
box on 
% Make title:
title('\rm (b) Zoom-in Thwaites Glacier','FontSize',titlefontsize2)  

% Plot zoom-in Siple coast
nexttile(9,[1,1]) 
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Drag/1e3)
% Set up color scale and colorbar:
colormap(gca,dragcmap)
caxis(draglims)
% Set up axes:
xlim(xlims_sip)
ylim(ylims_sip)
set(gca,'LineWidth',axeslinewidth2)
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])
set(gca,'XTick',[])
set(gca,'YTick',[])
box on 
% Make title:
title('\rm (c) Zoom-in Siple Coast','FontSize',titlefontsize2)  

% Plot basal drag 
nexttile(1,[2,3])  
% Draw main figure:
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Drag/1e3)
hold on 
rectangle('Position',[xlims_pi(1) ylims_pi(1) xlims_pi(2)-xlims_pi(1) ylims_pi(2)-ylims_pi(1)],'LineWidth',axeslinewidth)
text(-1740,-270,['(',char(96+1),')'],'color','k','FontSize',titlefontsize,'HorizontalAlignment','left','FontWeight', 'bold')
rectangle('Position',[xlims_thw(1) ylims_thw(1) xlims_thw(2)-xlims_thw(1) ylims_thw(2)-ylims_thw(1)],'LineWidth',axeslinewidth)
text(-1570,-560,['(',char(96+2),')'],'color','k','FontSize',titlefontsize,'HorizontalAlignment','left','FontWeight', 'bold')
rectangle('Position',[xlims_sip(1) ylims_sip(1) xlims_sip(2)-xlims_sip(1) ylims_sip(2)-ylims_sip(1)],'LineWidth',axeslinewidth)
text(-970,-410,['(',char(96+3),')'],'color','k','FontSize',titlefontsize,'HorizontalAlignment','left','FontWeight', 'bold')

% Set up color scale and colorbar:
colormap(gca,dragcmap)
caxis(draglims)
hc=colorbar;
set(get(hc,'Ylabel'),'String','kPa','FontSize',colorbarfontsize)
set(hc,'YTick',dragtick*[ceil(draglims(1)/dragtick):1:floor(draglims(2)/dragtick)])
set(hc,'Ticks',[0 50 100 150 200],'TickLabels',{'0','50','100','150','\geq 200'})
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize,'LineWidth',axeslinewidth)
% Set up axes:
xlim(xlims/1000)
ylim(ylims/1000)
set(gca,'XTick',xtick*[ceil(xlims(1)/(1e3*xtick)):1:floor(xlims(2)/(1e3*xtick))])
set(gca,'YTick',ytick*[ceil(ylims(1)/(1e3*ytick)):1:floor(ylims(2)/(1e3*ytick))])
xlabel('x (km)','FontSize',axesfontsize)
ylabel('y (km)','FontSize',axesfontsize)
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'FontSize',axesfontsize,'LineWidth',axeslinewidth)
box on 
% Make title:
title('\rm (d) Best basal drag map, m=3, Budd, N_{CUAS}','FontSize',titlefontsize) 

% Save figure:
saveas(FigH,figname,'png');
exportgraphics(FigH,figname,'resolution',resolution)

%% FIGURE Best Drag Coeff Map  
FigH2 = figure('Position', get(0, 'Screensize'));
tiledlayout(3,3,'TileSpacing',TileSpacing,'Padding',Padding);

% plot zoom-in Pine Island 
nexttile(7,[1,1])
% Draw main figure:
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',(ModelDragCoeff.^2/(secondsperyear^(1/(2*m)))))

% Set up color scale and colorbar:
colormap(gca,dragcoeffcmap)
caxis(dragcoeffclims)
% Set up axes:
xlim(xlims_pi)
ylim(ylims_pi)
set(gca,'LineWidth',axeslinewidth2)
set(gca,'XTick',[])
set(gca,'YTick',[])
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])
box on 
title('\rm (a) Zoom-in Pine Island Glacier','FontSize',titlefontsize2)  

% plot zoom-in Thwaites
nexttile(8,[1,1])
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',(ModelDragCoeff.^2/(secondsperyear^(1/(2*m)))))
% Set up color scale and colorbar:
colormap(gca,dragcoeffcmap)
caxis(dragcoeffclims)
% Set up axes:
xlim(xlims_thw)
ylim(ylims_thw)
set(gca,'LineWidth',axeslinewidth2)
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])
set(gca,'XTick',[])
set(gca,'YTick',[])
box on 
title('\rm (b) Zoom-in Thwaites Glacier','FontSize',titlefontsize2)  

% Plot zoom-in Siple coast
nexttile(9,[1,1]) 
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',(ModelDragCoeff.^2/(secondsperyear^(1/(2*m)))))
% Set up color scale and colorbar:
colormap(gca,dragcoeffcmap)
caxis(dragcoeffclims)
% Set up axes:
xlim(xlims_sip)
ylim(ylims_sip)
set(gca,'LineWidth',axeslinewidth2)
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])
set(gca,'XTick',[])
set(gca,'YTick',[])
box on 
% Make title:
title('\rm (c) Zoom-in Siple Coast','FontSize',titlefontsize2) 

% Plot best drag coefficient map 
nexttile(1,[2,3])  
% Draw main figure:
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',(ModelDragCoeff.^2/(secondsperyear^(1/(2*m)))))
hold on 
rectangle('Position',[xlims_pi(1) ylims_pi(1) xlims_pi(2)-xlims_pi(1) ylims_pi(2)-ylims_pi(1)],'LineWidth',axeslinewidth)
text(-1740,-270,['(',char(96+1),')'],'color','k','FontSize',titlefontsize,'HorizontalAlignment','left','FontWeight', 'bold')
rectangle('Position',[xlims_thw(1) ylims_thw(1) xlims_thw(2)-xlims_thw(1) ylims_thw(2)-ylims_thw(1)],'LineWidth',axeslinewidth)
text(-1570,-560,['(',char(96+2),')'],'color','k','FontSize',titlefontsize,'HorizontalAlignment','left','FontWeight', 'bold')
rectangle('Position',[xlims_sip(1) ylims_sip(1) xlims_sip(2)-xlims_sip(1) ylims_sip(2)-ylims_sip(1)],'LineWidth',axeslinewidth)
text(-970,-410,['(',char(96+3),')'],'color','k','FontSize',titlefontsize,'HorizontalAlignment','left','FontWeight', 'bold')

% Generate drag coeff units string:
if ~strcmp(nsource,'none')
    if r==1
        if m==1
            cunitsstring='(m/yr)^{-1}';
        else
            cunitsstring=['(m/yr)^{-1/',num2str(m),'}'];
        end
    else
        if m==1
            cunitsstring=['Pa^{',num2str(1-r),'}(m/yr)^{-1}'];
        else
            cunitsstring=['Pa^{',num2str(1-r),'}(m/yr)^{-1/',num2str(m),'}'];
        end
    end
else
    if m==1
        cunitsstring='Pa(m/yr)^{-1}';
    else
        cunitsstring=['Pa(m/yr)^{-1/',num2str(m),'}'];
    end
end

% Set up color scale and colorbar:
colormap(gca,dragcoeffcmap)
caxis(dragcoeffclims)
hc=colorbar;
%set(hc,'YTick',ccticks)
%set(hc,'YTickLabel',cclabels)
set(get(hc,'Ylabel'),'String',cunitsstring,'FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize,'LineWidth',axeslinewidth)
set(hc,'Ticks',[0 0.1 0.2 0.3 0.4 0.5 0.6],'TickLabels',{'0','0.1','0.2','0.3','0.4','0.5','\geq 0.6'})

% Set up axes:
xlim(xlims/1000)
ylim(ylims/1000)
set(gca,'XTick',xtick*[ceil(xlims(1)/(1e3*xtick)):1:floor(xlims(2)/(1e3*xtick))])
set(gca,'YTick',ytick*[ceil(ylims(1)/(1e3*ytick)):1:floor(ylims(2)/(1e3*ytick))])
%set(gca,'YTickLabel',[])
xlabel('x (km)','FontSize',axesfontsize)
ylabel('y (km)','FontSize',axesfontsize)
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'FontSize',axesfontsize,'LineWidth',axeslinewidth)
box on 
% Make title:
title('\rm (d) Best basal drag coefficient map, m=3, Budd, N_{CUAS}','FontSize',titlefontsize) 

% Save figure:
saveas(FigH2,figname_coeff,'png');
exportgraphics(FigH2,figname_coeff,'resolution',resolution)

%% FIGURE Sliding ratio 
FigH3 = figure('Position', get(0, 'Screensize'));
% Make figure:
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d,'EdgeColor','none','FaceColor','interp','FaceVertexCData',slipratio)

% Set up color scale and colorbar:
caxis(clims)  
colormap(gca,cmap)
hc=colorbar;
set(get(hc,'Ylabel'),'String','v_b/v_s (%)','FontSize',colorbarfontsize2)
set(hc,'TickDir','out')
set(hc,'Ticks',[0 20 40 60 80 100],'TickLabels',{'0','20','40','60','80','100'})
set(hc,'FontSize',colorbarfontsize2,'LineWidth',axeslinewidth)
% Set up axes:
xlim(xlims/1000)
ylim(ylims/1000)

set(gca,'XTick',xtick*[ceil(xlims(1)/(1e3*xtick)):1:floor(xlims(2)/(1e3*xtick))])
set(gca,'YTick',ytick*[ceil(ylims(1)/(1e3*ytick)):1:floor(ylims(2)/(1e3*ytick))])
xlabel('x (km)','FontSize',axesfontsize2)
ylabel('y (km)','FontSize',axesfontsize2)
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'FontSize',axesfontsize2,'LineWidth',axeslinewidth)
box on 
% Make title:
title('\rm Slip ratio, m=3, Budd, N_{CUAS}','FontSize',titlefontsize3) 

% Save figure:
saveas(FigH3,figname_ratio,'png');
exportgraphics(FigH3,figname_ratio,'resolution',resolution)

disp('Done!')
toc