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');

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,:,:,:,:)=UWND;
V_month(2,:,:,:,:)=VWND;

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

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

EAWMI_nostandard=zeros(length(year),1);

EAWMI_nostandard(1)=mean(EAMI_month([1,2],1),1);
for yearID=2:150
    EAWMI_nostandard(yearID)=(EAMI_month(12,yearID-1)+EAMI_month(1,yearID)+EAMI_month(2,yearID))/3;
end

EAWMI=zscore(EAWMI_nostandard); %standarded

% save mat file and nc file
save('EAWMI_IRUN.mat','EAWMI')

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



