clc;clear all;close all;
addpath('/work/yiwen89419/cesm_POP/matlab_toolbox/colorbrewer')
addpath('/work/yiwen89419/cesm_POP/matlab_toolbox/m_map')

filename='/work/yhtseng00/scratch/cesm2/archive/G_JRA55_test1_vdc013l/ocn/hist/G_JRA55_test1_vdc013l.pop.h.0001-01.nc'
%filename='test.nc'

LON=ncread(filename,'ULONG');LAT=ncread(filename,'ULAT');
lon=LON(:,1);lat=LAT(1,:);

uvel = ncread(filename,'UVEL');  %cm/s
vvel = ncread(filename,'VVEL');  %cm/s 
rot  = ncread(filename,'ANGLE'); %radius
%cal true velocity %
uPop =  uvel.*cos(rot) - vvel.*sin(rot);  
vPop =  uvel.*sin(rot) + vvel.*cos(rot);

uPop=uPop(:,:,1); %surface
vPop=vPop(:,:,1); %surface
%-------- SSH--> U&V --------%
SSH=ncread(filename,'SSH');
%ang=rot.*(180/pi);
%dx=ncread(filename,'DXT');
%dy=ncread(filename,'DYT');
%g = 9.8;
%f = 5e-5;

%fu=zeros(320,384);fv=zeros(320,384);
%for ii=2:320-1
% for jj=2:384-1
% D2(ii,jj) = (SSH(ii,jj+1)-SSH(ii,jj-1))...
%            ./ (dy(ii,jj-1)./2+dy(ii,jj)+ dy(ii,jj+1)./2 ) ; 
% D1(ii,jj) = (SSH(ii+1,jj)-SSH(ii-1,jj))...
%            ./ (dx(ii-1,jj)./2+dx(ii,jj)+ dx(ii+1,jj)./2 ) ;
%--------------------------%
% theta2=ang;
% theta1=90-ang;
% fu(ii,jj) = -g/f.*(sind(theta2(ii,jj)).*D1(ii,jj)+...
%                    sind(theta1(ii,jj)).*D2(ii,jj)); 
% fv(ii,jj) =  g/f.*(cosd(theta2(ii,jj)).*D1(ii,jj)+...
%                    cosd(theta1(ii,jj)).*D2(ii,jj)); 
% end
%end
%-----Interplation------%
lat_real=linspace(-80,90,384);
lon_real=linspace(0,360,320);
[xq,yq] = meshgrid(lon_real,lat_real);

SSH_N=griddata(LON,LAT,SSH,xq,yq,'linear');
Uvel_N=griddata(LON,LAT,uPop,xq,yq,'linear');
Vvel_N=griddata(LON,LAT,vPop,xq,yq,'linear');

SSH_N2=circshift(SSH_N,-180+20,2);
Uvel_N2=circshift(Uvel_N,-180+20,2);
Vvel_N2=circshift(Vvel_N,-180+20,2);

SSH_N2(:,160:162)=[];
Uvel_N2(:,160:162)=[];
Vvel_N2(:,160:162)=[];
%--------------------------------------------
lon2=linspace(-180,180,320-3);
lat2=linspace(-80,90,384);

%--range_total--%
lon_loc=lon2>=-180 & lon2<=180;
lat_loc=lat2>=-80 & lat2<=90;
[xx,yy]=meshgrid(lon2(lon_loc),lat2(lat_loc));

%--range_polar(1)--%
lat_loc2=lat2 >= 80;
[xx2,yy2]=meshgrid(lon2(lon_loc),lat2(lat_loc2));
UP=Uvel_N2(lat2>=80,:);
VP=Vvel_N2(lat2>=80,:);

%--range_polar(3)--%
lat_loc3=lat2 >= 86;
[xx3,yy3]=meshgrid(lon2(lon_loc),lat2(lat_loc3));
UP2=Uvel_N2(lat2>=86,:);
VP2=Vvel_N2(lat2>=86,:);

%--range_polar(4)--%
lat_loc4=lat2 >= 88;
[xx4,yy4]=meshgrid(lon2(lon_loc),lat2(lat_loc4));
UP3=Uvel_N2(lat2>=88,:);
VP3=Vvel_N2(lat2>=88,:);

%-----interp to finer grid ------%
%UP(isnan(UP))=0;VP(isnan(VP))=0;
%for j=1:size(UP,1)
% U_new(j,:) = interp1([1:1:nx],UP(j,:),nt2,'linear');
% V_new(j,:) = interp1([1:1:nx],VP(j,:),nt2,'linear');
%end

%xx_new=linspace(-180,180,3170);
%[X2,Y2]=meshgrid(lon2(lon_loc),lat2(lat_loc2));
%-----------------------%
vecx=[1:5:length(lon2(lon_loc))];
vecy=[1:5:length(lat2(lat_loc))];

vecx_po=[1:15:317];vecy_po=[1:5:23];
vecx_po2=[1:30:317];vecy_po2=[4:8:8];
vecx_po3=[1:180:317];vecy_po3=[1:2:5];


Uvel_N2(lat2>=80,:)=nan;
Vvel_N2(lat2>=80,:)=nan;
UP(14:23,:)=nan;VP(14:23,:)=nan;
UP2(6:10,:)=nan;VP2(6:10,:)=nan;
UP3(end,:)=nan;VP3(end,:)=nan;
%------------------------------------------%
set(0,'defaultfigurecolor','w')  
m_proj('stereographic','lat',90,'long',0,'radius',25);
%m_proj('Azimuthal Equal-area','latitude',90,'radius',45,'rectbox','on');
m_coast('patch',[.7 .7 .7],'edgecolor','k');
m_grid('xtick',12,'tickdir','out','ytick',[70 80 85],'Fontsize',16);hold on
m_contourf(xx,yy,SSH_N2(lat_loc,lon_loc),'LineColor','none');hold on      
%-----gets properties of colorbar-----%
colormap(flipud(brewermap(50, 'RdYlBu')));
c=colorbar('location','southoutside')
caxis([-80 80])
c.Title.String = "cm";
a=get(c);a.Position                          
%-----------------------------------------%
scale_factor=0.3
q1=m_quiver(xx(vecy,vecx),yy(vecy,vecx),Uvel_N2(vecy,vecx)*scale_factor...
           ,Vvel_N2(vecy,vecx)*scale_factor,'AutoScale','off','Linewidth',1.5);hold on

q2=m_quiver(xx2(vecy_po,vecx_po),yy2(vecy_po,vecx_po),UP(vecy_po,vecx_po)*scale_factor...
           ,VP(vecy_po,vecx_po)*scale_factor,'AutoScale','off','Linewidth',1.5);hold on

q3=m_quiver(xx3(vecy_po2,vecx_po2),yy3(vecy_po2,vecx_po2),UP2(vecy_po2,vecx_po2)*scale_factor...
           ,VP2(vecy_po2,vecx_po2)*scale_factor,'AutoScale','off','Linewidth',1.5);hold on

q4=m_quiver(xx4(vecy_po3,vecx_po3),yy4(vecy_po3,vecx_po3),UP3(vecy_po3,vecx_po3)*scale_factor...
           ,VP3(vecy_po3,vecx_po3)*scale_factor,'AutoScale','off','Linewidth',1.5);hold on

%-----------------------------------------%
q5=m_quiver(90,70,20*scale_factor,0,'AutoScale','off','Linewidth',3,'ShowArrowHead','on')
q1.Color = 'k';q2.Color = 'k';q3.Color = 'k';q4.Color = 'k';
q5.Color = 'w';
text(.25,.03,'20 cm/s','Fontsize',12,'fontweight', 'bold') 




