format compact;clear all

en_number = 10;

cases = {'b.e21.BSSP245smbb.f09_g17.001';...
         'b.e21.BSSP245smbb.f09_g17.002';...
         'b.e21.BSSP245smbb.f09_g17.003';...
         'b.e21.BSSP245smbb.f09_g17.004';...
         'b.e21.BSSP245smbb.f09_g17.005';...
         'b.e21.BSSP245smbb.f09_g17.006';...
         'b.e21.BSSP245smbb.f09_g17.007';...
         'b.e21.BSSP245smbb.f09_g17.009';...
         'b.e21.BSSP245smbb.f09_g17.010';...
         'b.e21.BSSP245smbb.f09_g17.011';... 
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.001';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.002';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.003';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.004';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.005';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.006';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.007';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.009';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.010';...
         'b.e21.BSSP245smbb.f09_g17.MCB-050PCT.011'};

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

var = 'TREFHT';

vars = {'LANDFRAC';var};

for i=1:length(cases);
for j=1:length(vars);

    
    if( i<=10 )
      f = ['../../' char(cases(i)) '/atm/proc/tseries/month_1/' char(cases(i)) '.cam.h0.' char(vars(j)) '.201501-206412.nc']
	f2 = ['../../' char(cases(i)) '/atm/proc/tseries/month_1/' char(cases(i)) '.cam.h0.' char(vars(j)) '.206501-210012.nc'];
      date1 = ncread(f,'date');
      date2 = ncread(f2,'date');
      data(i,j).date = [ncread(f,'date');ncread(f2,'date')];
    %if(j>=length(vars)-1)
    %  xxx = ncread(f,char(vars(j)));
    %  data(i,j).value = squeeze(xxx(:,:,32,:));
    %  clear xxx
    %else
      data(i,j).value = ncread(f,char(vars(j)));
    %end
      data(i,j).value(:,:,length(date1)+1:length(date1)+length(date2)) = ncread(f2,char(vars(j)));
      clear date1 date2
    else 
      f = ['../' char(cases(i)) '/atm/proc/tseries/month_1/' char(cases(i)) '.cam.h0.' char(vars(j)) '.' char(dates(i)) '.nc']
      data(i,j).date = ncread(f,'date');
    %if(j>=length(vars)-1)
    %  xxx = ncread(f,char(vars(j)));
    %  data(i,j).value = squeeze(xxx(:,:,32,:));
    %  clear xxx
    %else
      data(i,j).value = ncread(f,char(vars(j)));
    %end
    end

    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).gw = ncread(f,'gw');
    [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

j=2;
 
for i=1:10
    x1(i,:,:,:) = data(i,j).value(:,:,61:300);
    x2(i,:,:,:) = data(i,j).value(:,:,421:660);
end

for i=11:20
    x3(i-10,:,:,:) = data(i,j).value(:,:,181:420);
end

X1 = mean(x1,4);
X2 = mean(x2,4);
X3 = mean(x3,4);

x1_JJA = zeros(10,length(data(1,1).lon),length(data(1,1).lat),60);
x2_JJA = x1_JJA;
x3_JJA = x1_JJA;
x1_DJF = x1_JJA;
x2_DJF = x1_JJA;
x3_DJF = x1_JJA;

count = 0;   
for year = 1:20
count = count+1;
x1_JJA(:,:,:,count) = x1(:,:,:,(year-1)*12+6);
x1_DJF(:,:,:,count) = x1(:,:,:,(year-1)*12+1);
x2_JJA(:,:,:,count) = x2(:,:,:,(year-1)*12+6);
x2_DJF(:,:,:,count) = x2(:,:,:,(year-1)*12+1);
x3_JJA(:,:,:,count) = x3(:,:,:,(year-1)*12+6);
x3_DJF(:,:,:,count) = x3(:,:,:,(year-1)*12+1);
count = count+1;
x1_JJA(:,:,:,count) = x1(:,:,:,(year-1)*12+7);
x1_DJF(:,:,:,count) = x1(:,:,:,(year-1)*12+2);
x2_JJA(:,:,:,count) = x2(:,:,:,(year-1)*12+7);
x2_DJF(:,:,:,count) = x2(:,:,:,(year-1)*12+2);
x3_JJA(:,:,:,count) = x3(:,:,:,(year-1)*12+7);
x3_DJF(:,:,:,count) = x3(:,:,:,(year-1)*12+2);
count = count+1;
x1_JJA(:,:,:,count) = x1(:,:,:,(year-1)*12+8);
x1_DJF(:,:,:,count) = x1(:,:,:,(year-1)*12+12);
x2_JJA(:,:,:,count) = x2(:,:,:,(year-1)*12+8);
x2_DJF(:,:,:,count) = x2(:,:,:,(year-1)*12+12);
x3_JJA(:,:,:,count) = x3(:,:,:,(year-1)*12+8);
x3_DJF(:,:,:,count) = x3(:,:,:,(year-1)*12+12);
end

X1_JJA = mean(x1_JJA,4);
X2_JJA = mean(x2_JJA,4);
X3_JJA = mean(x3_JJA,4);

X1_DJF = mean(x1_DJF,4);
X2_DJF = mean(x2_DJF,4);
X3_DJF = mean(x3_DJF,4);

% coast
load coastlines;
lon2=coastlon;
lat2=coastlat;

lon1=coastlon;
lat1=coastlat;
for i=1:length(lon1)
  if(lon1(i)<0)
    lon1(i)=lon1(i)+360;
  end
  if((lon1(i)<=1)||(lon1(i)>=359))
    lon1(i)=NaN;
  end
end
clear long lat lon2 lat2
% coast

lon = data(1,1).lon;
lat = data(1,1).lat;

im = length(lon);
jm = length(lat);

level = 41;n = ceil(level/2);
cmap1 = [linspace(0, 1, n); linspace(0, 1, n); linspace(1, 1, n)]';
cmap2 = [linspace(1, 1, n); linspace(1, 0, n); linspace(1, 0, n)]';
ctb = [cmap1; cmap2(2:end, :)];

save([char(var) '.mat'],'-v7.3')