format compact;clear all


cases = {'b.e21.BSSP245cmip6.f09_g17.CMIP6-baseline.000';...
         'b.e21.BSSP245cmip6.f09_g17.CMIP6-MCB-025PCT.000';...
         'b.e21.BSSP245cmip6.f09_g17.CMIP6-MCB-050PCT.000';...
         'b.e21.BSSP245cmip6.f09_g17.CMIP6-MCB-075PCT.000';...
         'b.e21.BSSP245cmip6.f09_g17.CMIP6-MCB-125PCT.000'};
         
vars = {'LANDFRAC';'TREFHT'};

dates = {'201501-206912';...
         '203501-206912';...
         '203501-206912';...
         '203501-206912';...
         '203501-206912'};
         

data = struct('casename',{},'varname',{},'lon',{},'lat',{},'date',{},'year',{},'month',{},'gw',{},'value',{},'global_avg',{},'annual_avg',{},'annual_avg_year',{});


for i=1:length(cases);
for j=1:length(vars)+2;
    
    if( j<=length(vars) )
      f = ['../' char(cases(i)) '/atm/proc/tseries/month_1/' char(cases(i)) '.cam.h0.' char(vars(j)) '.' char(dates(i)) '.nc']
    
    
      data(i,j).casename = char(cases(i));
      data(i,j).varname = char(vars(j));
      data(i,j).lon = ncread(f,'lon');
      data(i,j).lat = ncread(f,'lat');
      data(i,j).date = ncread(f,'date');
      data(i,j).gw = ncread(f,'gw');
      data(i,j).value = ncread(f,char(vars(j)));

    else
      data(i,j).casename = data(i,j-1).casename;
      data(i,j).lon = data(i,1).lon;
      data(i,j).lat = data(i,1).lat;
      data(i,j).date = data(i,1).date;
      data(i,j).gw = data(i,1).gw;   
      data(i,j).value = data(i,2).value;
      if( j==length(vars)+1 )
        for jj=1:length(data(1,1).lat)
		 data(i,j).value(:,jj,:) = data(i,j).value(:,jj,:).*sin(data(1,1).lat(jj)/180*pi);
        end
      else
        for jj=1:length(data(1,1).lat)
		 data(i,j).value(:,jj,:) = data(i,j).value(:,jj,:).*0.5.*(3*sin(data(1,1).lat(jj)/180*pi).*sin(data(1,1).lat(jj)/180*pi)-1);
        end 
      end
    end
    [data(i,j).year data(i,j).month] = set_yr_mon(data(i,j).date);

    for n=1:length(data(i,j).date)
        data(i,j).global_avg(:,n) =  global_avg(squeeze(data(i,j).value(:,:,n)),data(i,j).gw,squeeze(data(i,1).value(:,:,n)));
    end 
    
    for m=1:3
       [data(i,j).annual_avg(m,:) data(i,j).annual_avg_year] = annual_avg_1d(squeeze(data(i,j).global_avg(m,:)),data(i,j).year,data(i,j).month);
    end
end
end


figure(1);orient tall;orient landscape
j=2;
m=1;
subplot('position',[.1 .7 .7 .25])
hold on

l1=plot([2015:1:2069],data(1,j).annual_avg(m,:),'k-');
l2=plot([2035:1:2069],data(2,j).annual_avg(m,:),'g-');
l3=plot([2035:1:2069],data(3,j).annual_avg(m,:),'r-');
l4=plot([2035:1:2069],data(4,j).annual_avg(m,:),'b-');
l5=plot([2035:1:2069],data(5,j).annual_avg(m,:),'m-');

set(l1,'linewidth',1)
set(l2,'linewidth',1)
set(l3,'linewidth',1)
set(l4,'linewidth',1)
set(l5,'linewidth',1)

set(title('(a) Average surface temperature (K) global'),'fontsize',10)
set(gca,'xlim',[2020 2070])
plot([2020 2039],[mean(data(1,j).annual_avg(m,6:25)) mean(data(1,j).annual_avg(m,6:25))],'k--')

bb = mean(data(2,j).annual_avg(m,16:35));
cc = mean(data(3,j).annual_avg(m,16:35));
dd = mean(data(4,j).annual_avg(m,16:35));
ee = mean(data(5,j).annual_avg(m,16:35));

plot([2050 2069],[bb bb],'g--')
plot([2050 2069],[cc cc],'r--')
plot([2050 2069],[dd dd],'b--')
plot([2050 2069],[ee ee],'m--')

clear a aa bb cc dd ee
box on
set(gca,'fontsize',8)
ylabel('Temperature (K)')

legend([l1 l2 l3 l4 l5],{'CESM SSP2.4-5','MCB 2.5%','MCB 5%','MCB 7.5%','MCB 12.5%'},'Location','Northwest')
legend('boxoff')

m=2;
subplot('position',[.1 .4 .7 .25])
hold on

l1=plot([2015:1:2069],data(1,j).annual_avg(m,:),'k-');
l2=plot([2035:1:2069],data(2,j).annual_avg(m,:),'g-');
l3=plot([2035:1:2069],data(3,j).annual_avg(m,:),'r-');
l4=plot([2035:1:2069],data(4,j).annual_avg(m,:),'b-');
l5=plot([2035:1:2069],data(5,j).annual_avg(m,:),'m-');

set(l1,'linewidth',1)
set(l2,'linewidth',1)
set(l3,'linewidth',1)
set(l4,'linewidth',1)
set(l5,'linewidth',1)

set(title('(b) Average surface temperature (K) over land'),'fontsize',10)
set(gca,'xlim',[2020 2070])
plot([2020 2039],[mean(data(1,j).annual_avg(m,6:25)) mean(data(1,j).annual_avg(m,6:25))],'k--')

bb = mean(data(2,j).annual_avg(m,16:35));
cc = mean(data(3,j).annual_avg(m,16:35));
dd = mean(data(4,j).annual_avg(m,16:35));
ee = mean(data(5,j).annual_avg(m,16:35));

plot([2050 2069],[bb bb],'g--')
plot([2050 2069],[cc cc],'r--')
plot([2050 2069],[dd dd],'b--')
plot([2050 2069],[ee ee],'m--')

clear a aa bb cc dd ee
box on
set(gca,'fontsize',8)
ylabel('Temperature (K)')

m=3;
subplot('position',[.1 .1 .7 .25])
hold on

l1=plot([2015:1:2069],data(1,j).annual_avg(m,:),'k-');
l2=plot([2035:1:2069],data(2,j).annual_avg(m,:),'g-');
l3=plot([2035:1:2069],data(3,j).annual_avg(m,:),'r-');
l4=plot([2035:1:2069],data(4,j).annual_avg(m,:),'b-');
l5=plot([2035:1:2069],data(5,j).annual_avg(m,:),'m-');

set(l1,'linewidth',1)
set(l2,'linewidth',1)
set(l3,'linewidth',1)
set(l4,'linewidth',1)
set(l5,'linewidth',1)

set(title('(c) Average surface temperature (K) over ocean'),'fontsize',10)
set(gca,'xlim',[2020 2070])
plot([2020 2039],[mean(data(1,j).annual_avg(m,6:25)) mean(data(1,j).annual_avg(m,6:25))],'k--')

bb = mean(data(2,j).annual_avg(m,16:35));
cc = mean(data(3,j).annual_avg(m,16:35));
dd = mean(data(4,j).annual_avg(m,16:35));
ee = mean(data(5,j).annual_avg(m,16:35));

plot([2050 2069],[bb bb],'g--')
plot([2050 2069],[cc cc],'r--')
plot([2050 2069],[dd dd],'b--')
plot([2050 2069],[ee ee],'m--')

clear a aa bb cc dd ee
box on
set(gca,'fontsize',8)
ylabel('Temperature (K)')
xlabel('Time (year)')