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;

%========== EXP & CTRL ===========%
pathi1='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/S_ISOPYCNAL/'
pathi2='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_linear/S_ISOPYCNAL/'

%======== last 10-yrs ========%
load([pathi1,'isopycnal_depth_274.mat']);
S_sigma_test1=S_sigma{6};D_sigma_test1=D_sigma{6};
load([pathi2,'isopycnal_depth_274.mat']);
S_sigma_kpp=S_sigma{6};D_sigma_kpp=D_sigma{6};
%-----%
levIN=21; % !! Choose: 27.6 
Depth_276=D_sigma_test1(:,:,levIN); 
 
%--------- TSR data ---------%
for Model_TSR=1
pathi1='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/';
pathi2='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_linear/'
file1=[pathi1,'TSRyr_000101-006112.nc']
file2=[pathi2,'TSRyr_000101-006112.nc']

salt_test1=ncread(file1,'SALT');
temp_test1=ncread(file1,'TEMP');
pd_test1=sw_dens0(salt_test1,temp_test1);

salt_kpp=ncread(file2,'SALT');
temp_kpp=ncread(file2,'TEMP');
pd_kpp=sw_dens0(salt_kpp,temp_kpp);

n=52:61;
SALT={salt_test1 salt_kpp };
TEMP={temp_test1 temp_kpp};
PD={pd_test1 pd_kpp};

for ii=1:size(SALT,2)
S1{ii}=squeeze(nanmean(SALT{ii}(:,:,:,n),4));
T1{ii}=squeeze(nanmean(TEMP{ii}(:,:,:,n),4));
PD1{ii}=squeeze(nanmean(PD{ii}(:,:,:,n),4));
end
end


S_Var={S1{1} S1{2}}
T_Var={T1{1} T1{2}}
PD_Var={PD1{1} PD1{2}}

temp=ncread(file1,'TEMP');
mask=temp(:,:,:,1);mask(~isnan(mask))=1; 
nx = 320;ny = 384;nz=60;
depth_3d = zeros(nx, ny, nz);
for i = 1:nx
 for j = 1:ny
  depth_3d(i, j, :) = depth_m;
 end
end
depth_3d=depth_3d.*mask;

nx=320;ny=384;
S_sigma=cell(1,size(S_Var,2));
T_sigma=cell(1,size(S_Var,2));
PD_sigma=cell(1,size(S_Var,2));

Depth_IN={Depth_276 Depth_276 }; %CTRL ,KPP 

for ii=1:size(S_Var,2)
 for i=1:nx
    for j=1:ny
  depth0{ii} = squeeze(Depth_IN{ii}(i,j));
  %----------- choose depth of same sigma_lev ------------%
  [~,loc1{ii}]=min(abs(depth_m-depth0{ii}));  %find the nearest value
    
  if numel(loc1{ii})>=1
  S_var{ii} = S_Var{ii}(i,j,loc1{ii});
  T_var{ii} = T_Var{ii}(i,j,loc1{ii});
  R_var{ii} = PD_Var{ii}(i,j,loc1{ii});

  S_sigma{ii}(i,j) =S_var{ii};
  T_sigma{ii}(i,j) =T_var{ii};
  PD_sigma{ii}(i,j)=R_var{ii};
  else
  S_sigma{ii}(i,j) = nan;
  T_sigma{ii}(i,j) = nan;
  PD_sigma{ii}(i,j)= nan;
  end
%----------------------------------%  
    end
 end
end

%======= cal DIFF =======%
S_sigma_diff= S_sigma{2}-S_sigma{1};
T_sigma_diff= T_sigma{2}-T_sigma{1};
PD_sigma_EXP= PD_sigma{2};
%====== 1X1 ======%
lat_real=-80:1:90;
lon_real=0:1:360;
[xq,yq] = meshgrid(lon_real,lat_real);

SALT_N=griddata(LON,LAT,S_sigma_diff,xq,yq); 
TEMP_N=griddata(LON,LAT,T_sigma_diff,xq,yq);
PD_N=griddata(LON,LAT,PD_sigma_EXP,xq,yq);
Depth_N=griddata(LON,LAT,Depth_IN{2},xq,yq);
%--------------------------------------%
SALT_N2=circshift(SALT_N,-180,2);
TEMP_N2=circshift(TEMP_N,-180,2);
PD_N2=circshift(PD_N,-180,2);
Depth_N2=circshift(Depth_N,-180,2);

SALT_N2(:,181:183)=[];
TEMP_N2(:,181:183)=[];
PD_N2(:,181:183)=[];
Depth_N2(:,181:183)=[];

%--------- choose Location ----------%
lon2=linspace(-180,180,size(SALT_N2,2));
lat2=linspace(-80,90,size(SALT_N2,1));
[xx,yy]=meshgrid(lon2,lat2);
%------------------------------------------------%
Xw=[-55 -55 -37.5 -37.5 -55];
Yw=[47.5 57.5 57.5 47.5 47.5];

%======== along Greenland =========%
for line=1
%======== east-line ==========%
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];

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]
end

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

figure;
time_labels = {'KPP_{ICEL} RHO at depth of[PHC2]','KPP_{ICEL} RHO at depth of[CTRL]'};

m_proj('miller','lat',[45,80],'lon',[-100,30]);
m_pcolor(xx,yy,PD_N2-1000);shading interp;hold on;
m_contourf(xx,yy,Depth_N2,[20 40 60 80 100 150],'LineWidth',1.0,'Fill','off','Showtext','on');
% ------------------------------------------------------------------------------%
m_coast('patch','w','edgecolor','k');
m_grid('xtick',6,'ytick',[50 60 70 75],'linestyle','none','tickdir',...
   'out','Linewidth',3,'Fontsize',16);
m_line(Xw,Yw,'linewi',2,'color','k');hold on       %study region
m_line(DS_x1,DS_y1,'linewi',2,'color','r');hold on %DS_East
m_line(DS_x2,DS_y2,'linewi',2,'color','r');hold on %DS_West
m_line(LS_x,LS_y,'linewi',2,'color','b');hold on   %LS
% -------------------------------------------------------------------%
colormap(flipud(brewermap(20, 'RdBu')));
title({[time_labels{2}],['27.6 kg/m^{3}']},'Interpreter','tex');
caxis([27.35 27.85]);
set(gca,'Fontsize',16,'Linewidth',2,'TickDir','out');ytickformat('%.2f')
c=colorbar('location','eastoutside','Position', [0.93 .12 .02 .80],...
       'colormap', colormap(flipud(brewermap(30, 'RdBu'))))
c.Title.String = "kg/m^3"  




