load ~/NC/grid
year=1:150;

SST_temp=ncread('~/Aerosol_Variaiton/IRUN/read_nc/ncfile/read_nc_TS.nc','TS');
OCNFRAC=ncread('/Users/yangyang/Aerosol_Variaiton/IRUN/read_nc/ncfile/cesm.nc','OCNFRAC');

SST_ANN=zeros(144,96,150);
SST_DJF=zeros(144,96,150);
SST_JJA=zeros(144,96,150);

SST_ANN(:,:,:)=mean(SST_temp,3);
SST_JJA(:,:,:)=mean(SST_temp(:,:,[6,7,8],:),3);
SST_DJF(:,:,1)=mean(SST_temp(:,:,[1,2],1),3);
for yearID=2:150
    SST_DJF(:,:,yearID)=(SST_temp(:,:,12,yearID-1)+SST_temp(:,:,1,yearID)+SST_temp(:,:,2,yearID))/3;
end

%Set fillValues
for i=1:144
    for j=1:96
        if OCNFRAC(i,j)<1
            SST_ANN(i,j,:)=-999;
            SST_DJF(i,j,:)=-999;
            SST_JJA(i,j,:)=-999;
        end
    end
end
save('SST.mat','SST_ANN','SST_DJF','SST_JJA')

nccreate('SST.nc','Year','Dimensions',{'time',Inf},'Format','classic')
ncwrite('SST.nc','Year', year)
ncwriteatt('SST.nc','Year','long_name','Year')
nccreate('SST.nc','lat','Dimensions',{'lat',96},'Format','classic')
ncwrite('SST.nc','lat', lat)
ncwriteatt('SST.nc','lat','units','degrees_north')
ncwriteatt('SST.nc','lat','long_name','Latitude')
nccreate('SST.nc','lon','Dimensions',{'lon',144},'Format','classic')
ncwrite('SST.nc','lon', lon)
ncwriteatt('SST.nc','lon','units','degrees_east')
ncwriteatt('SST.nc','lon','long_name','Longitude')
nccreate('SST.nc','SST_ANN','Dimensions',{'lon',144,'lat',96,'year',length(year)},'Format','classic')
ncwrite('SST.nc','SST_ANN', SST_ANN)
nccreate('SST.nc','SST_DJF','Dimensions',{'lon',144,'lat',96,'year',length(year)},'Format','classic')
ncwrite('SST.nc','SST_DJF', SST_DJF)
nccreate('SST.nc','SST_JJA','Dimensions',{'lon',144,'lat',96,'year',length(year)},'Format','classic')
ncwrite('SST.nc','SST_JJA', SST_JJA)
ncdisp('SST.nc')