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

%========== Model output ===========%
filename="/data3/yiwen/cesm2/JGR_ocean/Model_output/G_JRA55_test1.pop.h.0031-04.nc"
LAT=ncread(filename,'TLAT');LON=ncread(filename,'TLONG');
depth=ncread(filename,'z_t');depth=double(depth);  %unit=cm
depth_m=depth./100;

%==========(1) T&S&PD ==========%
pathi1='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_linear/'
pathi2='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/';
file1=[pathi1,'TSRyr_000101-006112.nc']
file2=[pathi2,'TSRyr_000101-006112.nc']

for loading_TSR=1
%====== EXP ======%
salt_vdc=ncread(file1,'SALT');
temp_vdc=ncread(file1,'TEMP');
pd_vdc=sw_dens0(salt_vdc,temp_vdc);

%====== CTRL =======%
salt_test1=ncread(file2,'SALT');
temp_test1=ncread(file2,'TEMP');
pd_test1=sw_dens0(salt_test1,temp_test1);

%------ last 10-yrs -------%
n=52:61;
S1=squeeze(nanmean(salt_vdc(:,:,:,n),4));
S2=squeeze(nanmean(salt_test1(:,:,:,n),4));
T1=squeeze(nanmean(temp_vdc(:,:,:,n),4));
T2=squeeze(nanmean(temp_test1(:,:,:,n),4));
PD1=squeeze(nanmean(pd_vdc(:,:,:,n),4));
PD2=squeeze(nanmean(pd_test1(:,:,:,n),4));
%===== diff[EXP- CTRL] ======%
PD_diff=PD1-PD2; 
SALT_diff=S1-S2;
end

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

 for k=1:60    
    SALT1_N(:,:,k)=griddata(LON,LAT,SALT_diff(:,:,k),xq,yq); 
    PD_L(:,:,k)=griddata(LON,LAT,PD2(:,:,k),xq,yq);   %CTRL
    PD_VDC(:,:,k)=griddata(LON,LAT,PD1(:,:,k),xq,yq); %VDC_R / KPP_ICE0
 end

SALT1_N2=circshift(SALT1_N,-180,2);
PD_L2=circshift(PD_L,-180,2);
PD_VDC2=circshift(PD_VDC,-180,2);

SALT1_N2(:,181:183,:)=[];
PD_L2(:,181:183,:)=[];
PD_VDC2(:,181:183,:)=[];
end


lat2=linspace(-80,90,size(SALT1_N2,1));
lon2=linspace(-180,180,size(SALT1_N2,2));
%===== point of incline line =====%
for Incline_Line=1
%======== east current ========%
DS_x1=[-13 -14 -15 -17 -18 -18 -19 -20 -22 -24 -25 -26 ...
       -27 -29 -31 -32 -35 -37 -39 -39 -40 -41 -41 -42];

DS_y1=[75 74 73 72 71 70 69 68 67 67 67 67 ...
       66 66 66 65 64 63 62 61 61 60 59 58];

%======== west current ========%
DS_x2=[-42 -44 -46 -49 -50 -52 -56];
DS_y2=[58 59 59 60 61 63 65];

for ii=1:length(DS_y1)
[~,DS_latIN1(ii)]=min(abs(lat2-DS_y1(ii)));
[~,DS_lonIN1(ii)]=min(abs(lon2-DS_x1(ii)));
end
for ii=1:length(DS_y2)
[~,DS_latIN2(ii)]=min(abs(lat2-DS_y2(ii)));
[~,DS_lonIN2(ii)]=min(abs(lon2-DS_x2(ii)));
end

%------(2)------%
LS_x=[-56 -56 -57 -58 -58 -57 -57 -56 -55 -53 -52 -50 ...
      -50 -50 -51 -50 -49]
LS_y=[65 64 63 62 61 60 59 58 57 56 56 56 55 54 52 51 50]

for ii=1:length(LS_y)
[~,LS_latIN(ii)]=min(abs(lat2-LS_y(ii)));
[~,LS_lonIN(ii)]=min(abs(lon2-LS_x(ii)));
end

%----------------------------------------------------------%

 for x=1:length(DS_lonIN1)
 SALT_depth_M(x,:)=squeeze(SALT1_N2(DS_latIN1(x),DS_lonIN1(x),:)); 
 PD_Line(x,:)     =squeeze(PD_L2(DS_latIN1(x),DS_lonIN1(x),:));   %east
 end

 for x=1:length(DS_lonIN2)
 SALT_depth_MW(x,:)=squeeze(SALT1_N2(DS_latIN2(x),DS_lonIN2(x),:));     
 PD_Line_VDC(x,:) =squeeze(PD_VDC2(DS_latIN2(x),DS_lonIN2(x),:)); %west
 end

 for x=1:length(LS_lonIN)
 SALT_depth_M2(x,:)=squeeze(SALT1_N2(LS_latIN(x),LS_lonIN(x),:)); 
 PD_Line2(x,:)=squeeze(PD_L2(LS_latIN(x),LS_lonIN(x),:)); 
 PD_Line_VDC2(x,:)=squeeze(PD_VDC2(LS_latIN(x),LS_lonIN(x),:));
 end
end

%%
%========================= [DS] =============================%
set(0,'defaultfigurecolor','w')
[xx1,yy1]=meshgrid(lat2(DS_latIN1),depth_m); %east
[xx2,yy2]=meshgrid(lat2(DS_latIN2),depth_m); %west
D1=300; D2=3000;
depth_loc1=depth_m<=D1;
depth_loc2=depth_m>=0 & depth_m<=D2;
time_labels = {'Difference [2009-2018]'};

function set_fig_common()
    set(gca,'Fontsize',14,'Linewidth',2,'TickDir','out',...
        'XMinorTick','on','YMinorTick','on','YDir','reverse');
    ytickformat('%.2f');
    colormap(flipud(brewermap(70, 'RdBu')));
    caxis([-.2 .2]);
end

figure(1); clf;
set(gcf, 'Position', [100, 100, 850, 650]);
t = tiledlayout(2, 2, 'TileSpacing','none', 'Padding','loose');

% ===== 左上 (WEST, shallow) =====
nexttile(1)
pcolor(xx2(depth_loc1,:), yy2(depth_loc1,:), SALT_depth_MW(:,depth_loc1)'); shading flat; hold on
[C2,H2]=contourf(xx2(depth_loc1,:), yy2(depth_loc1,:), PD_Line_VDC(:,depth_loc1)'-1000, [27.6 27.95],...
    'LineColor','g','LineWidth',2.0,'Fill','off','Showtext','on','LabelSpacing',400);
clabel(C2,H2,'Fontsize',14,'fontweight','bold','Color',[0.4660 0.6740 0.1880])
xlim([60 65]); ylim([0, D1-10]);
set(gca,'XDir','reverse','xticklabel',{[]})
set_fig_common()
title('Difference [2009-2018]')
% ylabel('Depth(m)')

% ===== 右上 (EAST, shallow) =====
nexttile(2)
pcolor(xx1(depth_loc1,:), yy1(depth_loc1,:), SALT_depth_M(:,depth_loc1)'); shading flat; hold on
[C2,H2]=contourf(xx1(depth_loc1,:), yy1(depth_loc1,:), PD_Line(:,depth_loc1)'-1000, [27.6 27.95],...
    'LineColor','g','LineWidth',2.0,'Fill','off','Showtext','on','LabelSpacing',400);
clabel(C2,H2,'Fontsize',14,'fontweight','bold','Color',[0.4660 0.6740 0.1880])
xlim([60 75]); ylim([0, D1-10]);
set(gca,'xticklabel',{[]})
set_fig_common()
set(gca,'yticklabel',{[]})
title('Difference [2009-2018]')

% ===== 左下 (WEST, deep) =====
nexttile(3)
pcolor(xx2(depth_loc2,:), yy2(depth_loc2,:), SALT_depth_MW(:,depth_loc2)'); shading flat; hold on
[C2,H2]=contourf(xx2(depth_loc2,:), yy2(depth_loc2,:), PD_Line_VDC(:,depth_loc2)'-1000, [27.6 27.95],...
    'LineColor','g','LineWidth',2.0,'Fill','off','Showtext','on');
clabel(C2,H2,'Fontsize',14,'fontweight','bold','Color',[0.4660 0.6740 0.1880])
xlim([60 65]); ylim([D1, D2]);
set(gca,'XDir','reverse')
xlabel('lat')
set(gca,'xticklabel',{[]}, 'yticklabel',{[]})
set_fig_common()
text(65, 2500, {'(WEST)';'Denmark';'Strait'}, 'FontSize',16,...
    'FontWeight','bold','Color','r')

% ===== 右下 (EAST, deep) =====
nexttile(4)
pcolor(xx1(depth_loc2,:), yy1(depth_loc2,:), SALT_depth_M(:,depth_loc2)'); shading flat; hold on
[C2,H2]=contourf(xx1(depth_loc2,:), yy1(depth_loc2,:), PD_Line(:,depth_loc2)'-1000, [27.6 27.95],...
    'LineColor','g','LineWidth',2.0,'Fill','off','Showtext','on');
clabel(C2,H2,'Fontsize',14,'fontweight','bold','Color',[0.4660 0.6740 0.1880])
xlim([60 75]); ylim([D1, D2]);
xlabel('lat')
set(gca,'xticklabel',{[]})
set_fig_common()
set(gca,'yticklabel',{[]})
text(61, 2500, {'(EAST)';'Denmark';'Strait'}, 'FontSize',16,...
    'FontWeight','bold','Color','r')

% ===== 手動調整 axes 位置 =====
left_margin  = 0.20;
right_margin = 0.20;
mid_gap      = 0;
total_width  = 1 - left_margin - right_margin;
W_width = total_width * (5/20);   % West 5度
E_width = total_width * (15/20);  % East 15度

top_margin    = 0.08;
bottom_margin = 0.08;
row_gap       = 0;
total_height  = 1 - top_margin - bottom_margin;
upper_height  = total_height * 0.40;  % 上排 (0~300m)
lower_height  = total_height * 0.60;  % 下排 (300~3000m)

left_W    = left_margin;
left_E    = left_margin + W_width + mid_gap;
bot_upper = bottom_margin + lower_height + row_gap;
bot_lower = bottom_margin;

axs = findall(gcf, 'type', 'axes');
% findall 順序是反向的: axs(1)=右下, axs(2)=左下, axs(3)=右上, axs(4)=左上
set(axs(4), 'Position', [left_W, bot_upper, W_width, upper_height]);  % 左上
set(axs(3), 'Position', [left_E, bot_upper, E_width, upper_height]);  % 右上
set(axs(2), 'Position', [left_W, bot_lower, W_width, lower_height]);  % 左下
set(axs(1), 'Position', [left_E, bot_lower, E_width, lower_height]);  % 右下


% ylabel
annotation('textbox', [0.045, 0.46, 0.06, 0.08], ...
    'String', 'Depth(m)', ...
    'FontSize', 14, ...
    'FontWeight', 'bold', ...
    'EdgeColor', 'none', ...   
    'HorizontalAlignment', 'center', ...
    'VerticalAlignment', 'middle', ...
    'Rotation', 90);
% ===== Colorbar =====
c = colorbar('location','eastoutside','TicksMode','auto',...
    'Position',[0.92 0.08 0.018 0.84],...
    'colormap', colormap(flipud(brewermap(50,'RdBu'))));
c.Title.String = "psu";
c.Title.FontSize = 14;
c.LineWidth = 1.5;


%%
%========================= [LS] =============================%

[xx1,yy1]=meshgrid(lat2(LS_latIN),depth_m);
time_labels = {'Difference [2009-2018]'};
%---------------------------------------------%

figure(2)
for Labrador_Sea=1
t = tiledlayout(2,1,'TileSpacing','none');
nexttile

pcolor(xx1(depth_loc1,:),yy1(depth_loc1,:),SALT_depth_M2(:,depth_loc1)');shading flat;hold on
[C2,H2]=contourf(xx1(depth_loc1,:),yy1(depth_loc1,:),PD_Line_VDC2(:,depth_loc1)'-1000,[27.6:.1:27.7],...
            'LineColor','g','LineWidth',2.0,'Fill','off','Showtext','off');
clabel(C2,H2,'Fontsize',12,'fontweight', 'bold','Color',[0.4660 0.6740 0.1880])
ylim([0, D1-10]);set(gca,'xticklabel',{[]})
set(gca,'Fontsize',14,'Linewidth',2,'TickDir','out',...
    'XMinorTick','on','YMinorTick','on','YDir','reverse');
ytickformat('%.2f');
colormap(flipud(brewermap(70, 'RdBu')));
axis square;
xlim([50 65]);
caxis([-.2 .2])
title(time_labels);
%-------------------------------------------------%
nexttile
pcolor(xx1(depth_loc2,:),yy1(depth_loc2,:),SALT_depth_M2(:,depth_loc2)');shading flat;hold on
[C2,H2]=contourf(xx1(depth_loc2,:),yy1(depth_loc2,:),PD_Line2(:,depth_loc2)'-1000,[27.6:.1:27.7],...
            'LineColor','g','LineWidth',2.0,'Fill','off','Showtext','off');
xlabel('lat');ylim([D1, D2]);

label_h = ylabel('Depth(m)');
label_h.Position(1) = 45;   
label_h.Position(2) = 250;  
set(gca,'Fontsize',14,'Linewidth',2,'TickDir','out',...
    'XMinorTick','on','YMinorTick','on','YDir','reverse');
ytickformat('%.2f');
colormap(flipud(brewermap(70, 'RdBu')));
axis square;
xlim([50 65]);
caxis([-.2 .2])
text(51, 2700,{'Labrador' ; 'Sea'},'FontSize',20,...
     'FontWeight','bold','Color','b','Interpreter','none');hold on

c=colorbar('location','eastoutside','TicksMode','auto','Position', [.65 .1 .015 .8],...
           'colormap',colormap(flipud(brewermap(50, 'RdBu'))))
c.Title.String = "psu";
c.Title.FontSize=14;
c.LineWidth = 1.5;
end

