clear all;format compact

% 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

f2 = '../seeding_masks/b.e21.BSSP245cmip6.f09_g17.CMIP6-MCB-cntl.000.SWCF.nc';
f1 = '../seeding_masks/b.e21.BSSP245cmip6.f09_g17.CMIP6-baseline.000.SWCF.nc';

tmp2 = ncread(f2,'SWCF');
tmp1 = ncread(f1,'SWCF');

nn = 240;

SWCF1 = tmp1(:,:,1:nn);
SWCF2 = tmp2(:,:,1:nn);

swcf = SWCF2-SWCF1;

clear tmp1 tmp2 SWCF1 SWCF2

tmp2 = ncread(f2,'FSNTOA');
tmp1 = ncread(f1,'FSNTOA');

FSNTOA1 = tmp1(:,:,1:nn);
FSNTOA2 = tmp2(:,:,1:nn);

fsntoa = FSNTOA2-FSNTOA1;

clear tmp1 tmp2 FSNTOA1 FSNTOA2

tmp2 = ncread(f2,'FSNTOAC');
tmp1 = ncread(f1,'FSNTOAC');

FSNTOAC1 = tmp1(:,:,1:nn);
FSNTOAC2 = tmp2(:,:,1:nn);

fsntoac = FSNTOAC2-FSNTOAC1;

clear tmp1 tmp2 FSNTOAC1 FSNTOAC2

tmp2 = ncread(f2,'FLUTC');
tmp1 = ncread(f1,'FLUTC');

FLUTC1 = tmp1(:,:,1:nn);
FLUTC2 = tmp2(:,:,1:nn);

flutc = FLUTC2-FLUTC1;

clear tmp1 tmp2 FLUTC1 FLUTC2

tmp2 = ncread(f2,'FLUT');
tmp1 = ncread(f1,'FLUT');

FLUT1 = tmp1(:,:,1:nn);
FLUT2 = tmp2(:,:,1:nn);

flut = FLUT2-FLUT1;

clear tmp1 tmp2 FLUT1 FLUT2



tmp = ncread(f2,'OCNFRAC');
ocnfrac = tmp(:,:,1:nn);
clear tmp

lon = ncread(f1,'lon');im=length(lon);
lat = ncread(f1,'lat');jm=length(lat);

gw = ncread(f2,'gw');

%r = [5:5:50]/100;
%r=[.025 .05 .075 .1 .125 .15 .175 .2];
%r = [.05:.05:1];
r = linspace(0,1,101);

for i=1:nn
    year = floor((i-0.1)/12);
    mon = i-year*12;
    year = year + 1;

    SWCF(:,:,mon,year) = squeeze(swcf(:,:,i));
    FSNTOA(:,:,mon,year) = squeeze(fsntoa(:,:,i));
    FSNTOAC(:,:,mon,year) = squeeze(fsntoac(:,:,i));
    FLUTC(:,:,mon,year) = squeeze(flutc(:,:,i));
    FLUT(:,:,mon,year) = squeeze(flut(:,:,i));
    OCNFRAC(:,:,mon,year) = squeeze(ocnfrac(:,:,i));
end

mask = zeros(im,jm,12,length(r));


forcing = zeros(12,length(r));
f1 = zeros(12,length(r));
f2 = zeros(12,length(r));
f3 = zeros(12,length(r));
f4 = zeros(12,length(r));

for mon = 1:12

    total_area = 0;
    pts = 0;
    for i=1:im;for j=1:jm
        if( mean(OCNFRAC(i,j,mon,:),4)>=0.95 ) 
         total_area = total_area+gw(j);
         pts = pts+1;
         S(pts,1) = squeeze(mean(SWCF(i,j,mon,:),4));
         S(pts,2) = i;
         S(pts,3) = j;
         F1(pts,1) = squeeze(mean(FSNTOA(i,j,mon,:),4));
         F2(pts,1) = squeeze(mean(FSNTOAC(i,j,mon,:),4));
         F3(pts,1) = squeeze(mean(FLUTC(i,j,mon,:),4));
         F4(pts,1) = squeeze(mean(FLUT(i,j,mon,:),4));
        end
    end;end   

    [SS,II]=sort(S(:,1),'ascend');
    area = 0;
    for i=1:pts
      area = area+gw(S(II(i),3));
      for n=1:length(r)
       if( area<=r(n)*total_area ) 
        mask(S(II(i),2),S(II(i),3),mon,n) = 1;

        forcing(mon,n) = forcing(mon,n)+gw(S(II(i),3))*S(II(i),1)/2/im;

        f1(mon,n) = f1(mon,n)+gw(S(II(i),3))*F1(II(i),1)/2/im;
        f2(mon,n) = f2(mon,n)+gw(S(II(i),3))*F2(II(i),1)/2/im;
        f3(mon,n) = f3(mon,n)+gw(S(II(i),3))*F3(II(i),1)/2/im;
        f4(mon,n) = f4(mon,n)+gw(S(II(i),3))*F4(II(i),1)/2/im;
       end
      end
    end 

    clear S SS II
end

figure;orient tall;orient landscape
subplot('position',[.1 .1 .35 .4])
ANN = mean(forcing,1);
set(plot(r*100,ANN,'k-'),'linewidth',1)
xlabel('percentage of ocean area')
ylabel('\DeltaSWCF (W/m2)')
title('(a) Radiative forcing of MCB-375')
axis([0 100 -18 0])

FF1 = mean(f1,1);
FF2 = mean(f2,1);
FF3 = mean(f3,1);
FF4 = mean(f4,1);

subplot('position',[.5 .1 .35 .4])
set(plot(r*100,FF1,'r-'),'linewidth',1)
hold on
set(plot(r*100,FF2,'g-'),'linewidth',1)
set(plot(r*100,FF3,'b-'),'linewidth',1)
set(plot(r*100,FF4,'m-'),'linewidth',1)
xlabel('percentage of ocean area')
ylabel('radiative fluxes (W/m2)')
title('(b) Radiative fluxes at TOA')
legend('FSNTOA','FSNTOAC','FLUT','FLUTC')
legend('boxoff')
axis([0 100 -18 0])