% contourf_insolation_timexlon
clc;clear
new_jet=[0,0.400000000000000,0.400000000000000;0,0.600000000000000,0.600000000000000;0,0.800000000000000,0.800000000000000;0,1,1;0.200000000000000,1,1;0.600000000000000,1,1;0.796078431372549,1,1;0.898039215686275,1,1;1,0.898039215686275,0.796078431372549;1,0.792156862745098,0.600000000000000;1,0.678431372549020,0.396078431372549;1,0.556862745098039,0.200000000000000;1,0.431372549019608,0;0.800000000000000,0.333333333333333,0;0.600000000000000,0.239215686274510,0;0.400000000000000,0.152941176470588,0];
new_jet=[0,0,0.588235294117647;0,0,0.803921568627451;0.254901960784314,0.411764705882353,0.882352941176471;0.117647058823529,0.564705882352941,1;0,0.749019607843137,1;0.627450980392157,0.823529411764706,1;0.823529411764706,0.960784313725490,1;1,1,0.784313725490196;1,0.882352941176471,0.196078431372549;1,0.666666666666667,0;1,0.431372549019608,0;1,0,0;0.784313725490196,0,0;0.627450980392157,0.137254901960784,0.137254901960784];
new_jet=[0,0,128;0,0,133;0,0,138;0,0,143;0,0,148;0,0,153;0,0,158;0,0,164;0,0,169;0,0,174;0,0,179;0,0,184;0,0,189;0,0,194;0,0,199;0,0,204;0,0,209;0,0,214;0,0,219;0,0,225;0,0,230;0,0,235;0,0,240;0,0,245;0,0,250;0,0,255;10,10,255;20,20,255;31,31,255;41,41,255;51,51,255;61,61,255;71,71,255;82,82,255;92,92,255;102,102,255;112,112,255;122,122,255;133,133,255;143,143,255;153,153,255;163,163,255;173,173,255;184,184,255;194,194,255;204,204,255;214,214,255;224,224,255;235,235,255;245,245,255;255,255,255;255,245,245;255,235,235;255,224,224;255,214,214;255,204,204;255,194,194;255,184,184;255,173,173;255,163,163;255,153,153;255,143,143;255,133,133;255,122,122;255,112,112;255,102,102;255,92,92;255,82,82;255,71,71;255,61,61;255,51,51;255,41,41;255,31,31;255,20,20;255,10,10;255,0,0;250,0,0;245,0,0;240,0,0;235,0,0;230,0,0;225,0,0;219,0,0;214,0,0;209,0,0;204,0,0;199,0,0;194,0,0;189,0,0;184,0,0;179,0,0;174,0,0;169,0,0;164,0,0;158,0,0;153,0,0;148,0,0;143,0,0;138,0,0;133,0,0;128,0,0]./255;

lon=ncread('../R115ka_dkrz/R115ka_dkrz_echam6_echam_3900-3999.ymonmeanfixcalendar.mermean.nc','lon');
lat=ncread('../R115ka_dkrz/R115ka_dkrz_echam6_echam_3900-3999.ymonmean.inso.zonmean.nc','lat');
[yy,tt]=meshgrid(lat,1:12);

slm=ncread('../Rruns/slm.nc','slm');

[yy2,tt2]=meshgrid(lat,1:366);

titlename={'115ka','pre-industrial','115ka-pi'};
expname={'Rpi','R115ka'};
calendar_conv={'yes','no',''};
varname={'srad0d','temp2','seaice','tsurf','siced'};
mask={'land','ocn','globe'};
region={'Arctic','Antarc'};
region2={'Arctic','Antarctic'};

cnlevels=[-10:1:1,-0.3:0.1:0.3,1:1:10];

hold on

% inso
cnlevels=cnlevels*5;
for vv=1:1

    subplot(3,2,1)
    filename{1}=['../Rruns/calendar_conversion/',expname{1},'_',varname{vv},'_no.nc'];
    draw_value=ncread(filename{1},varname{vv});
    draw(:,:,1) = squeeze(nanmean(draw_value,1));
    filename{2}=['../Rruns/calendar_conversion/',expname{2},'_',varname{vv},'_no.nc'];
    draw_value=ncread(filename{2},varname{vv});
    draw(:,:,2) = squeeze(nanmean(draw_value,1));
    [C,h]=contourf(tt',yy',draw(:,:,2)-draw(:,:,1),cnlevels,'LineStyle','none');
    h.Parent.XTick=1:1:12;
    h.Parent.XTickLabel={'J','F','M','A','M','J','J','A','S','O','N','D'};
    h.Parent.YTick=-60:30:60;
    h.Parent.YTickLabel={'60S','30S','0','30N','60N'};
    ylim([0 88.5])
    colormap(new_jet)
    h0=colorbar;
    h0.Ticks=[-8:4:8]*5;
    title('(a) insolation, classical')
    caxis([-55 55])
    set(gca,'FontSize',15)

    clear draw_value draw
    subplot(3,2,2)
    filename{1}=['../Rruns/calendar_conversion/',expname{1},'_',varname{vv},'_yes.nc'];
    draw_value=ncread(filename{1},varname{vv});
    draw(:,:,1) = squeeze(nanmean(draw_value,1));
    filename{2}=['../Rruns/calendar_conversion/',expname{2},'_',varname{vv},'_yes.nc'];
    draw_value=ncread(filename{2},varname{vv});
    draw(:,:,2) = squeeze(nanmean(draw_value,1));
    [C,h]=contourf(tt',yy',draw(:,:,2)-draw(:,:,1),cnlevels,'LineStyle','none');
    h.Parent.XTick=1:1:12;
    h.Parent.XTickLabel={'J','F','M','A','M','J','J','A','S','O','N','D'};
    h.Parent.YTick=-60:30:60;
    h.Parent.YTickLabel={'60S','30S','0','30N','60N'};
    ylim([0 88.5])
    colormap(new_jet)
    h0=colorbar;
    h0.Ticks=[-8:4:8]*5;
    title('(b) insolation, angular')
    caxis([-55 55])
    set(gca,'FontSize',15)


end

% temp2 over land
for vv=2:2

    subplot(3,2,3)
    filename{1}=['../Rruns/calendar_conversion/',expname{1},'_',varname{vv},'_no.nc'];
    draw_value=ncread(filename{1},varname{vv});
    % land
    draw_value=draw_value.*slm;
    draw(:,:,1) = squeeze(nanmean(draw_value,1));
    filename{2}=['../Rruns/calendar_conversion/',expname{2},'_',varname{vv},'_no.nc'];
    draw_value=ncread(filename{2},varname{vv});
    % land
    draw_value=draw_value.*slm;
    draw(:,:,2) = squeeze(nanmean(draw_value,1));
    [C,h]=contourf(tt',yy',draw(:,:,2)-draw(:,:,1),20,'LineStyle','none');
    h.Parent.XTick=1:1:12;
    h.Parent.XTickLabel={'J','F','M','A','M','J','J','A','S','O','N','D'};
    h.Parent.YTick=-60:30:60;
    h.Parent.YTickLabel={'60S','30S','0','30N','60N'};
    ylim([0 88.5])
    %         colormap(new_jet)
    h0=colorbar;
    h0.Ticks=-4:2:4;
    title('(c) ground temperature, classical')
    caxis([-6 6])
    set(gca,'FontSize',15)

    clear draw_value draw
    subplot(3,2,4)
    filename{1}=['../Rruns/calendar_conversion/',expname{1},'_',varname{vv},'_yes.nc'];
    draw_value=ncread(filename{1},varname{vv});
    % land
    draw_value=draw_value.*slm;
    draw(:,:,1) = squeeze(nanmean(draw_value,1));
    filename{2}=['../Rruns/calendar_conversion/',expname{2},'_',varname{vv},'_yes.nc'];
    draw_value=ncread(filename{2},varname{vv});
    % land
    draw_value=draw_value.*slm;
    draw(:,:,2) = squeeze(nanmean(draw_value,1));
    [C,h]=contourf(tt',yy',draw(:,:,2)-draw(:,:,1),20,'LineStyle','none');
    h.Parent.XTick=1:1:12;
    h.Parent.XTickLabel={'J','F','M','A','M','J','J','A','S','O','N','D'};
    h.Parent.YTick=-60:30:60;
    h.Parent.YTickLabel={'60S','30S','0','30N','60N'};
    ylim([0 88.5])
    colormap(new_jet)
    h0=colorbar;
    h0.Ticks=-4:2:4;
    title('(d) ground temperature, angular')
    caxis([-6 6])
    set(gca,'FontSize',15)


end

% temp2 over sea
% cnlevels=[-8:2:1,-0.6:0.2:0.6,1:2:8];
% cnlevels=-9:.3:9;
cnlevels=[-8:1:1,-0.9:0.1:0.9,1:1:8];
for vv=2:2

    subplot(3,2,5)
    filename{1}=['../Rruns/calendar_conversion/',expname{1},'_',varname{vv},'_no.nc'];
    draw_value=ncread(filename{1},varname{vv});
    % ocn
    draw_value=draw_value.*(1-slm);
    draw(:,:,1) = squeeze(nanmean(draw_value,1));
    filename{2}=['../Rruns/calendar_conversion/',expname{2},'_',varname{vv},'_no.nc'];
    draw_value=ncread(filename{2},varname{vv});
    % ocn
    draw_value=draw_value.*(1-slm);
    draw(:,:,2) = squeeze(nanmean(draw_value,1));
    [C,h]=contourf(tt',yy',draw(:,:,2)-draw(:,:,1),cnlevels,'LineStyle','none');
    h.Parent.XTick=1:1:12;
    h.Parent.XTickLabel={'J','F','M','A','M','J','J','A','S','O','N','D'};
    h.Parent.YTick=-60:30:60;
    h.Parent.YTickLabel={'60S','30S','0','30N','60N'};
    ylim([0 88.5])
    %         colormap(new_jet)
    h0=colorbar;
    h0.Ticks=-4:2:4;
    title('(e) ocean temperature, classical')
    caxis([-6 6])
    set(gca,'FontSize',15)

    clear draw_value draw
    subplot(3,2,6)
    filename{1}=['../Rruns/calendar_conversion/',expname{1},'_',varname{vv},'_yes.nc'];
    draw_value=ncread(filename{1},varname{vv});
    % ocn
    draw_value=draw_value.*(1-slm);
    draw(:,:,1) = squeeze(nanmean(draw_value,1));
    filename{2}=['../Rruns/calendar_conversion/',expname{2},'_',varname{vv},'_yes.nc'];
    draw_value=ncread(filename{2},varname{vv});
    % ocn
    draw_value=draw_value.*(1-slm);
    draw(:,:,2) = squeeze(nanmean(draw_value,1));
    [C,h]=contourf(tt',yy',draw(:,:,2)-draw(:,:,1),cnlevels,'LineStyle','none');
    h.Parent.XTick=1:1:12;
    h.Parent.XTickLabel={'J','F','M','A','M','J','J','A','S','O','N','D'};
    h.Parent.YTick=-60:30:60;
    h.Parent.YTickLabel={'60S','30S','0','30N','60N'};
    ylim([0 88.5])
    colormap(new_jet)
    h0=colorbar;
    h0.Ticks=-4:2:4;
    title('(f) ocean temperature, angular')
    caxis([-6 6])
    set(gca,'FontSize',15)

end

set(gcf,'Position',[360 64 646 633])
print('-dpng',['A11_inso_vs_temp_land_sea','.png'])
% close all
