clc;close all;clear all;
addpath('/home/yiwen/matlab_test/tight_subplot')
addpath('/data2/summer2021/cbrewer')

cd '/data3/yiwen/cesm_final_exp/plot_test/MLD/'

ncTOPO ='/data3/yiwen/cesm_final_exp/GTJRAv2.10c_v2/control_run/hist/TOPO.nc';
ncid=netcdf.open(ncTOPO,'NC_NOWRITE');
varid = netcdf.inqVarID(ncid,'mask'); mask=netcdf.getVar(ncid,varid,'double'); % size=320x288x45
varid = netcdf.inqVarID(ncid,'lon'); lon=netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'lat'); lat=netcdf.getVar(ncid,varid,'double');
zlevel=ncread(ncTOPO,'lev_c');
zlevel=zlevel/100;  % turn cm tSo m
netcdf.close(ncid)

load ('GSHHS_COAST220HL1.mat'); %grid::

lon=ncread('mld_DR003_c1m_reg2.0.nc','lon');
lat=ncread('mld_DR003_c1m_reg2.0.nc','lat');
mld=ncread('mld_DR003_c1m_reg2.0.nc','mld');
MLDobs=nanmean(mld,3);
MLDobs(MLDobs>2000)=nan;
MLDobs(MLDobs<0)=nan;
[XX,YY]=meshgrid(lon,lat);

contour=[0 20 60 100 200 500 1500 2200]; %this is desired colorbar
pseudo_contour=[1:length(contour)];
% meanMLDa_new=interp1(contour,pseudo_contour,meanMLDa,'pchip'); %scale data to pseudo colorbar
% meanMLDb_new=interp1(contour,pseudo_contour,meanMLDb,'pchip'); %scale data to pseudo colorbar
meanMLDobs_new_TOPO=interp1(contour,pseudo_contour,MLDobs,'pchip');

%=====================================================
%loading EXP variables
%=====================================================
mask(mask==0)=nan;
pathi='/data3/yiwen/cesm_final_exp/plot_test/check_diff_cycle/WIND_EXP/MLD/';
load([pathi,'MLD_MAR_control.mat'])
A=squeeze(nanmean(MLDb,3));
A=mask(:,:,1).*A;
meanMLDb_control=interp1(contour,pseudo_contour,A,'pchip'); %scale data to pseudo colorbar

load([pathi,'MLD_MAR_v3.mat'])
B=squeeze(nanmean(MLDb,3));
B=mask(:,:,1).*B;
meanMLDb_v3=interp1(contour,pseudo_contour,B,'pchip'); %scale data to pseudo colorbar

load([pathi,'MLD_MAR_v4.mat'])
C=squeeze(nanmean(MLDb,3));
C=mask(:,:,1).*C;
meanMLDb_v4=interp1(contour,pseudo_contour,C,'pchip'); %scale data to pseudo colorbar


% diff3=C-A;%C=MAR
% diff9=B-A;
diffv3= B-A;
diffv4= C-A;
%% <plot control run> 
mask(mask==0)=nan;
lat_idx=lat>=0;
lon_idx=lon>270&lon<360;
[xx ,yy]=meshgrid(lon(lon_idx),lat(lat_idx));

%======================================
% Y=log2(x) --> 2^Y=X
%======================================
control=log2(A);
control(control == -inf) = 0;

m_proj('miller','lat',[0,85],'lon',[270,360]);
m_pcolor(xx,yy,control(lon_idx,lat_idx)');shading interp
m_coast('patch','w','edgecolor','k');
m_grid('box','fancy','linestyle','none');
title({'control MLD(MAR) [m]'},'Fontsize',12)
colormap(brewermap(30, 'OrRd'));
c = colorbar('location','eastoutside');
caxis([0 10])
c.Label.String = 'log2(m)';
c.Label.VerticalAlignment = 'top';
%% plot wind0.1 & wind1.3 original data
%======================================
% Y=log2(x) --> 2^Y=X
%======================================
V4=log2(C);V3=log2(B);
V4(V4 == -inf) = 0;
V3(V3 == -inf) = 0;


[ha, hfigure]= tight_subplot(2, 1, [0.13 0.15], [0.11 0.11], [0.1 0.1]);%, 20, 30);

for ii = 1:2
    axes(ha(ii));
    switch ii
      case 1
      m_proj('miller','lat',[0,85],'lon',[270,360]);
      m_pcolor(xx,yy,V4(lon_idx,lat_idx)');shading interp
      m_coast('patch','w','edgecolor','k');
      m_grid('box','fancy','linestyle','none');
      title({'wind0.1 MLD(MAR) [m]'},'Fontsize',12)
      colormap(brewermap(30, 'OrRd'));
      caxis([0 12])
      
      case 2
      m_proj('miller','lat',[0,85],'lon',[270,360]);
      m_pcolor(xx,yy,V3(lon_idx,lat_idx)');shading interp
      m_coast('patch','w','edgecolor','k');
      m_grid('box','fancy','linestyle','none');
      title({'wind1.3 MLD(MAR) [m]'},'Fontsize',12)
      colormap(brewermap(30, 'OrRd'));
      caxis([0 10])
             
    end
end

c=colorbar('location','eastoutside','Position', [0.62 0.207 0.01 0.615],...
         'colormap', colormap(brewermap(30, 'OrRd')))   
caxis([0 10])
c.Label.String = 'log2(m)';
c.Label.VerticalAlignment = 'top';
%% plot EXP-control MLD value
[ha, hfigure]= tight_subplot(2, 1, [0.13 0.15], [0.11 0.11], [0.1 0.1]);%, 20, 30);

for ii = 1:2
    axes(ha(ii));
    switch ii    
      case 1
      m_proj('miller','lat',[0,85],'lon',[270,360]);
      m_pcolor(xx,yy,diffv4(lon_idx,lat_idx)');shading interp
      m_coast('patch','w','edgecolor','k');
      m_grid('box','fancy','linestyle','none');
      title({'\Delta wind0.1 MLD(MAR) [m]'},'Fontsize',12)
      colormap(flipud(brewermap(30, 'RdBu')));
      caxis([-150 150])
      
      case 2
      m_proj('miller','lat',[0,85],'lon',[270,360]);
      m_pcolor(xx,yy,diffv3(lon_idx,lat_idx)');shading interp
      m_coast('patch','w','edgecolor','k');
      m_grid('box','fancy','linestyle','none');
      title({'\Delta wind1.3 MLD(MAR) [m]'},'Fontsize',12)
      colormap(flipud(brewermap(30, 'RdBu')));
      caxis([-150 150])
      
    end
end

c=colorbar('location','eastoutside','Position', [0.62 0.207 0.01 0.615],...
         'colormap', colormap(flipud(brewermap(30, 'RdBu'))))   
caxis([-150 150])
c.Label.String = '(m)';
c.Label.VerticalAlignment = 'top';
     