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

UWND=ncread('~/Aerosol_Variaiton/IRUN/read_nc/ncfile/read_nc_UWND.nc','UWND');
VWND=ncread('~/Aerosol_Variaiton/IRUN/read_nc/ncfile/read_nc_VWND.nc','VWND');
year=1950:2000;

V_month=zeros(2,144,96,12,length(year));
dV_month=zeros(2,144,96,12,length(year));
dV_month_fan=zeros(144,96,12,length(year));
V_fan=ones(144,96);
Sita_month=zeros(144,96,12,length(year));
Vw=zeros(2,144,96);
Vs=zeros(2,144,96);


V_month(1,:,:,:,:)=uwnd850;
V_month(2,:,:,:,:)=vwnd850;

Vw(:,:,1:48)=mean(V_month(:,:,1:48,7,:),5);
Vw(:,:,49:96)=mean(V_month(:,:,49:96,1,:),5);
Vs(:,:,1:48)=mean(V_month(:,:,1:48,1,:),5);
Vs(:,:,49:96)=mean(V_month(:,:,49:96,7,:),5);
V=(Vw+Vs)/2;



for yearID=1:length(year)
for monthID=1:12
dV_month(:,:,:,monthID,yearID)=Vw(:,:,:)-V_month(:,:,:,monthID,yearID);
end
end

lat = lat*pi/180;

for i=2:143
for j=2:95
for yearID=1:length(year)
for monthID=1:12
dV_month_fan(i,j,monthID,yearID)= sqrt((dV_month(1,i-1,j,monthID,yearID)^2+dV_month(2,i-1,j,monthID,yearID)^2+4*(dV_month(1,i,j,monthID,yearID)^2+dV_month(2,i,j,monthID,yearID)^2)+dV_month(1,i+1,j,monthID,yearID)^2+dV_month(2,i+1,j,monthID,yearID)^2)*cos(lat(j))+(dV_month(1,i,j-1,monthID,yearID)^2+dV_month(2,i,j-1,monthID,yearID)^2)*cos(lat(j-1))+(dV_month(1,i,j+1,monthID,yearID)^2+dV_month(2,i,j+1,monthID,yearID)^2)*cos(lat(j+1)));
end
end
V_fan(i,j)= sqrt((V(1,i-1,j)^2+V(2,i-1,j)^2+4*(V(1,i,j)^2+V(2,i,j)^2)+V(1,i+1,j)^2+V(2,i+1,j)^2)*cos(lat(j))+(V(1,i,j-1)^2+V(2,i,j-1)^2)*cos(lat(j-1))+(V(1,i,j+1)^2+V(2,i,j+1)^2)*cos(lat(j+1)));
end
end

for yearID=1:length(year)
for monthID=1:12
Sita_month(:,:,monthID,yearID)=dV_month_fan(:,:,monthID,yearID)./V_fan(:,:) -2;
end
end

EASMI_month(1:12,1:length(year)) = mean(mean(Sita_month(45:57,54:70,:,:))); %110E-140E 10N-40N

EASMI_nostandard=zeros(51,1);
EASMI_nostandard(:)=mean(EASMI_month(6:8,:),1);

EASMI=zscore(EASMI_nostandard); %standarded

% save mat file and nc file
save('EASMI.mat','EASMI_All','EASMI_Direct','EASMI_Ctrl')

nccreate('EASMI.nc','Year','Dimensions',{'time',Inf},'Format','classic')
ncwrite('EASMI.nc','Year', year)
ncwriteatt('EASMI.nc','Year','long_name','Year 1950-2000')
nccreate('EASMI.nc','EASMI_All','Dimensions',{'time',Inf})
ncwrite('EASMI.nc','EASMI_All', EASMI)  
ncdisp('EASMI.nc')



