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');    %unit=cm
depth_m=depth./100;

%==========(1) T&S&PD ==========%
pathi1='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/';      %CTRL
pathi2='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_linear/' %KPP_ICEL
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_vdc=ncread(file2,'SALT');
temp_vdc=ncread(file2,'TEMP');
pd_vdc=sw_dens0(salt_vdc,temp_vdc);

%======== every 10-yrs ========%
n={1:10 11:20 21:30 31:40 41:50 52:61}

%------ year mean -------%
for ii=1:size(n,2)
S1{ii}=squeeze(nanmean(salt_vdc(:,:,:,n{ii}),4));
S2{ii}=squeeze(nanmean(salt_test1(:,:,:,n{ii}),4));
T1{ii}=squeeze(nanmean(temp_vdc(:,:,:,n{ii}),4));
T2{ii}=squeeze(nanmean(temp_test1(:,:,:,n{ii}),4));

PD1{ii}=squeeze(nanmean(pd_vdc(:,:,:,n{ii}),4));
PD2{ii}=squeeze(nanmean(pd_test1(:,:,:,n{ii}),4));

PD_diff{ii}=PD1{ii}(:,:,11)-PD2{ii}(:,:,11);
SALT_diff{ii}=S1{ii}(:,:,11)-S2{ii}(:,:,11);
end

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

for i=1:size(n,2)      
SALT_N{i}=griddata(LON,LAT,SALT_diff{i},xq,yq); 
PD_N{i}=griddata(LON,LAT,PD_diff{i},xq,yq);
%--------------------------------------%
SALT_N2{i}=circshift(SALT_N{i},-180,2);
PD_N2{i}=circshift(PD_N{i},-180,2);

SALT_N2{i}(:,181:183)=[];
PD_N2{i}(:,181:183)=[];
end

%--------- choose Location ----------%
lon2=linspace(-180,180,size(SALT_N2{1},2));
lat2=linspace(-80,90,size(SALT_N2{1},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];

DS_x=[-13 -14 -15 -17 -18 -18 -19 -20 -22 -24 -25 -26 ...
      -27 -29 -31 -32 -35 -37 -39 -39 -40 -41 -41 -42 ...
      -44 -46 -49 -50 -52 -56];

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

%========[DS] & [LS] ========%
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]


set(0,'defaultfigurecolor','w')
[ha, hfigure]= tight_subplot(2, 3, [.06 .05], [.05 .08], [.08 .08]);

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

for ii = 1:6
    axes(ha(ii));
m_proj('miller','lat',[30,80],'lon',[-80,30]);
m_pcolor(xx,yy,SALT_N2{ii});shading interp;hold on;
% ------------------------------------------------------------------------------%
m_coast('patch','w','edgecolor','k');
m_grid('xtick',3,'ytick',[50 65 70 75],'linestyle','none','tickdir',...
   'out','Linewidth',3,'Fontsize',14);
m_line(Xw,Yw,'linewi',2,'color','k');hold on   %study region
m_line(DS_x,DS_y,'linewi',2,'color','r');hold on %DS
m_line(LS_x,LS_y,'linewi',2,'color','b');hold on %LS

% -------------------------------------------------------------------%
colormap(flipud(brewermap(20, 'RdBu')));
title({'SALT Diff.: at depth of',['105m' time_labels{ii}]},'Interpreter','tex');
caxis([-.2 .2]); 

set(gca,'Fontsize',12,'Linewidth',2,'TickDir','out');ytickformat('%.2f')
end

c=colorbar('location','eastoutside','Position', [0.94 0.13 .01 .8],...
          'colormap', colormap(flipud(brewermap(20, 'RdBu'))))
c.FontSize = 14;

