%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB Script to run and compare results of the MISMIP+ experiment 
% (Cornford et al. 2020) at different spatial resolutions, and with and
% without applying a flux condition to calculate grounding line migrations
%
% Author: Violaine Coulon (violaine.coulon@ulb.be)
% Date: October 2023
%
% Data Source:
% The data used in this script are available at 10.5281/zenodo.8398772
%
% MATLAB Version and Toolboxes:
% This script is designed to work with MATLAB R20XX or later and does not
% require any additional toolboxes beyond the standard MATLAB functions.
% Some functions used in the script might be available in earlier versions
% of MATLAB as well.
%
% Please make sure you have the necessary MATLAB version installed to run
% this script.
%
% Usage:
% 1. Ensure you have the necessary MATLAB version and toolboxes installed.
% 2. Download the paper and the required data from the provided sources.
% 3. Run this script 
%
% Citation:
% If you use this script or the reproduced figure in your work, please
% consider citing the original paper as follows:
%
% Coulon, V., Klose, A. K., Kittel, C., Edwards, T., Turner, F., Winkelmann, R., and Pattyn, F.: 
% Disentangling the drivers of future Antarctic ice loss with a historically-calibrated ice-sheet model, 
% EGUsphere [preprint], https://doi.org/10.5194/egusphere-2023-1532, 2023.
%
% Disclaimer:
% This script is provided as-is without any warranties. The author is not
% responsible for any inaccuracies or errors in the reproduced figure.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

run_mismip=1;
plot_fig=1;

if run_mismip==1
    MISMIPplus(1,0)
    close all
    MISMIPplus(2,0)
    close all
    MISMIPplus(2,1)
    close all
    MISMIPplus(4,0)
    close all
    MISMIPplus(4,1)
    close all
    MISMIPplus(8,0)
    close all
    MISMIPplus(8,1)
    close all
end

if plot_fig==1

    % plot GL position

    name='1km_schoof0';
    ctr.delta=1e3;
    ctr.imax=43; % need number of cells + 2
    ctr.jmax=702; % need number of cells + 1
    Li=(ctr.imax-2)*ctr.delta;
    Lj=(ctr.jmax-2)*ctr.delta;
    [X,Y]=meshgrid(-ctr.delta:ctr.delta:Lj,-ctr.delta:ctr.delta:Li);

    figure('Units','centimeters','Position', [5 5 25 20])

    load(['mismiptest1_',name]);
    Li=(ctr.imax-2)*ctr.delta/1e3;
    Lj=(ctr.jmax-2)*ctr.delta/1e3;
    [Xo,Yo] = meshgrid(-ctr.delta/1e3:ctr.delta/1e3:Lj,-ctr.delta/1e3:ctr.delta/1e3:Li);
    contour(Xo,-Yo,MASK,[1],'k','linewidth',4.5);
    hold on;
    load(['mismiptest1r_',name]);
    contour(Xo,-Yo,MASK,[1],'b','linewidth',4.5,'HandleVisibility','off');
    load(['mismiptest1rr_',name]);
    contour(Xo,-Yo,MASK,[1],'color',[255 165 0]/255,'linewidth',4.5,'HandleVisibility','off');
    load(['mismiptest1ra_',name]);
    contour(Xo,-Yo,MASK,[1],'r','linewidth',4.5,'HandleVisibility','off');
    xlim([300 560]);
    xlabel('x (km)');
    ylabel('y (km)');
    grid on;
    %legend('Ice0','Ice1r','Ice1rr','Ice1ra');

    name='2km_schoof0';

    ctr.delta=2e3;
    ctr.imax=23; % need number of cells + 2
    ctr.jmax=352; % need number of cells + 1
    Li=(ctr.imax-2)*ctr.delta;
    Lj=(ctr.jmax-2)*ctr.delta;
    [X,Y]=meshgrid(-ctr.delta:ctr.delta:Lj,-ctr.delta:ctr.delta:Li);

    load(['mismiptest1_',name]);
    Li=(ctr.imax-2)*ctr.delta/1e3;
    Lj=(ctr.jmax-2)*ctr.delta/1e3;
    [Xo,Yo] = meshgrid(-ctr.delta/1e3:ctr.delta/1e3:Lj,-ctr.delta/1e3:ctr.delta/1e3:Li);
    contour(Xo,-Yo,MASK,[1],'k','linewidth',3);
    hold on;
    load(['mismiptest1r_',name]);
    contour(Xo,-Yo,MASK,[1],'b','linewidth',3,'HandleVisibility','off');
    load(['mismiptest1rr_',name]);
    contour(Xo,-Yo,MASK,[1],'color',[255 165 0]/255,'linewidth',3,'HandleVisibility','off');
    load(['mismiptest1ra_',name]);
    contour(Xo,-Yo,MASK,[1],'r','linewidth',3,'HandleVisibility','off');
    xlim([300 560]);
    xlabel('x (km)');
    ylabel('y (km)');
    grid on;

    name='2km_schoof1';

    load(['mismiptest1_',name]);
    Li=(ctr.imax-2)*ctr.delta/1e3;
    Lj=(ctr.jmax-2)*ctr.delta/1e3;
    [Xo,Yo] = meshgrid(-ctr.delta/1e3:ctr.delta/1e3:Lj,-ctr.delta/1e3:ctr.delta/1e3:Li);
    contour(Xo,-Yo,MASK,[1],'k','linewidth',2,'LineStyle','--');
    hold on;
    load(['mismiptest1r_',name]);
    contour(Xo,-Yo,MASK,[1],'b','linewidth',2,'LineStyle','--','HandleVisibility','off');
    load(['mismiptest1rr_',name]);
    contour(Xo,-Yo,MASK,[1],'color',[255 165 0]/255,'linewidth',2,'LineStyle','--','HandleVisibility','off');
    load(['mismiptest1ra_',name]);
    contour(Xo,-Yo,MASK,[1],'r','linewidth',2,'LineStyle','--','HandleVisibility','off');
    xlim([300 560]);
    xlabel('x (km)');
    ylabel('y (km)');
    grid on;

    name='4km_schoof0';
    ctr.delta=4e3;
    ctr.imax=13; % need number of cells + 2
    ctr.jmax=177; % need number of cells + 1
    Li=(ctr.imax-2)*ctr.delta;
    Lj=(ctr.jmax-2)*ctr.delta;
    [X,Y]=meshgrid(-ctr.delta:ctr.delta:Lj,-ctr.delta:ctr.delta:Li);

    load(['mismiptest1_',name]);
    Li=(ctr.imax-2)*ctr.delta/1e3;
    Lj=(ctr.jmax-2)*ctr.delta/1e3;
    [Xo,Yo] = meshgrid(-ctr.delta/1e3:ctr.delta/1e3:Lj,-ctr.delta/1e3:ctr.delta/1e3:Li);
    contour(Xo,-Yo,MASK,[1],'k','linewidth',2);
    hold on;
    load(['mismiptest1r_',name]);
    contour(Xo,-Yo,MASK,[1],'b','linewidth',2,'HandleVisibility','off');
    load(['mismiptest1rr_',name]);
    contour(Xo,-Yo,MASK,[1],'color',[255 165 0]/255,'linewidth',2,'HandleVisibility','off');
    load(['mismiptest1ra_',name]);
    contour(Xo,-Yo,MASK,[1],'r','linewidth',2,'HandleVisibility','off');
    xlim([300 560]);
    xlabel('x (km)');
    ylabel('y (km)');
    grid on;

    name='4km_schoof1';

    load(['mismiptest1_',name]);
    Li=(ctr.imax-2)*ctr.delta/1e3;
    Lj=(ctr.jmax-2)*ctr.delta/1e3;
    [Xo,Yo] = meshgrid(-ctr.delta/1e3:ctr.delta/1e3:Lj,-ctr.delta/1e3:ctr.delta/1e3:Li);
    contour(Xo,-Yo,MASK,[1],'k','linewidth',2,'LineStyle',':');
    load(['mismiptest1r_',name]);
    contour(Xo,-Yo,MASK,[1],'b','linewidth',2,'LineStyle',':','HandleVisibility','off');
    load(['mismiptest1rr_',name]);
    contour(Xo,-Yo,MASK,[1],'color',[255 165 0]/255,'linewidth',2,'LineStyle',':','HandleVisibility','off');
    load(['mismiptest1ra_',name]);
    contour(Xo,-Yo,MASK,[1],'r','linewidth',2,'LineStyle',':','HandleVisibility','off');
    xlim([300 560]);
    xlabel('x (km)','FontSize',24);
    ylabel('y (km)','FontSize',24);
    grid on;

    name='8km_schoof0';

    ctr.delta=8e3;
    ctr.imax=8; % need number of cells + 2
    ctr.jmax=89; % need number of cells + 1
    Li=(ctr.imax-2)*ctr.delta;
    Lj=(ctr.jmax-2)*ctr.delta;
    [X,Y]=meshgrid(-ctr.delta:ctr.delta:Lj,-ctr.delta:ctr.delta:Li);

    load(['mismiptest1_',name]);
    Li=(ctr.imax-2)*ctr.delta/1e3;
    Lj=(ctr.jmax-2)*ctr.delta/1e3;
    [Xo,Yo] = meshgrid(-ctr.delta/1e3:ctr.delta/1e3:Lj,-ctr.delta/1e3:ctr.delta/1e3:Li);
    contour(Xo,-Yo,MASK,[1],'k','linewidth',1,'LineStyle','-');
    load(['mismiptest1r_',name]);
    contour(Xo,-Yo,MASK,[1],'b','linewidth',1,'LineStyle','-','HandleVisibility','off');
    load(['mismiptest1rr_',name]);
    contour(Xo,-Yo,MASK,[1],'color',[255 165 0]/255,'linewidth',1,'LineStyle','-','HandleVisibility','off');
    load(['mismiptest1ra_',name]);
    contour(Xo,-Yo,MASK,[1],'r','linewidth',1,'LineStyle','-','HandleVisibility','off');
    xlim([300 560]);
    xlabel('x (km)','FontSize',18);
    ylabel('y (km)','FontSize',18);
    grid on;

    name='8km_schoof1';

    load(['mismiptest1_',name]);
    Li=(ctr.imax-2)*ctr.delta/1e3;
    Lj=(ctr.jmax-2)*ctr.delta/1e3;
    [Xo,Yo] = meshgrid(-ctr.delta/1e3:ctr.delta/1e3:Lj,-ctr.delta/1e3:ctr.delta/1e3:Li);
    contour(Xo,-Yo,MASK,[1],'k','linewidth',1,'LineStyle','-.');
    load(['mismiptest1r_',name]);
    contour(Xo,-Yo,MASK,[1],'b','linewidth',1,'LineStyle','-.','HandleVisibility','off');
    load(['mismiptest1rr_',name]);
    contour(Xo,-Yo,MASK,[1],'color',[255 165 0]/255,'linewidth',1,'LineStyle','-.','HandleVisibility','off');
    load(['mismiptest1ra_',name]);
    contour(Xo,-Yo,MASK,[1],'r','linewidth',1,'LineStyle','-.','HandleVisibility','off');
    xlim([300 560]);
    xlabel('x (km)','FontSize',18);
    ylabel('y (km)','FontSize',18);
    grid on;

    ylim([-40 0])

    legend({'1 km', '2 km', '2km - flux condition', '4 km', '4 km - flux condition', '8 km','8 km - flux condition'},'Location','SouthWest','EdgeColor','none','Color','none','FontSize',18)
    set(gca,'Units','centimeters','Position', [2 2 22.3 16.64],'FontSize',18)

    % Plot GL vs time

    name='1km_schoof0';
    ctr.delta=1e3;
    ctr.imax=45; % need number of cells + 2
    ctr.jmax=701; % need number of cells + 1
    Li=(ctr.imax-2)*ctr.delta;
    Lj=(ctr.jmax-2)*ctr.delta;
    [X,Y]=meshgrid(-ctr.delta:ctr.delta:Lj,-ctr.delta:ctr.delta:Li);

    figure('Units','centimeters','Position', [5 5 25 20])
    % Ice0
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1_ice0_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        if i==1
            plot(time(i),X(1,z)/1000,'s','color','k','MarkerFaceColor','k','MarkerSize',8); hold on
        else
            plot(time(i),X(1,z)/1000,'s','color','k','HandleVisibility','off','MarkerFaceColor','k','MarkerSize',8); hold on
        end
    end
    % 1rr
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1rr_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'s','color',[255 165 0]/255,'HandleVisibility','off','MarkerFaceColor',[255 165 0]/255,'MarkerSize',8); hold on
    end
    % 1r
    time=0:10:100;
    for i=1:11
        load(['mismiptest1r_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'s','color','b','HandleVisibility','off','MarkerFaceColor','b','MarkerSize',8); hold on
    end
    % 1ra
    time=[100:10:200];
    for i=1:11
        load(['mismiptest1ra_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'s','color','r','HandleVisibility','off','MarkerFaceColor','r','MarkerSize',8); hold on
    end
    ylim([320 460])

    name='2km_schoof0';

    ctr.delta=2e3;
    ctr.imax=23; % need number of cells + 2
    ctr.jmax=352; % need number of cells + 1
    Li=(ctr.imax-2)*ctr.delta;
    Lj=(ctr.jmax-2)*ctr.delta;
    [X,Y]=meshgrid(-ctr.delta:ctr.delta:Lj,-ctr.delta:ctr.delta:Li);

    % Ice0
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1_ice0_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        if i==1
            plot(time(i),X(1,z)/1000,'o','color','k','MarkerFaceColor','k','MarkerSize',8); hold on
        else
            plot(time(i),X(1,z)/1000,'o','color','k','HandleVisibility','off','MarkerFaceColor','k','MarkerSize',8); hold on
        end
    end
    % 1rr
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1rr_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'o','color',[255 165 0]/255,'HandleVisibility','off','MarkerFaceColor',[255 165 0]/255,'MarkerSize',8); hold on
    end
    % 1r
    time=0:10:100;
    for i=1:11
        load(['mismiptest1r_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'o','color','b','HandleVisibility','off','MarkerFaceColor','b','MarkerSize',8); hold on
    end
    % 1ra
    time=[100:10:200];
    for i=1:11
        load(['mismiptest1ra_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'o','color','r','HandleVisibility','off','MarkerFaceColor','r','MarkerSize',8); hold on
    end
    ylim([320 460])

    name='2km_schoof1';

    % Ice0
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1_ice0_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        if i==1
            plot(time(i),X(1,z)/1000,'o','color','k','MarkerSize',8,'LineWidth',1); hold on
        else
            plot(time(i),X(1,z)/1000,'o','color','k','HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
        end
    end
    % 1rr
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1rr_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'o','color',[255 165 0]/255,'HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
    end
    % 1r
    time=0:10:100;
    for i=1:11
        load(['mismiptest1r_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'o','color','b','HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
    end
    % 1ra
    time=[100:10:200];
    for i=1:11
        load(['mismiptest1ra_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'o','color','r','HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
    end
    ylim([320 460])

    name='4km_schoof0';
    ctr.delta=4e3;
    ctr.imax=13; % need number of cells + 2
    ctr.jmax=177; % need number of cells + 1
    Li=(ctr.imax-2)*ctr.delta;
    Lj=(ctr.jmax-2)*ctr.delta;
    [X,Y]=meshgrid(-ctr.delta:ctr.delta:Lj,-ctr.delta:ctr.delta:Li);

    % Ice0
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1_ice0_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        if i==1
            plot(time(i),X(1,z)/1000,'^','color','k','MarkerFaceColor','k','MarkerSize',8); hold on
        else
            plot(time(i),X(1,z)/1000,'^','color','k','HandleVisibility','off','MarkerFaceColor','k','MarkerSize',8); hold on
        end
    end
    % 1rr
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1rr_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'^','color',[255 165 0]/255,'HandleVisibility','off','MarkerFaceColor',[255 165 0]/255,'MarkerSize',8); hold on
    end
    % 1r
    time=0:10:100;
    for i=1:11
        load(['mismiptest1r_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'^','color','b','HandleVisibility','off','MarkerFaceColor','b','MarkerSize',8); hold on
    end
    % 1ra
    time=[100:10:200];
    for i=1:11
        load(['mismiptest1ra_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'^','color','r','HandleVisibility','off','MarkerFaceColor','r','MarkerSize',8); hold on
    end
    ylim([320 460])

    name='4km_schoof1';

    % Ice0
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1_ice0_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        if i==1
            plot(time(i),X(1,z)/1000,'^','color','k','MarkerSize',8,'LineWidth',1); hold on
        else
            plot(time(i),X(1,z)/1000,'^','color','k','HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
        end
    end
    % 1rr
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1rr_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'^','color',[255 165 0]/255,'HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
    end
    % 1r
    time=0:10:100;
    for i=1:11
        load(['mismiptest1r_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'^','color','b','HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
    end
    % 1ra
    time=[100:10:200];
    for i=1:11
        load(['mismiptest1ra_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'^','color','r','HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
    end
    ylim([320 460])

    name='8km_schoof0';
    ctr.delta=8e3;
    ctr.imax=8; % need number of cells + 2
    ctr.jmax=89; % need number of cells + 1
    Li=(ctr.imax-2)*ctr.delta;
    Lj=(ctr.jmax-2)*ctr.delta;
    [X,Y]=meshgrid(-ctr.delta:ctr.delta:Lj,-ctr.delta:ctr.delta:Li);

    % Ice0
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1_ice0_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        if i==1
            plot(time(i),X(1,z)/1000,'pentagram','color','k','MarkerFaceColor','k','MarkerSize',8); hold on
        else
            plot(time(i),X(1,z)/1000,'pentagram','color','k','HandleVisibility','off','MarkerFaceColor','k','MarkerSize',8); hold on
        end
    end
    % 1rr
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1rr_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'pentagram','color',[255 165 0]/255,'HandleVisibility','off','MarkerFaceColor',[255 165 0]/255,'MarkerSize',8); hold on
    end
    % 1r
    time=0:10:100;
    for i=1:11
        load(['mismiptest1r_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'pentagram','color','b','HandleVisibility','off','MarkerFaceColor','b','MarkerSize',8); hold on
    end
    % 1ra
    time=[100:10:200];
    for i=1:11
        load(['mismiptest1ra_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'pentagram','color','r','HandleVisibility','off','MarkerFaceColor','r','MarkerSize',8); hold on
    end
    ylim([320 460])


    name='8km_schoof1';

    % Ice0
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1_ice0_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        if i==1
            plot(time(i),X(1,z)/1000,'pentagram','color','k','MarkerSize',8,'LineWidth',1); hold on
        else
            plot(time(i),X(1,z)/1000,'pentagram','color','k','HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
        end
    end
    % 1rr
    time=[0:10:200];
    for i=1:21
        load(['mismiptest1rr_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'pentagram','color',[255 165 0]/255,'HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
    end
    % 1r
    time=0:10:100;
    for i=1:11
        load(['mismiptest1r_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'pentagram','color','b','HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
    end
    % 1ra
    time=[100:10:200];
    for i=1:11
        load(['mismiptest1ra_',name,'_',num2str(i-1,'%03i')])
        z=find(glMASK(1,:)==2);
        plot(time(i),X(1,z)/1000,'pentagram','color','r','HandleVisibility','off','MarkerSize',8,'LineWidth',1); hold on
    end
    ylim([320 460])

    xlabel('Time (a)','FontSize',24)
    ylabel('Grounding line position, x_g(y=0) (km) ','FontSize',24)
    legend({'1 km', '2 km', '2km - flux condition', '4 km', '4 km - flux condition', '8 km', '8 km - flux condition'},'Location','SouthWest','EdgeColor','none','Color','none','Fontsize',18)
    set(gca,'Units','centimeters','Position', [2 2 22.3 16.64],'FontSize',18)
    set(gca,'color','none')

end