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{6};Uvel_C=U_sigma{6};Vvel_C=V_sigma{6};
load([pathi2 'isopycnal_depth_278.mat']);
D_sigma_kpp=D_sigma{6};Uvel_kpp=U_sigma{6};Vvel_kpp=V_sigma{6};

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

D_sigma={D_sigma_test1(:,:,L) D_sigma_kpp(:,:,L)}
U_sigma={Uvel_C(:,:,L) Uvel_kpp(:,:,L)}
V_sigma={Vvel_C(:,:,L) Vvel_kpp(:,:,L)}

%------- interplation -------%
lat_real=-80:1:90;
lon_real=0:1:360;
[xq,yq] = meshgrid(lon_real,lat_real);

for ii=1:size(D_sigma,2) 
Uvel_N{ii}=griddata(LON,LAT,U_sigma{ii},xq,yq);
Vvel_N{ii}=griddata(LON,LAT,V_sigma{ii},xq,yq);
D_sigma_N{ii}=griddata(LON,LAT,D_sigma{ii},xq,yq);

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

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

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


%================================
% Masking [filter]
%================================
%--- difference[lat*lon] ----%
D_sigma_diff=D_sigma_N2{2}-D_sigma_N2{1};  
U_diff=Uvel_N2{2}-Uvel_N2{1};
V_diff=Vvel_N2{2}-Vvel_N2{1};


%----set mask value [1,0]----%
mask=D_sigma_diff;
mask(~isnan(mask))=1;mask(isnan(mask))=0;
D=D_sigma_diff;D(isnan(D))=0;

%--- difference [EXP-CTRL] ---%
U=U_diff.*mask;
V=V_diff.*mask;

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

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

%========================
% polar vector
%========================
lon2=linspace(-180,180,size(Uvel_diff,2));
lat2=linspace(-80,90,size(Uvel_diff,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 vector=1
%---- set varibles ----%
U_value=Uvel_diff;
V_value=Vvel_diff;
%--range_polar(1)--%
lat_loc2=lat2 >= 80;
[xx2,yy2]=meshgrid(lon2(lon_loc),lat2(lat_loc2));
UP=U_value(lat2>=80,:);VP=V_value(lat2>=80,:);

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

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

%--- repeat value == nan ---%
Uvel_diff(lat2>=80,:)=nan;
Vvel_diff(lat2>=80,:)=nan;
UP(7:11,:)=nan;VP(7:11,:)=nan;
UP2(3:5,:)=nan;VP2(3:5,:)=nan;
UP3(end,:)=nan;VP3(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')
 
m_proj('azimuthal equal-area','latitude',80,'radius',45,'rectbox','on');% 27.95
m_contourf(xx,yy,smoothed_D_diff(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',[45 60 70 80],'Fontsize',16);
m_line(Xw,Yw,'linewi',2,'color','b');hold on
colormap(flipud(brewermap(50, 'YlGn')));

for vec=1
scale_factor=12.0; %27.95
%------------- Vectors -------------%
q1=m_quiver(xx(vecy,vecx),yy(vecy,vecx),Uvel_diff(vecy,vecx)*scale_factor...
           ,Vvel_diff(vecy,vecx)*scale_factor,'AutoScale','off','Linewidth',0.5)

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',0.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',0.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',0.5);hold on
       
%======= 27.95 ========%
q5=m_quiver([90,90],[67.5,67.5],[.5,.5]*scale_factor,[0,0],...
            'AutoScale','off','Linewidth',1.5,'ShowArrowHead','on','MaxHeadSize',2)
text(.38,.1,{'0.5'; 'cm/s'},'Fontsize',14,'fontweight', 'bold') 

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

%==== KPP_ICE ====%
caxis([-350 0]) %27.95
%-----gets properties of colorbar-----%
c=colorbar('location','eastoutside','Position', [0.78 0.12 .015 .75],...
          'colormap', colormap(flipud(brewermap(50, 'YlGn'))))
c.Title.String = {"[27.95 kg/m^3]";"m"};
c.FontSize = 14;






