clear all;format compact

load S5

XX = zeros(length(cases),length(vars)-4,192,55);

for i=1:length(cases);

for j=5:length(vars)

    if(i<=10)
      tmp = squeeze(mean(data(i,j).value(:,:,1:660),1));
      for n=1:55
            XX(i,j-4,:,n) = squeeze(mean(tmp(:,(n-1)*12+1:n*12),2));
      end
    else
      tmp = squeeze(mean(data(i,j).value,1));
      for n=21:55
              XX(i,j-4,:,n) = squeeze(mean(tmp(:,(n-21)*12+1:(n-20)*12),2));
      end
    end


    clear tmp
end
end

for j=1:length(vars)-4
X1(1:10,j,:) = squeeze(mean(XX(1:10,j,:,6:25),4));
X2(1:10,j,:) = squeeze(mean(XX(1:10,j,:,36:55),4));
X3(1:10,j,:) = squeeze(mean(XX(11:20,j,:,36:55),4));
end


panels = [0.07 .7 .12 .25;.21 .7 .12 .25;.35 .7 .12 .25;...
          0.49 .7 .12 .25;0.63 .7 .12 .25;.77 .7 .12 .25;...
          0.07 .4 .12 .25;.21 .4 .12 .25;.35 .4 .12 .25;...
          0.49 .4 .12 .25;0.63 .4 .12 .25;.77 .4 .12 .25;...
          0.07 .1 .12 .25;.21 .1 .12 .25;.35 .1 .12 .25;...
          0.49 .1 .12 .25;0.63 .1 .12 .25;.77 .1 .12 .25];

font = 6;

figure;orient tall;orient landscape;
subplot('position',panels(1,:))
aa = squeeze(X3(:,2,:)-X1(:,2,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,2,:)-X1(:,2,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
set(ylabel('latitude'),'fontsize',font)
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(a) \DeltaAllsky LW flux at model top'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(2,:))
aa = squeeze(X3(:,4,:)-X1(:,4,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,4,:)-X1(:,4,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(b) \DeltaClearsky LW flux at model top'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(3,:))
aa = -squeeze(X3(:,6,:)-X1(:,6,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(-squeeze(mean(X3(:,6,:)-X1(:,6,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(c) \DeltaLW cloud forcing'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(4,:))
aa = squeeze(X3(:,3,:)-X1(:,3,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,3,:)-X1(:,3,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(d) \DeltaAllsky SW flux at model top'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(5,:))
aa = squeeze(X3(:,5,:)-X1(:,5,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,5,:)-X1(:,5,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(e) \DeltaClearsky SW flux at model top'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(6,:))
aa = squeeze(X3(:,7,:)-X1(:,7,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,7,:)-X1(:,7,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(f) \DeltaSW cloud forcing'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(7,:))
aa = squeeze(X3(:,20,:)+2*X3(:,21,:)-X1(:,20,:)-2*X1(:,21,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,20,:)+2*X3(:,21,:)-X1(:,20,:)-2*X1(:,21,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
set(ylabel('latitude'),'fontsize',font)
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(g) \DeltaTotal LW flux at surface'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(8,:))
aa = squeeze(X3(:,20,:)+X3(:,21,:)-X1(:,20,:)-X1(:,21,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,20,:)+X3(:,21,:)-X1(:,20,:)-X1(:,21,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(h) \DeltaUpward LW flux at surface'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(9,:))
aa = squeeze(X3(:,21,:)-X1(:,21,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,21,:)-X1(:,21,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(i) \DeltaDownward LW flux at surface'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(10,:))
aa = squeeze(2*X3(:,15,:)-X3(:,14,:)-2*X1(:,15,:)+X1(:,14,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(2*X3(:,15,:)-X3(:,14,:)-2*X1(:,15,:)+X1(:,14,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(j) \DeltaTotal SW flux at surface'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(11,:))
aa = squeeze(X3(:,15,:)-X3(:,14,:)-X1(:,15,:)+X1(:,14,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,15,:)-X3(:,14,:)-X1(:,15,:)+X1(:,14,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(k) \DeltaUpward SW flux at surface'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(12,:))
aa = squeeze(X3(:,15,:)-X1(:,15,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,15,:)-X1(:,15,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('radiative forcing (W/m2)'),'fontsize',font)
set(title('(l) \DeltaDownward SW flux at surface'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(13,:))
aa = squeeze(X3(:,10,:)-X1(:,10,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,10,:)-X1(:,10,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
set(ylabel('latitude'),'fontsize',font)
%set(gca,'yticklabel',[])
set(xlabel('cloud fraction(%)'),'fontsize',font)
set(title('(m) \DeltaLow cloud fraction'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(14,:))
aa = squeeze(X3(:,9,:)-X1(:,9,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,9,:)-X1(:,9,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('cloud fraction(%)'),'fontsize',font)
set(title('(n) \DeltaMid cloud fraction'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(15,:))
aa = squeeze(X3(:,8,:)-X1(:,8,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,8,:)-X1(:,8,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('cloud fraction(%)'),'fontsize',font)
set(title('(o) \DeltaHigh cloud fraction'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(16,:))
aa = squeeze(X3(:,12,:)-X1(:,12,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh; hold on
l1 = plot(squeeze(mean(X3(:,12,:)-X1(:,12,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('cloud fraction(%)'),'fontsize',font)
set(title('(p) \DeltaTotal cloud fraction'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(17,:))
aa = squeeze(X3(:,11,:)-X1(:,11,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh aa; hold on
l1 = plot(squeeze(mean(X3(:,11,:)-X1(:,11,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('temperature (K)'),'fontsize',font)
set(title('(q) \DeltaSurface temperature'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)

subplot('position',panels(18,:))
aa = squeeze(X3(:,13,:)-X1(:,13,:));
yy1 = mean(aa,1)-2*std(aa,1);
yy2 = mean(aa,1)+2*std(aa,1);
hhh = patch([yy1 fliplr(yy2)],[data(1,1).lat' fliplr(data(1,1).lat')],[0 0 0.8]);
clear aa yy1 yy2 hhh aa; hold on
l1 = plot(squeeze(mean(X3(:,13,:)-X1(:,13,:),1)),data(1,1).lat,'b-');
%set(l1,'linewidth',1);
plot([0 0],[-90 90],'k--')
set(gca,'ylim',[-90 90],'ytick',[-90:30:90],'yminortick','on')
%set(ylabel('latitude'),'fontsize',font)
set(gca,'yticklabel',[])
set(xlabel('sea-ice fraction (%)'),'fontsize',font)
set(title('(r) \DeltaSea-ice fraction'),'fontsize',font)
set(gca,'fontsize',font)
box on
alpha(0.2)