% ThermalFigure_v1

% Lea-Sophie Hoeyns, 11/30/2023

% This script makes a figure of the 1D steady-state thermal model output of  
% both depth-averaged rheological stiffness (Bbar) and the basal melt rates for publication. 

% 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

% a. Depth-averaged rheology Bbar; b. Basal melt rates. 

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:
gridfile= '/albedo/home/lhoeyns/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/AntarcticaCombinedGrids_v4.mat';
boundaryfile='~/1_Models/WAIS_FrictionOptimization_v2/Data/outline/BoundaryPicking_WAIS-with-iceshelves_50km_BedMachineGrid_new_pick_bearbeitet.exp';
% SSA run with Nocp or Ncuas and 1DTemp
modelfile='/albedo/scratch/user/lhoeyns/1_ISSM_execution/Zwischenspeicher/OldLcurveRange/Run_1/WAIS.dragInv_HO_m3_Ncuas_1DTemp_Lcurve_17-Jan-2024.mat';
thermalfile='~/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/ThermalOutput.mat';
nfile='~/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/spinup10years_1000m_ROI_WAIS_last_cropped.nc';
%nfile='/albedo/work/projects/p_modsem23/shared_inputs/Ant1DMelt_spinup40years_4km_Sy0.4_Tvostok0.1_output_cropped_last.nc';
figname='/albedo/home/lhoeyns/1_Models/WAIS_FrictionOptimization_v2/Scripts_MW/Plotscripts/Figurescripts-Manuscript/Zenodo/ThermalFigure_v1_final.png'; 

%% PARAMETERS 
% Map limits:
xlims=[-2000,600];                     % km
ylims=[-1400,600];                     % km
% Map cropping limits: 
xlims_crop=[-2000,600]*1e3;            % [1x2] m
ylims_crop=[-1400,600]*1e3;            % [1x2] m
% Axes limits:s
xtick=500;                             % km
ytick=500;                             % km
% Caxis limits:
bbarlims=[0.4,1];                      % [1x2] MPa*yr^-1/3
bbartick=0.1;                          
meltlims=[10^-3,5.001];                % [1x2] m/yr
meltticks=0.5;

% Aesthetic settings
% Colormaps:
bbarcmap=flipud(slanCM('romao'));
meltcmap=crameri('-davos',10);
% Colors:
boundarylinecolor='k';                 % valid color string or RGB
gray=[.7 .7 .7];
maskcolor='k';                         % valid color string or RGB

% Fontsizes:
axesfontsize=16;                       % points
titlefontsize=16;                      % points
colorbarfontsize=16;                   % points
% Linewidths: 
axeslinewidth=2;                       % points
boundarylinewidth=3;                   % points
masklinewidth=1.2;                     % points

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

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

%% PREPARATION 
% Load domain boundary:
DomainBoundary=expread(boundaryfile);
% Load grids:
load(gridfile,'X','Y','BedElev','SurfElev','IceThick','Mask','stereoprojection')

% Identify grid bounds within the figure:
ind1_grid=find(Y>=ylims(1)*1e3,1,'first');
ind2_grid=find(Y<=ylims(2)*1e3,1,'last');
ind3_grid=find(X>=xlims(1)*1e3,1,'first');
ind4_grid=find(X<=xlims(2)*1e3,1,'last');

% Rename coordinate variables:
X_grid=X;
Y_grid=Y;

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

% Evaluate grid size:
[ysize,xsize]=size(Mask);
dx=X_grid(2)-X_grid(1); % m

% Generate contour structure for mask:
%C_main=contourc(X/1000,Y/1000,Mask,[.5,1.5,2.5]);
C_main=contourc(X_grid(ind3_grid:ind4_grid)/1000,Y_grid(ind1_grid:ind2_grid)/1000,Mask(ind1_grid:ind2_grid,ind3_grid:ind4_grid),[.5,1.5]);%,2.5]);

% Load Model: 
thisfile=load(modelfile,'BasicModel');
md=thisfile.BasicModel;
md=FlattenModelStructure(md,'expand');

% 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);
% 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;

%% LOAD DATA 

% Load temperature and basal melt rates: 
% load(thermalfile);
% Temperature=temp(1:numvertices2d)-273.15;
% Meltrate=meltrate(1:numvertices2d);

% Load basal melt rates from nfile:
loadnetcdf(nfile,'x','y','watersource'); 
watersource=watersource';
MeltRate_vertices=MultiWavelengthInterpolator(x,y,watersource,x2d,y2d,elements2d);
MeltRate_elements=mean(MeltRate_vertices(elements2d),2);
MeltRate_vertices=MeltRate_vertices*secondsperyear;
MeltRateRate_vertices(MeltRate_vertices<10^-4)=10^-4;

% Load rheology B and compute depth-averaged Bbar
InitialBbar_vertices=DepthAverage(md,md.materials.rheology_B);
InitialBbar_vertices=LaplacianInterpolation_mesh(x2d,y2d,elements2d,InitialBbar_vertices);
md=FlattenModelStructure(md,'expand');
md=collapse(md);
md.materials.rheology_B=InitialBbar_vertices;
Rheology=md.materials.rheology_B;

%% MAKE FIGURE
FigH = figure('Position', get(0, 'Screensize'));
tiledlayout(1,2,'TileSpacing',TileSpacing,'Padding',Padding)

%% MAKE SUBFIGURE temperature
% t1=nexttile;
% % Plot Temperature:
% IceFraction_elements=mean(md.mask.ice_levelset(elements2d)<0,2);
% patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceFraction_elements>.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Temperature)
% 
% %patch('Vertices',[md.mesh.x/1000,md.mesh.y/1000],'Faces',md.mesh.elements(IceFraction_elements>.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Temperature)
% hold on %OceanAreaFraction_elements<.5
% % % Plot grounding line and ice front:
% % startind=1;
% % done=0;
% % while done==0
% %     stopind=startind+C_main(2,startind);
% %     plot3(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),(topz+2)*ones(size(C_main(1,startind+1:stopind))),'Color',maskcolor,'LineWidth',masklinewidth)
% %     startind=stopind+1;
% %     if startind>size(C_main,2)
% %         done=1;
% %     end
% % end
% % Plot grounding line and ice front:
% startind=1;
% done=0;
% while done==0
%     stopind=startind+C_main(2,startind);
%     plot(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),'Color',maskcolor,'LineWidth',masklinewidth)
%     startind=stopind+1;
%     if startind>size(C_main,2)
%         done=1;
%     end
% end
% % Plot domain outline:
% plot(DomainBoundary.x/1000,DomainBoundary.y/1000,'Color',boundarylinecolor,'LineWidth',boundarylinewidth)
% % % Plot domain outline:
% % plot3(DomainBoundary.x/1000,DomainBoundary.y/1000,(topz+3)*ones(size(DomainBoundary.x)),'Color',boundarylinecolor,'LineWidth',boundarylinewidth)
% 
% % Set up color scale and colorbar:
% colormap(t1,tempcmap)
% caxis(templims)
% hc=colorbar;
% set(hc,'YTick',temptick*[ceil(templims(1)/temptick):1:floor(templims(2)/temptick)])
% %set(hc,'YTick',[0 10 20 30])
% set(get(hc,'Ylabel'),'String','C','FontSize',colorbarfontsize)
% 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)/xtick):1:floor(xlims(2)/xtick)])
% set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
% %set(gca,'XTickLabel',[])
% %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)
% set(gca,'ydir','normal')
% %set(gca,'TickDir','out')
% box on 
% % Make title:
% title('\rm (a) Basal temperature','FontSize',titlefontsize)

%% MAKE SUBFIGURE depth-averaged rheology (Bbar)
t1=nexttile;

% Plot Rheology:
IceFraction_elements=mean(md.mask.ice_levelset(elements2d)<0,2);
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceFraction_elements>.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Rheology*1e-6*(secondsperyear^-(1/3))) % MPa

hold on 
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),'Color',maskcolor,'LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Plot domain outline:
plot(DomainBoundary.x/1000,DomainBoundary.y/1000,'Color',boundarylinecolor,'LineWidth',boundarylinewidth)

% Set up color scale and colorbar:
colormap(t1,bbarcmap)
caxis(bbarlims)
hc=colorbar;
%set(hc,'YTick',bbartick*[ceil(bbarlims(1)/bbartick):1:floor(bbarlims(2)/bbartick)]
hc.Ticks = [0.4,0.6,0.8,1,1.2];
set(get(hc,'Ylabel'),'String','MPa yr^{-1/3}','FontSize',colorbarfontsize)
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)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
%set(gca,'XTickLabel',[])
%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)
set(gca,'ydir','normal')
%set(gca,'TickDir','out')
box on 
% Make title:
title('\rm (a) Rheological stiffness (B)','FontSize',titlefontsize)

%% MAKE SUBFIGURE effective pressure Ncuas 
t2=nexttile;

% Plot basal melt rates:
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',MeltRate_vertices)
hold on
patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(OceanAreaFraction_elements>.5,:),'EdgeColor','none','FaceColor',gray,'FaceVertexCData',MeltRate_vertices)
hold on 
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),'Color','k','LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Plot domain outline:
plot(DomainBoundary.x/1000,DomainBoundary.y/1000,'Color',boundarylinecolor,'LineWidth',boundarylinewidth)

% Set up color scale and colorbar:
thiscmap=colormap(t2,meltcmap);
set(gca,'colorscale','log');
caxis(meltlims)
% Make colorbar:
hc=colorbar;
set(get(hc,'Ylabel'),'String','m/yr','FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize,'LineWidth',axeslinewidth)
%hc.Ticks=[10^-3,2.39*10^-3,5.5*10^-3,1.3*10^-2,3*10^-2,7*10^-2,0.17,0.38,0.93,2.2,5];
%hc.TickLabels={10^-3,2*10^-3,5*10^-3,1.5*10^-2,5.0*10^-2,0.100,0.200,0.400,0.900,2.000,5.000};
hc.Ticks=[10^-3,5.5*10^-3,3*10^-2,0.17,0.93,5];
hc.TickLabels={10^-3,5*10^-3,5.0*10^-2,0.200,0.900,5.000};
hc.Ruler.MinorTick = 'off';
% Set up axes:
xlim(xlims/1000)
ylim(ylims/1000)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
% set(gca,'XTickLabel',[])
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)
set(gca,'ydir','normal')
% set(gca,'TickDir','out')
box on 
% Make title:
title('\rm (b) Basal melt rate','FontSize',titlefontsize)

%% SAVE FIGURE 
saveas(FigH,figname,'png');
exportgraphics(FigH,figname,'resolution',resolution)

disp('Done!')
toc