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

%========== Model output ===========%
filename="/data3/yiwen/cesm2/JGR_ocean/Model_output/G_JRA55_test1.pop.h.0031-04.nc"
LON=ncread(filename,'ULONG');LAT=ncread(filename,'ULAT');
depth=ncread(filename,'z_t'); %cm

pathi1='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/';
pathi2='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_linear/'

load([pathi1 'isopycnal_depth_278.mat']);
D_sigma_test1=D_sigma;Uvel_C=U_sigma;Vvel_C=V_sigma;
load([pathi2 'isopycnal_depth_278.mat']);
D_sigma_vdc=D_sigma;Uvel_VDC=U_sigma;Vvel_VDC=V_sigma;

%--------choose pd=27.90(kg/m^3) level --------%
sigma_lev(16)
L=16;  %27.95 

for ii=1:size(Uvel_C,2)
 Uvel{ii}={Uvel_C{ii}(:,:,L) Uvel_VDC{ii}(:,:,L)}
 Vvel{ii}={Vvel_C{ii}(:,:,L) Vvel_VDC{ii}(:,:,L)}
 D_sigma{ii}={D_sigma_test1{ii}(:,:,L) D_sigma_vdc{ii}(:,:,L)}
end

%------- interplation -------%
lat_real=-80:1:90;
lon_real=0:1:360;
[xq,yq] = meshgrid(lon_real,lat_real);
for ii=1:size(Uvel_C,2)
  for Ca=1:2  
Uvel_N{ii}{Ca}=griddata(LON,LAT,Uvel{ii}{Ca},xq,yq);
Vvel_N{ii}{Ca}=griddata(LON,LAT,Vvel{ii}{Ca},xq,yq);
D_sigma_N{ii}{Ca}=griddata(LON,LAT,D_sigma{ii}{Ca},xq,yq);

%----------------------------------------------%
Uvel_N2{ii}{Ca}=circshift(Uvel_N{ii}{Ca},-180,2);
Vvel_N2{ii}{Ca}=circshift(Vvel_N{ii}{Ca},-180,2);
D_sigma_N2{ii}{Ca}=circshift(D_sigma_N{ii}{Ca},-180,2);

%-------1X1-------%
Uvel_N2{ii}{Ca}(:,181:183)=[];
Vvel_N2{ii}{Ca}(:,181:183)=[];
D_sigma_N2{ii}{Ca}(:,181:183)=[];

%-----------------%
Uvel_N2{ii}{Ca}(isnan(Uvel_N2{ii}{Ca}))=0;
Vvel_N2{ii}{Ca}(isnan(Vvel_N2{ii}{Ca}))=0;
 end
end


%================================
% Masking [filter]
%================================
for ii=1:size(Uvel_C,2) %every 10yrs
%--- difference[lat*lon] ----%
D_sigma_diff{ii}=D_sigma_N2{ii}{2}-D_sigma_N2{ii}{1};  
%----set mask value [1,0]----%
mask{ii}=D_sigma_diff{ii};
mask{ii}(~isnan(mask{ii}))=1;mask{ii}(isnan(mask{ii}))=0;
D{ii}=D_sigma_diff{ii};D{ii}(isnan(D{ii}))=0;

%--- difference [VDC-CTRL] ---%
U{ii}=(Uvel_N2{ii}{2}-Uvel_N2{ii}{1}).*mask{ii};
V{ii}=(Vvel_N2{ii}{2}-Vvel_N2{ii}{1}).*mask{ii};

    for smooth=1
        %setup the filter
        filter_size = 3; % 3X3
        filter = ones(filter_size) / (filter_size^2);

        mask{ii}=D_sigma_diff{ii};mask{ii}(~isnan(mask{ii}))=1; %same as original D_diff
        smoothed_D_diff{ii} = imfilter(D{ii}, filter, 'replicate').*mask{ii};
        Uvel_diff{ii} = imfilter(U{ii}, filter, 'replicate');
        Vvel_diff{ii} = imfilter(V{ii}, filter, 'replicate');
    end

end


%========================
% polar vector
%========================
lon2=linspace(-180,180,size(Uvel_diff{1},2));
lat2=linspace(-80,90,size(Uvel_diff{1},1));
%--range_total--%
lon_loc=lon2>=-180 & lon2<=180;
lat_loc=lat2>=-80 & lat2<=90;
[xx,yy]=meshgrid(lon2(lon_loc),lat2(lat_loc));

for ii=1:size(Uvel_C,2)
%---- set varibles ----%
U_value{ii}=Uvel_diff{ii};
V_value{ii}=Vvel_diff{ii};
%--range_polar(1)--%
lat_loc2=lat2 >= 80;
[xx2,yy2]=meshgrid(lon2(lon_loc),lat2(lat_loc2));
UP{ii}=U_value{ii}(lat2>=80,:);VP{ii}=V_value{ii}(lat2>=80,:);

%--range_polar(3)--%
lat_loc3=lat2 >= 86;
[xx3,yy3]=meshgrid(lon2(lon_loc),lat2(lat_loc3));
UP2{ii}=U_value{ii}(lat2>=86,:);VP2{ii}=V_value{ii}(lat2>=86,:);

%--range_polar(4)--%
lat_loc4=lat2 >= 88;
[xx4,yy4]=meshgrid(lon2(lon_loc),lat2(lat_loc4));
UP3{ii}=U_value{ii}(lat2>=88,:);VP3{ii}=V_value{ii}(lat2>=88,:);

%--- repeat value == nan ---%
Uvel_diff{ii}(lat2>=80,:)=nan;
Vvel_diff{ii}(lat2>=80,:)=nan;
UP{ii}(7:11,:)=nan;VP{ii}(7:11,:)=nan;
UP2{ii}(3:5,:)=nan;VP2{ii}(3:5,:)=nan;
UP3{ii}(end,:)=nan;VP3{ii}(end,:)=nan;
end
%====== study region ======%
Xw=[-55 -55 -37.5 -37.5 -55]
Yw=[47.5 57.5 57.5 47.5 47.5]


for vector_factor=1
vecx=[1:2:length(lon2(lon_loc))];vecy=[1:2:length(lat2(lat_loc))]; 
vecx_po=[1:10:length(lon2(lon_loc))];
vecx_po2=[1:20:length(lon2(lon_loc))];
vecx_po3=[1:50:length(lon2(lon_loc))];

vecy_po=[1:2:length(lat2(lat_loc2))];
vecy_po2=[1:2:length(lat2(lat_loc3))];
vecy_po3=[1:2:length(lat2(lat_loc4))];
end

set(0,'defaultfigurecolor','w')

time_labels = {'[1958-1967]','[1968-1977]','[1978-1987]'...
              ,'[1988-1997]','[1998-2007]','[2009-2018]'};

[ha, hfigure]= tight_subplot(2, 3, [.06 .03], [.06 .06], [.03 .1]);

for ii=1:6
  axes(ha(ii));  
m_proj('azimuthal equal-area','latitude',80,'radius',35,'rectbox','on');%rad:45 [27.95]
m_contourf(xx,yy,smoothed_D_diff{ii}(lat_loc,lon_loc),[-700:1:0],'LineColor','none');hold on         
m_coast('patch',[.7 .7 .7],'edgecolor','k');
m_grid('xtick',9,'tickdir','out','ytick',[60 70 80],'Fontsize',12);
m_line(Xw,Yw,'linewi',2,'color','b');hold on
colormap(flipud(brewermap(50, 'YlGn')));
caxis([-350 0]) %27.95
y=text(-.6,-.3,time_labels{ii},'Fontsize',20,'fontweight', 'bold')


set(y,'Rotation',90)
for vec=1
%------------- Vectors -------------%
scale_factor=8.0; %27.95

q1=m_quiver(xx(vecy,vecx),yy(vecy,vecx),Uvel_diff{ii}(vecy,vecx)*scale_factor...
           ,Vvel_diff{ii}(vecy,vecx)*scale_factor,'AutoScale','off','Linewidth',0.5)

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

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

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

%======= 27.95 ========%
q5=m_quiver([90,91],[67.5,67.5],[1.0,1.0]*scale_factor,[0,0],...
            'AutoScale','off','Linewidth',1.5,'ShowArrowHead','on','MaxHeadSize',2)
text(.35,.05,{'1.0'; 'cm/s'},'Fontsize',10,'fontweight', 'bold')  %27.95

q1.Color = 'k';q2.Color = 'k';q3.Color = 'k';q4.Color = 'k';
q5.Color = 'k';
end  

end

%-----gets properties of colorbar-----%
c=colorbar('location','eastoutside','Position', [0.93 0.08 .015 .8],...
          'colormap', colormap(flipud(brewermap(50, 'YlGn'))))
c.Title.String = {"[27.95 kg/m^3]";"m"};
c.FontSize = 14;

