clc;clear all;close all;
addpath('/data3/yiwen/cesm_final_exp/matlab_toolbox/colorbrewer')
addpath('/data3/yiwen/cesm_final_exp/matlab_toolbox/m_map')
addpath('/home/yiwen/matlab_test/toolbox/tight_subplot')

pathi1='/data3/yiwen/matlab_coawst/particle_tracking/0926_tracking/dt_10s/'
pathi2='/data3/yiwen/matlab_coawst/particle_tracking/0926_tracking/'
load([pathi2,'spun_up_10days'])

u_avg1=nanmean(u6(:,:,1152),3);v_avg1=nanmean(v6(:,:,1152),3); 
u_avg2=nanmean(u6(:,:,1153:1212),3);v_avg2=nanmean(v6(:,:,1153:1212),3); 
u_avg3=nanmean(u6(:,:,1213:1272),3);v_avg3=nanmean(v6(:,:,1213:1272),3); 
u_avg4=nanmean(u6(:,:,1273:1332),3);v_avg4=nanmean(v6(:,:,1273:1332),3); 

V1=sqrt(u_avg1(:,2:end).^2+v_avg1(2:end,:).^2);V1(isnan(V1))=999; 
V2=sqrt(u_avg2(:,2:end).^2+v_avg2(2:end,:).^2);V2(isnan(V2))=999;
V3=sqrt(u_avg3(:,2:end).^2+v_avg3(2:end,:).^2);V3(isnan(V3))=999;
V4=sqrt(u_avg4(:,2:end).^2+v_avg4(2:end,:).^2);V4(isnan(V4))=999;


u_avg1(isnan(u_avg1))=999;u_avg2(isnan(u_avg2))=999;
u_avg3(isnan(u_avg3))=999;u_avg4(isnan(u_avg4))=999;
v_avg1(isnan(v_avg1))=999;v_avg2(isnan(v_avg2))=999;
v_avg3(isnan(v_avg3))=999;v_avg4(isnan(v_avg4))=999;

file='Bitou_ocean_his.nc'
%------ check for wind vector ------%
Uwind=ncread(file,'Uwind');Vwind=ncread(file,'Vwind');
Uwind=nanmean(Uwind,3);Vwind=nanmean(Vwind,3);
V=sqrt(Uwind.^2+Vwind.^2);
vecx=Uwind(1:10:end,1:10:end);vecy=Vwind(1:10:end,1:10:end);
vecx(vecx>1.0)=nan;vecy(vecy>1.0)=nan;
scale_factor=0.003;

pcolor(lon,lat,Uwind);shading flat;hold on;colorbar
% q1=quiver(lon(1:10:end,1:10:end),lat(1:10:end,1:10:end),...
%     vecx*scale_factor,vecy*scale_factor,'AutoScale','off');
%%
% Velocity vector (U&V)
vecx1=u_avg1(1:10:end,1:10:end);vecy1=v_avg1(1:10:end,1:10:end);
vecx1(vecx1>1.0)=nan;vecy1(vecy1>1.0)=nan;
vecx2=u_avg2(1:10:end,1:10:end);vecy2=v_avg2(1:10:end,1:10:end);
vecx2(vecx2>1.0)=nan;vecy2(vecy2>1.0)=nan;
vecx3=u_avg3(1:10:end,1:10:end);vecy3=v_avg3(1:10:end,1:10:end);
vecx3(vecx3>1.0)=nan;vecy3(vecy3>1.0)=nan;
vecx4=u_avg4(1:10:end,1:10:end);vecy4=v_avg4(1:10:end,1:10:end);
vecx4(vecx4>1.0)=nan;vecy4(vecy4>1.0)=nan;

V={V1 V2 V3 V4};
vecx={vecx1 vecx2 vecx3 vecx4}
vecy={vecy1 vecy2 vecy3 vecy4}

set(0,'defaultfigurecolor','w')
[ha, hfigure]= tight_subplot(1, 4, [0.005 0.005], [0.13 0.05], [0.05 0.01]);
scale_factor=0.003;
xx=lon(2:end,2:end);yy=lat(2:end,2:end);


time_labels = {'4/9','0hr-6hr',...
               '6hr-12hr','12hr-18hr'};


for ii = 1:4
axes(ha(ii));
pcolor(lon(2:end,2:end),lat(2:end,2:end),V{ii});shading flat;hold on
q1=quiver(lon(1:10:end,1:10:end),lat(1:10:end,1:10:end),vecx{ii}*scale_factor,vecy{ii}*scale_factor,'AutoScale','off');
q2=quiver(121.913,25.123,0.05*scale_factor,0,'AutoScale','off');
text(121.913,25.1228,['0.05 m/s'],'Color',[1 1 1])
q1.Color = 'b';q2.Color = 'b';
title({['U&V vector:' time_labels{ii}]}, 'FontSize', 18);
caxis([0 1.])
axis square;
 if ii>1
    set(gca,'yticklabel',[])
 end
end

c=colorbar('location','southoutside','Position', [0.23 0.1 0.6 0.02],...
 'colormap', colormap(brewermap(20, 'Reds')));c.Title.String = "m/s"; 

%==== save to jpg =====%
set(gcf, 'PaperPositionMode' , 'auto' ) 
print( 'UV' , '-djpeg' , '-r600' )



