clear all
filenames=dir('MERRA2_400.inst3_3d_aer_Nv.2011*')
def='./'

l=length(filenames)
k1=0;

ncdisp([def filenames(1).name])
ncid = netcdf.open([def filenames(1).name],'NOWRITE');

varid= netcdf.inqVarID(ncid,'lat'); 
LAT= netcdf.getVar(ncid,varid);
XLAT=(mean(LAT(:,:,1),3)) ;


varid= netcdf.inqVarID(ncid,'lon'); 
LONG= netcdf.getVar(ncid,varid);
XLONG=(mean(LONG(:,:,1),3)) ;

varid= netcdf.inqVarID(ncid,'lev'); 
lev= netcdf.getVar(ncid,varid);
level=(mean(lev(:,:,1),3)) ;


for k=1:l
  
ncid = netcdf.open([def filenames(k).name],'NOWRITE');
k
% 4D variables
%k1=k1+1;

% varid= netcdf.inqVarID(ncid,'DUSMASS'); 
% dustsurmass= netcdf.getVar(ncid,varid);
% Sur_DM=(mean(dustsurmass(:,:,1),3)) ;

varid= netcdf.inqVarID(ncid,'DU003'); 
DU003= netcdf.getVar(ncid,varid);
DU3(:,:,k)=(mean(DU003(25:35,:,:),1)) ;

varid= netcdf.inqVarID(ncid,'DU002'); 
DU002= netcdf.getVar(ncid,varid);
DU2(:,:,k)=(mean(DU002(25:35,:,:),1)) ;

varid= netcdf.inqVarID(ncid,'DU001'); 
DU001= netcdf.getVar(ncid,varid);
DU1(:,:,k)=(mean(DU001(25:35,:,:),1)) ;

varid= netcdf.inqVarID(ncid,'DU004'); 
DU004= netcdf.getVar(ncid,varid);
DU4(:,:,k)=(mean(DU004(25:35,:,:),1)) ;

varid= netcdf.inqVarID(ncid,'DU005'); 
DU005= netcdf.getVar(ncid,varid);
DU5(:,:,k)=(mean(DU005(25:35,:,:),1)) ;

varid= netcdf.inqVarID(ncid,'BCPHILIC'); 
BCPHILIC= netcdf.getVar(ncid,varid);
BC1(:,:,k)=(mean(BCPHILIC(25:35,:,:),1)) ;

varid= netcdf.inqVarID(ncid,'BCPHOBIC'); 
BCPHOBIC= netcdf.getVar(ncid,varid);
BC2(:,:,k)=(mean(BCPHOBIC(25:35,:,:),1)) ;

end


for k=1:l
  
ncid = netcdf.open([def filenames(k).name],'NOWRITE');
k
% 4D variables
%k1=k1+1;

% varid= netcdf.inqVarID(ncid,'DUSMASS'); 
% dustsurmass= netcdf.getVar(ncid,varid);
% Sur_DM=(mean(dustsurmass(:,:,1),3)) ;

varid= netcdf.inqVarID(ncid,'DU003'); 
DU003= netcdf.getVar(ncid,varid);
DU3_TP(:,:,k)=(mean(DU003(33:47,:,:),1)) ;

varid= netcdf.inqVarID(ncid,'DU002'); 
DU002= netcdf.getVar(ncid,varid);
DU2_TP(:,:,k)=(mean(DU002(33:47,:,:),1)) ;

varid= netcdf.inqVarID(ncid,'DU001'); 
DU001= netcdf.getVar(ncid,varid);
DU1_TP(:,:,k)=(mean(DU001(33:47,:,:),1)) ;

varid= netcdf.inqVarID(ncid,'DU004'); 
DU004= netcdf.getVar(ncid,varid);
DU4_TP(:,:,k)=(mean(DU004(33:47,:,:),1)) ;

varid= netcdf.inqVarID(ncid,'DU005'); 
DU005= netcdf.getVar(ncid,varid);
DU5_TP(:,:,k)=(mean(DU005(33:47,:,:),1)) ;


varid= netcdf.inqVarID(ncid,'BCPHILIC'); 
BCPHILIC= netcdf.getVar(ncid,varid);
BC1_TP(:,:,k)=(mean(BCPHILIC(33:47,:,:),1)) ;


varid= netcdf.inqVarID(ncid,'BCPHOBIC'); 
BCPHOBIC= netcdf.getVar(ncid,varid);
BC2_TP(:,:,k)=(mean(BCPHOBIC(33:47,:,:),1)) ;

end

%%
clear data
file=dir('/Volumes/Untitled/snowdata/HMAkarldata/2011/V6modscag_dalbedo_WRFlatlon*.mat')

names{1}=file(11).name;
names{2}=file(10).name;
names{3}=file(9).name;
names{4}=file(3).name;
names{5}=file(4).name;
names{6}=file(12).name;
names{7}=file(7).name;
names{8}=file(1).name;
names{9}=file(8).name;
names{10}=file(6).name;
names{11}=file(5).name;
names{12}=file(2).name;
def='F:\snowdata\HMAkarldata\2011\'

  mnth=[1 32 62 93 124 152 183 213 244 274 305 336 366]


data=zeros(210,150,l);

for j=6:10
data1=load([def names{j}]);
srt=mnth(j-1);
en=mnth(j);
if j==6
data(:,:,srt-mnth(5)+1:en-1-mnth(5)+1)=(data1.data1_hr);
else
data(:,:,srt-mnth(5)+1:en-1-mnth(5)+1)=(data1.data1_hr);
end
end



wrffilenames=dir('C:\Users\sara936\Desktop\HMA\evaluation\aerosol\wrfout_d01_2014-09-30_00%3A00%3A00')
dr='C:\Users\sara936\Desktop\HMA\evaluation\aerosol\'
%ncdisp('./wrf3hr_d01_2014-05-01_00%3A00%3A00')
ncid = netcdf.open([dr wrffilenames(1).name],'NOWRITE');

varid = netcdf.inqVarID(ncid,'XLONG'); 
lon_wrf = netcdf.getVar(ncid,varid);

varid = netcdf.inqVarID(ncid,'XLAT'); 
lat_wrf = netcdf.getVar(ncid,varid);



lon1=squeeze(lon_wrf(:,:,1));
lat1=squeeze(lat_wrf(:,:,1));

S=shaperead('C:\Users\sara936\Desktop\HMA\evaluation\paper1_analysis\00_rgi60_regions\00_rgi60_O2Regions.shp');

t=DU1+DU2+DU3+DU4+DU5;
t=t*1000000;

BC=BC1+BC2;
BC=BC*1000000000;

t_TP=DU1_TP+DU2_TP+DU3_TP+DU4_TP+DU5_TP;
t_TP=t_TP*1000000;

BC_TP=BC1_TP+BC2_TP;
BC_TP=BC_TP*1000000000;

for k=1:l

EALd(k)=(sum(mean(t(25:30,45:63,k),1),2));
EALb(k)=(sum(mean(BC(25:30,45:63,k),1),2));
EALd_TP(k)=(sum(mean(t_TP(20:24,35:63,k),1),2));
EALb_TP(k)=(sum(mean(BC_TP(19:24,40:62,k),1),2));
end

lon1=squeeze(lon_wrf(:,:,1));
lat1=squeeze(lat_wrf(:,:,1));

for k=1:l
deltsnowred=data(:,:,k);
in = inpolygon(lon1,lat1,S(62).X,S(62).Y);
DA_WH(k)=nanmean(nonzeros(deltsnowred(in)));

in = inpolygon(lon1,lat1,S(61).X,S(61).Y);
DA_KK(k)=nanmean(nonzeros(deltsnowred(in)));

in = inpolygon(lon1,lat1,S(60).X,S(60).Y);
DA_HK(k)=nanmean(nonzeros(deltsnowred(in)));

in = inpolygon(lon1,lat1,S(63).X,S(63).Y);
DA_CH(k)=nanmean(nonzeros(deltsnowred(in)));

in = inpolygon(lon1,lat1,S(64).X,S(64).Y);
DA_EH(k)=nanmean(nonzeros(deltsnowred(in)));

in = inpolygon(lon1,lat1,S(55).X,S(55).Y);
DA_TP(k)=nanmean(nonzeros(deltsnowred(in)));

end

save('2011data.mat','EALd','EALb','EALd_TP','EALb_TP','DA_WH','DA_EH','DA_HK','DA_KK','DA_TP')



%%
clear d *bin
d=prctile(EALd(:),[5:5:95]);
for i=1:length(d)-1
    id=EALd>d(i) & EALd<d(i+1)
    EAldbin(i)=nanmean(EALd(id));
     DA_TPbin(i)=nanmean(DA_TP(id));
     DA_WHbin(i)=nanmean(DA_WH(id));
     DA_EHbin(i)=nanmean(DA_EH(id));
     DA_HKbin(i)=nanmean(DA_HK(id));
     DA_KKbin(i)=nanmean(DA_KK(id));
     
     DA_CHbin(i)=nanmean(DA_CH(id));
      stdDA_TPbin(i)=nanstd(DA_TP(id));
     stdDA_WHbin(i)=nanstd(DA_WH(id));
     stdDA_EHbin(i)=nanstd(DA_EH(id));
     stdDA_HKbin(i)=nanstd(DA_HK(id));
     stdDA_KKbin(i)=nanstd(DA_KK(id));
     stdDA_CHbin(i)=nanstd(DA_CH(id));
end
%C=varycolor(4);


figure(1)
scatter(EALd(60:210),DA_WH(60:210),'MarkerFaceColor','r','MarkerEdgeColor','r',...
    'MarkerFaceAlpha',.3,'MarkerEdgeAlpha',.3)
hold on
errorbar(EAldbin,DA_WHbin,stdDA_WHbin,'rx')
hold on
[p2 S]= polyfit(EAldbin,DA_WHbin,1)
 [y delta] = polyval(p2,EAldbin,S);
hold on
plot(EAldbin,y,'--r','Linewidth',3)
box on
[px2 rmse] = corrcoef(EAldbin,DA_WHbin); 

  text(.05,.95,[ 'WH: R^2 = ',num2str(px2(1,2).^2,'%.2f') ],'Fontsize',15,'Units','normalized','color','r')

 %  text(.2,.95,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','r') 




hold on
scatter(EALd(60:210),DA_HK(60:210),'MarkerFaceColor','b','MarkerEdgeColor','b',...
    'MarkerFaceAlpha',.3,'MarkerEdgeAlpha',.3)
hold on
errorbar(EAldbin,DA_HKbin,stdDA_HKbin,'bx')
hold on
[p2 S]= polyfit(EAldbin,DA_HKbin,1)
 [y delta] = polyval(p2,EAldbin,S);
hold on
plot(EAldbin,y,'--b','Linewidth',3)
box on
[px2 rmse] = corrcoef(EAldbin,DA_HKbin); 
 text(.05,.9,[ 'HK: R^2 = ',num2str(px2(1,2).^2,'%.2f') ],'Fontsize',15,'Units','normalized','color','b')

  % text(.2,.9,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','b') 




hold on
scatter(EALd(60:210),DA_KK(60:210),'MarkerFaceColor','m','MarkerEdgeColor','m',...
    'MarkerFaceAlpha',.3,'MarkerEdgeAlpha',.3)
hold on
errorbar(EAldbin,DA_KKbin,stdDA_KKbin,'mx')
hold on
[p2 S]= polyfit(EAldbin,DA_KKbin,1)
 [y delta] = polyval(p2,EAldbin,S);
hold on
plot(EAldbin,y,'--m','Linewidth',3)
box on
[px2 rmse] = corrcoef(EAldbin,DA_KKbin); 
 text(.05,.85,[ 'KK: R^2 = ',num2str(px2(1,2).^2,'%.2f') ],'Fontsize',15,'Units','normalized','color','m')

  % text(.2,.85,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','m') 


%%



clear d *bin
d=prctile(EALd_TP(30:210),[0:4:100]);
for i=1:length(d)-1
    id=EALd_TP>d(i) & EALd_TP<d(i+1)
    EAldbin(i)=nanmean(EALd_TP(id));
     DA_TPbin(i)=nanmean(DA_TP(id));
     DA_WHbin(i)=nanmean(DA_WH(id));
     DA_EHbin(i)=nanmean(DA_EH(id));
     DA_HKbin(i)=nanmean(DA_HK(id));
     DA_KKbin(i)=nanmean(DA_KK(id));
     DA_CHbin(i)=nanmean(DA_CH(id));
      stdDA_TPbin(i)=nanstd(DA_TP(id));
     stdDA_WHbin(i)=nanstd(DA_WH(id));
     stdDA_EHbin(i)=nanstd(DA_EH(id));
     stdDA_HKbin(i)=nanstd(DA_HK(id));
     stdDA_KKbin(i)=nanstd(DA_KK(id));
     stdDA_CHbin(i)=nanstd(DA_CH(id));
end



hold on
scatter(EALd_TP(30:210),DA_TP(30:210),'MarkerFaceColor','k','MarkerEdgeColor','k',...
    'MarkerFaceAlpha',.3,'MarkerEdgeAlpha',.3)
hold on
errorbar(EAldbin,DA_TPbin,stdDA_TPbin,'kx')
hold on
[p2 S]= polyfit(EAldbin,DA_TPbin,1)
 [y delta] = polyval(p2,EAldbin,S);
hold on
plot(EAldbin,y,'--k','Linewidth',3)
box on
[px2 rmse] = corrcoef(EAldbin,DA_TPbin); 
 text(.05,.8,[ 'TP: R^2 = ',num2str(px2(1,2).^2,'%.2f') ],'Fontsize',15,'Units','normalized','color','k')

  % text(.2,.8,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','k') 

  xlabel('MERRA2-EAL_{dust}(mg/m^3)')
    ylabel('\Delta Albedo due to LAP (%)')
  set(gca,'fontsize',15)
% 
% hold on
% scatter(EALd_TP(30:210),(DA_CH(30:210)+DA_CH(30:210))/2,'MarkerFaceColor','g','MarkerEdgeColor','g',...
%     'MarkerFaceAlpha',.3,'MarkerEdgeAlpha',.3)
% hold on
% errorbar(EAldbin,(DA_CHbin+DA_CHbin)/2,stdDA_CHbin,'gx')
% hold on
% [p2 S]= polyfit(EAldbin,(DA_CHbin+DA_CHbin)/2,1)
%  [y delta] = polyval(p2,EAldbin,S);
% hold on
% plot(EAldbin,y,'--g','Linewidth',3)
% box on
% [px2 rmse] = corrcoef(EAldbin,(DA_CHbin+DA_CHbin)/2); 
%  text(.05,.75,[ 'R^2 = ',num2str(px2(1,2).^2,'%.2f') ],'Fontsize',15,'Units','normalized','color','g')
% 
%    text(.2,.75,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','g') 
% 




%%
%%
clear d *bin
d=prctile(EALb(30:210),[0:4:100]);
for i=1:length(d)-1
    id=EALb>d(i) & EALb<d(i+1)
    EALbbin(i)=nanmean(EALb(id));
     DA_TPbin(i)=nanmean(DA_TP(id));
     DA_WHbin(i)=nanmean(DA_WH(id));
     DA_EHbin(i)=nanmean(DA_EH(id));
     DA_HKbin(i)=nanmean(DA_HK(id));
     DA_KKbin(i)=nanmean(DA_KK(id));
     DA_CHbin(i)=nanmean(DA_CH(id));
      stdDA_TPbin(i)=nanstd(DA_TP(id));
     stdDA_WHbin(i)=nanstd(DA_WH(id));
     stdDA_EHbin(i)=nanstd(DA_EH(id));
     stdDA_HKbin(i)=nanstd(DA_HK(id));
     stdDA_KKbin(i)=nanstd(DA_KK(id));
     stdDA_CHbin(i)=nanstd(DA_CH(id));
end


figure(2)
scatter(EALb(30:210),DA_WH(30:210),'MarkerFaceColor','r','MarkerEdgeColor','r',...
    'MarkerFaceAlpha',.3,'MarkerEdgeAlpha',.3)
hold on
errorbar(EALbbin,DA_WHbin,stdDA_WHbin,'rx')
hold on
[p2 S]= polyfit(EALbbin,DA_WHbin,1)
 [y delta] = polyval(p2,EALbbin,S);
hold on
plot(EALbbin,y,'--r','Linewidth',3)
box on
[px2 rmse] = corrcoef(EALbbin,DA_WHbin); 

  text(.6,.85,[ 'r = ',num2str(px2(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','r')

   text(.8,.85,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','r') 




hold on
scatter(EALb(30:210),DA_HK(30:210),'MarkerFaceColor','b','MarkerEdgeColor','b',...
    'MarkerFaceAlpha',.3,'MarkerEdgeAlpha',.3)
hold on
errorbar(EALbbin,DA_HKbin,stdDA_HKbin,'bx')
hold on
[p2 S]= polyfit(EALbbin,DA_HKbin,1)
 [y delta] = polyval(p2,EALbbin,S);
hold on
plot(EALbbin,y,'--b','Linewidth',3)
box on
[px2 rmse] = corrcoef(EALbbin,DA_HKbin); 
 text(.6,.8,[ 'r = ',num2str(px2(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','b')

   text(.8,.8,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','b') 




hold on
scatter(EALb(30:210),DA_KK(30:210),'MarkerFaceColor','m','MarkerEdgeColor','m',...
    'MarkerFaceAlpha',.3,'MarkerEdgeAlpha',.3)
hold on
errorbar(EALbbin,DA_KKbin,stdDA_KKbin,'mx')
hold on
[p2 S]= polyfit(EALbbin,DA_KKbin,1)
 [y delta] = polyval(p2,EALbbin,S);
hold on
plot(EALbbin,y,'--m','Linewidth',3)
box on
[px2 rmse] = corrcoef(EALbbin,DA_KKbin); 
 text(.6,.75,[ 'r = ',num2str(px2(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','m')

   text(.8,.75,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','m') 


%%


clear d *bin
d=prctile(EALb_TP(30:210),[10:4:90]);
for i=1:length(d)-1
    id=EALb_TP>d(i) & EALb_TP<d(i+1)
    EALbbin(i)=nanmean(EALb_TP(id));
     DA_TPbin(i)=nanmean(DA_TP(id));
     DA_WHbin(i)=nanmean(DA_WH(id));
     DA_EHbin(i)=nanmean(DA_EH(id));
     DA_HKbin(i)=nanmean(DA_HK(id));
     DA_KKbin(i)=nanmean(DA_KK(id));
     DA_CHbin(i)=nanmean(DA_CH(id));
      stdDA_TPbin(i)=nanstd(DA_TP(id));
     stdDA_WHbin(i)=nanstd(DA_WH(id));
     stdDA_EHbin(i)=nanstd(DA_EH(id));
     stdDA_HKbin(i)=nanstd(DA_HK(id));
     stdDA_KKbin(i)=nanstd(DA_KK(id));
     stdDA_CHbin(i)=nanstd(DA_CH(id));
end


hold on
scatter(EALb_TP(30:210),DA_TP(30:210),'MarkerFaceColor','k','MarkerEdgeColor','k',...
    'MarkerFaceAlpha',.3,'MarkerEdgeAlpha',.3)
hold on
errorbar(EALbbin,DA_TPbin,stdDA_TPbin,'kx')
hold on
[p2 S]= polyfit(EALbbin,DA_TPbin,1)
 [y delta] = polyval(p2,EALbbin,S);
hold on
plot(EALbbin,y,'--k','Linewidth',3)
box on
[px2 rmse] = corrcoef(EALbbin,DA_TPbin); 
 text(.6,.7,[ 'r = ',num2str(px2(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','k')

   text(.8,.7,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','k') 



hold on
scatter(EALb_TP(30:210),(DA_CH(30:210)+DA_WH(30:210))/2,'MarkerFaceColor','g','MarkerEdgeColor','g',...
    'MarkerFaceAlpha',.3,'MarkerEdgeAlpha',.3)
hold on
errorbar(EALbbin,(DA_CHbin+DA_WHbin)/2,stdDA_EHbin,'gx')
hold on
[p2 S]= polyfit(EALbbin,(DA_CHbin+DA_WHbin)/2,1)
 [y delta] = polyval(p2,EALbbin,S);
hold on
plot(EALbbin,y,'--g','Linewidth',3)
box on
[px2 rmse] = corrcoef(EALbbin,(DA_CHbin+DA_WHbin)/2); 
 text(.6,.65,[ 'r = ',num2str(px2(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','g')

   text(.8,.65,[ 'p = ',num2str(rmse(1,2),'%.2f') ],'Fontsize',15,'Units','normalized','color','g') 


  xlabel('MERRA2 BC mass between 800-500 hPa (\mug/m^3)')
    ylabel('\Delta Albedo due to LAP (%)')
  set(gca,'fontsize',15)