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

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=[2.5:.5:20]/100;


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

    SWCF(:,:,mon,year) = squeeze(swcf(:,:,i));
    OCNFRAC(:,:,mon,year) = squeeze(ocnfrac(:,:,i));
end

mask = zeros(im,jm,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;
        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;
       end
      end
    end 

    clear S SS II
end

mon = 12;tt=[1:1:12];

    out_file = ['mask_0375_025-200_SSP245.nc'];
    ncid = netcdf.create(out_file,'NOCLOBBER');
    dimid1 = netcdf.defDim(ncid,'lon',im);
    dimid2 = netcdf.defDim(ncid,'lat',jm);
    dimid3 = netcdf.defDim(ncid,'time',mon);
    dimid4 = netcdf.defDim(ncid,'pct',length(r));
    varID3 = netcdf.defVar(ncid,'time','NC_DOUBLE',dimid3);
    varID1 = netcdf.defVar(ncid,'lon','NC_DOUBLE',dimid1);
    varID2 = netcdf.defVar(ncid,'lat','NC_DOUBLE',dimid2);
    varID4 = netcdf.defVar(ncid,'pct','NC_DOUBLE',dimid4);
    varID5 = netcdf.defVar(ncid,'mask','NC_DOUBLE',[dimid1 dimid2 dimid3 dimid4]);
    netcdf.endDef(ncid)
   
    netcdf.putVar(ncid,varID1,lon)
    netcdf.putVar(ncid,varID2,lat)
    netcdf.putVar(ncid,varID3,tt)
    netcdf.putVar(ncid,varID4,r)
    netcdf.putVar(ncid,varID5,mask);

    
  netcdf.close(ncid)