clc;clear all;close all;
addpath('wwseawater')

%========== 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');    %unit=cm
dz=ncread(filename,'dz')         %unit=cm

%========= EXP ==========%
for Model_load=1
pathi1='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/';      %CTRL
pathi3='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_linear/' %KPP_ICEL
pathi4='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_R/'      %KPP_R   

file1=[pathi1,'TSRyr_000101-006112.nc']
file3=[pathi3,'TSRyr_000101-006112.nc']
file4=[pathi4,'TSRyr_000101-006112.nc']

n=48:57; %2005-2014
salt_test1=ncread(file1,'SALT');salt_test1=salt_test1(:,:,:,n);
temp_test1=ncread(file1,'TEMP');temp_test1=temp_test1(:,:,:,n);

salt_ice=ncread(file3,'SALT');salt_ice=salt_ice(:,:,:,n);
temp_ice=ncread(file3,'TEMP');temp_ice=temp_ice(:,:,:,n);
salt_ice01=ncread(file4,'SALT');salt_ice01=salt_ice01(:,:,:,n);
temp_ice01=ncread(file4,'TEMP');temp_ice01=temp_ice01(:,:,:,n);
end

%========= OBS data(1) WOA =========%
WOA_s3='/data3/yiwen/cesm2/JGR_ocean/OBS_data/WOA23/woa23_A5B4_s00_01.nc'; %2005-2014
WOA_t3='/data3/yiwen/cesm2/JGR_ocean/OBS_data/WOA23/woa23_A5B4_t00_01.nc';
lon1=ncread(WOA_s3,'lon');
lat1=ncread(WOA_s3,'lat');
depth1=ncread(WOA_s3,'depth');%meter
depth_loc1=depth1<=650;       %choose depth reange
loc1=find(lat1>=75);
SALT_WOA05=ncread(WOA_s3,'s_an');
TEMP_WOA05=ncread(WOA_t3,'t_an');


for k=1:length(depth1)
TEMP_N(k)=nanmean(nanmean(TEMP_WOA05(:,loc1,k)));
SALT_N(k)=nanmean(nanmean(SALT_WOA05(:,loc1,k)));
end
RHO_N=sw_dens(SALT_N,TEMP_N,0); %kg/m^3

%========= OBS data(2) EN4 =========%
file='/data3/yiwen/cesm2/JGR_ocean/OBS_data/EN4/EN.4.2.2.f.analysis.g10.195801.nc'
pathi3='/data3/yiwen/cesm2/JGR_ocean/OBS_data/EN4/mat_files/'
load([pathi3 'EN4_TS'])
depth_OBS=ncread(file,'depth'); %meter
lat_OBS=ncread(file,'lat')
lon_OBS=ncread(file,'lon')

Var = {SALT_yr TEMP_yr};
periods = [48 57];%2005-2014
for ii = 1:length(Var)
  Var_parts{ii} = squeeze(nanmean(Var{ii}(:,:,:,periods(1):periods(2)), 4));
end

% set_variables
salt_OBS  = Var_parts{1};
temp_OBS  = Var_parts{2};

% 3) choose region
lon_idx = 1:length(lon_OBS);        % all lon
lat_idx = find(lat_OBS >= 75.0);    % >= 75 N
depth_idx = find(depth_OBS <= 750); % 0..750 m

% 4) volume weight（cos(lat)）
lat_sel = lat_OBS(lat_idx);
w_lat = cosd(lat_sel);              %weight ~ cos(latitude)
w_lat_norm = w_lat / nansum(w_lat); %normalize 

% 5) depth profile
temp_sel = temp_OBS(lon_idx, lat_idx, depth_idx);
salt_sel = salt_OBS(lon_idx, lat_idx, depth_idx);

% nanmean over dim 1 (lon) -> result dims 1 x lat x depth -> squeeze -> lat x depth
temp_lonmean = squeeze(nanmean(temp_sel, 1)); % lat x depth
salt_lonmean = squeeze(nanmean(salt_sel, 1)); % lat x depth

% weighted mean over lat
nd = size(temp_lonmean,2); % number of depths selected
TEMP_profile = nan(nd,1);
SALT_profile = nan(nd,1);
for k = 1:nd
    tvec = temp_lonmean(:,k); % lat x 1
    svec = salt_lonmean(:,k);
    valid_t = ~isnan(tvec);
    valid_s = ~isnan(svec);

    if any(valid_t)
        % ignore NaN： sum(w .* t) / sum(w for valid)
        w_used = w_lat .* valid_t;
        TEMP_profile(k) = nansum(tvec .* w_lat) / sum(w_used);
    else
        TEMP_profile(k) = NaN;
    end

    if any(valid_s)
        w_used = w_lat .* valid_s;
        SALT_profile(k) = nansum(svec .* w_lat) / sum(w_used);
    else
        SALT_profile(k) = NaN;
    end
end

% 6) 將 profile 與對應深度配對
depth_sel = depth_OBS(depth_idx); 

% 7) notice unit：sw_dens 要求溫度為 degC
TEMP_OBS = TEMP_profile; % degC
SALT_OBS = SALT_profile; % practical salinity (psu)
RHO_OBS = sw_dens(SALT_OBS, TEMP_OBS, zeros(size(TEMP_OBS))); 


%------------[Model data interp]---------------------%
Var={salt_test1 temp_test1... 
     salt_ice temp_ice...
     salt_ice01 temp_ice01};


for ii=1:size(Var,2)
Var_M{ii}=squeeze(nanmean(Var{ii},4)); %time_mean
end

%==============================%
%range of depth(cm->m)
%==============================%
depth_m=depth./100;
depth_loc=depth_m<=650;          %choose depth reange
loc1=find(LAT>=75);

lat_weights = cosd(LAT);         
area=ncread(filename,'TAREA');   %every layer has diff area[cm^2]
temp1=ncread(filename,'TEMP');
mask=temp1;mask(~isnan(mask))=1;


for k=1:60
areamask{k}=area.*mask(:,:,k).*dz(k);    %2D
areamask_IN{k}=areamask{k}(loc1).*lat_weights(loc1); %1D
sumarea(k)=squeeze(nansum(areamask_IN{k}));

%====== Model data =====%
for ii=1:size(Var_M,2)
Var_k{ii}{k}=Var_M{ii}(:,:,k);
Var_IN{ii}{k}=Var_k{ii}{k}(loc1);
Var_N{ii}(k)=nansum(nansum(Var_IN{ii}{k}.*areamask_IN{k}))./sumarea(k);

end
end

SALT={Var_N{1} Var_N{3} Var_N{5}} 
TEMP={Var_N{2} Var_N{4} Var_N{6}}

for ii=1:size(SALT,2),RHO_N2{ii}=sw_dens(SALT{ii},TEMP{ii},0);end %kg/m^3
%%
set(0, 'defaultfigurecolor', 'w');
C = colormap(lines(6));C=C(2:end,:);
L=linspace(-4,8,2);
L2=linspace(1025,1029,2);
D=[75 75];
D2=[250 250];


t=tiledlayout(1,2)
ax1=axes(t);
for ii=1:size(SALT,2)
line(TEMP{ii}(depth_loc),depth_m(depth_loc),'Color',C(ii,:),...
    'Linewidth',4,'Parent',ax1,'Linewidth',2.5);hold on
end
line(TEMP_N(depth_loc1),depth1(depth_loc1),'Color','g',...
    'Linewidth',4,'Parent',ax1,'Linewidth',3);
line(TEMP_OBS,depth_sel,'Color','c',...
    'Linewidth',4,'Parent',ax1,'Linewidth',3);
line(L,D,'LineStyle','-.','Color',[0 0 0 .45],'Linewidth',1.5,'Parent',ax1)
line(L,D2,'LineStyle','-.','Color',[0 0 0 .45],'Linewidth',1.5,'Parent',ax1)

set(ax1,'XColor','k','YColor','k');
xlabel('deg C','fontsize',16);ylabel('depth(m)');
xticklabel_fs = get(ax1,'XTickLabel');
set(gca, 'YDir', 'reverse', 'FontSize', 16, 'LineWidth', 2, 'FontWeight', 'bold', ...
    'XMinorTick', 'on', 'YMinorTick', 'on', 'TickDir', 'out');
ytickformat('%.2f');
grid on;axis square;hold on
xlim([-1.8 1.8])
ylim([0 600])
h=legend({'CTRL','KPP_{ICEL}','KPP_{R}','WOA','EN4'},'Location','southwest','Interpreter','tex');
legend('boxoff');set(h,'Fontsize',16);
%=============================================%
ax2=axes(t);
ax2.XAxisLocation='top';
ax2.YAxisLocation='right';
ax2.Color='none';
for ii=1:size(SALT,2)
line(SALT{ii}(depth_loc),depth_m(depth_loc),'LineStyle','--',...
     'Color',C(ii,:),'Parent',ax2,'Linewidth',2.5);hold on
end
line(SALT_N(depth_loc1),depth1(depth_loc1),'LineStyle','--',...
     'Color','g','Parent',ax2,'Linewidth',3);hold on
line(SALT_OBS,depth_sel,'LineStyle','--',...
     'Color','c','Parent',ax2,'Linewidth',3);
xlabel('psu','fontsize',16);
xticklabel_fs2 = get(ax2,'XTickLabel');
set(gca, 'YDir', 'reverse', 'FontSize', 16, 'LineWidth', 2, 'FontWeight', 'bold', ...
    'XMinorTick', 'on', 'YMinorTick', 'on', 'TickDir', 'out');
ytickformat('%.2f');
grid on;axis square;hold on
set(ax2,'YTick',[])
ylim([0 600])


nexttile;
for ii=1:size(SALT,2)
P=plot(RHO_N2{ii}(depth_loc),depth_m(depth_loc),'Linewidth',2.5,'color',C(ii,:));hold on
end
P=plot(RHO_N(depth_loc1),depth1(depth_loc1),'Linewidth',3,'color','g');hold on
P=plot(RHO_OBS,depth_sel,'Linewidth',3,'color','c');hold on
line(L2,D,'LineStyle','-.','Color',[0 0 0 .45],'Linewidth',1.5)
line(L2,D2,'LineStyle','-.','Color',[0 0 0 .45],'Linewidth',1.5)

h=legend({'CTRL','KPP_{ICEL}','KPP_{R}','WOA','EN4'},'Location','southwest','Interpreter','tex');
legend('boxoff');set(h,'Fontsize',16);
set(gca,'YDir','reverse')
xlabel('kg/m^3');ylabel('depth(m)');

set(gca,'Fontsize',16,'Linewidth',3,'XMinorTick','on',...
   'YMinorTick','on','TickDir','out','FontWeight','bold');ytickformat('%.2f')
axis square;grid on
xlim([1025 1028.3])
ylim([0 600])



