clc
clear all
close all

exp_name = 'IDS_fcst_CTRL_M09_2007_ens1';

gl_path = ['/discover/nobackup/elee15/GEOSldas_4_0/sims/IDS/fcst_ctrl/' exp_name '/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/'];
fz_path = '/discover/nobackup/projects/geoscm/fzeng/NOBACKUP/Catchment/SMAP_EASEv2_M09/e0004s_transientCO2_05/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/';

start_year = 2007;
stop_year = 2007;
start_month = 1;
total_months = (stop_year-start_year+1)*12;

start_year_str = num2str(start_year);
if start_month >= 10
   start_month_str = num2str(start_month);
else
   start_month_str = ['0' num2str(start_month)];
end

load('index_cell.mat')

gl_rst_file = [gl_path 'Y' start_year_str '/M' start_month_str '/' exp_name '.catchcn_internal_rst.' start_year_str start_month_str '01_0000'];
fz_rst_file = [fz_path 'Y' start_year_str '/M' start_month_str '/e0004s_transientCO2.ens0000.catchcn_ldas_rst.' start_year_str start_month_str '01_0000z'];

gl_poros_data = ncread(gl_rst_file,'POROS');
gl_poros_data = double(gl_poros_data);
fz_poros_data = ncread(fz_rst_file,'POROS');
fz_poros_data = double(fz_poros_data);

gl_tile_data = ncread(gl_rst_file,'TILE_ID');
gl_tile_data = double(gl_tile_data);
fz_tile_data = ncread(fz_rst_file,'TILE_ID');
fz_tile_data = double(fz_tile_data);

gl_ity_data = ncread(gl_rst_file,'ITY');
gl_ity_data = double(gl_ity_data);
fz_ity_data = ncread(fz_rst_file,'ITY');
fz_ity_data = double(fz_ity_data);

gl_tc_data = ncread(gl_rst_file,'TC');
gl_tc_data = double(gl_tc_data);
fz_tc_data = ncread(fz_rst_file,'TC');
fz_tc_data = double(fz_tc_data);

gl_rzexc_data = ncread(gl_rst_file,'RZEXC');
gl_rzexc_data = double(gl_rzexc_data);
fz_rzexc_data = ncread(fz_rst_file,'RZEXC');
fz_rzexc_data = double(fz_rzexc_data);

gl_catdef_data = ncread(gl_rst_file,'CATDEF');                                             
gl_catdef_data = double(gl_catdef_data);                                                   
fz_catdef_data = ncread(fz_rst_file,'CATDEF');                                             
fz_catdef_data = double(fz_catdef_data); 

gl_ghtcnt1_data = ncread(gl_rst_file,'GHTCNT1');                                             
gl_ghtcnt1_data = double(gl_ghtcnt1_data);                                                   
fz_ghtcnt1_data = ncread(fz_rst_file,'GHTCNT1');                                             
fz_ghtcnt1_data = double(fz_ghtcnt1_data); 

gl_htsnnn1_data = ncread(gl_rst_file,'HTSNNN1');
gl_htsnnn1_data = double(gl_htsnnn1_data);                                                              
fz_htsnnn1_data = ncread(fz_rst_file,'HTSNNN1');                                                        
fz_htsnnn1_data = double(fz_htsnnn1_data);  

gl_cnpft1_data = ncread(gl_rst_file,'CNPFT');                                         
gl_cnpft1_data = double(gl_cnpft1_data(:,864));                                                
fz_cnpft1_data = ncread(fz_rst_file,'CNPFT');                                                        
fz_cnpft1_data = double(fz_cnpft1_data(:,864));

gl_cnpft2_data = ncread(gl_rst_file,'CNPFT');   
gl_cnpft2_data = double(gl_cnpft2_data(:,865));                                            
fz_cnpft2_data = ncread(fz_rst_file,'CNPFT');                                            
fz_cnpft2_data = double(fz_cnpft2_data(:,865));

gl_cnpft3_data = ncread(gl_rst_file,'CNPFT');   
gl_cnpft3_data = double(gl_cnpft3_data(:,866));                                            
fz_cnpft3_data = ncread(fz_rst_file,'CNPFT');                                            
fz_cnpft3_data = double(fz_cnpft3_data(:,866));

gl_vgmax_data = ncread(gl_rst_file,'VGWMAX');
gl_vgmax_data = double(gl_vgmax_data);                                                   
fz_vgmax_data = ncread(fz_rst_file,'VGWMAX');
fz_vgmax_data = double(fz_vgmax_data);   

% read LDASsa coordinate data

coord_path = '/gpfsm/dnb05/projects/p144/fzeng/NOBACKUP/Catchment/SMAP_EASEv2_M09/e0004s_transientCO2_05/output/SMAP_EASEv2_M09_GLOBAL/rc_out/';

tileCoordFile = [coord_path 'e0004s_transientCO2.ldas_tilecoord.bin'];
[tile_coord ] = read_tilecoord( tileCoordFile, 0 );

ld_points = length(tile_coord.com_lat);
ld_lats = tile_coord.com_lat(1:ld_points);
ld_lons = tile_coord.com_lon(1:ld_points);
ld_tilenum = tile_coord.tile_id(1:ld_points);
ld_area = tile_coord.area(1:ld_points);

% read GL coordinate data
test_file = ['/discover/nobackup/elee15/GEOSldas_4_0/sims/IDS/fcst_ctrl/' exp_name '/output/SMAP_EASEv2_M36/cat/ens0000/Y' num2str(start_year) '/M0' num2str(start_month) '/' exp_name '.tavg24_1d_lnd_Nt.' num2str(start_year) '0' num2str(start_month) '01_1200z.nc4'];

lat_data = ncread(test_file,'lat');
lon_data = ncread(test_file,'lon');
lat_data = double(lat_data);
lon_data = double(lon_data);

nearest_tile = NaN*ones(length(lat_data),1);
nearest_index = NaN*ones(length(lat_data),1);
nearest_lat = NaN*ones(length(lat_data),1);
nearest_lon = NaN*ones(length(lat_data),1);

for i = 1:length(index_cell)

    dist_vec = sqrt((ld_lats(index_cell{i}) - lat_data(i)).^2 + (ld_lons(index_cell{i})-lon_data(i)).^2);
    min_index = find(dist_vec==min(dist_vec));
    if length(min_index)>1
       tmp = randi([1 length(min_index)],1,1);
       min_index = min_index(tmp);
    end
    nearest_tile(i) = ld_tilenum(index_cell{i}(min_index));
    nearest_index(i) = index_cell{i}(min_index);
    nearest_lat(i) = ld_lats(index_cell{i}(min_index));
    nearest_lon(i) = ld_lons(index_cell{i}(min_index));
end

nearest_tile_cn = NaN*ones(length(lat_data),1);
nearest_index_cn = NaN*ones(length(lat_data),1);
nearest_lat_cn = NaN*ones(length(lat_data),1);
nearest_lon_cn = NaN*ones(length(lat_data),1);

for i = 1:length(index_cell)
    
    tmp_index = index_cell{i};
    tmp_lat = ld_lats(tmp_index);
    tmp_lon = ld_lons(tmp_index);
    tmp_tilenum = ld_tilenum(tmp_index);
    tmp_ity = fz_ity_data(tmp_index);
    data_index = tmp_ity == gl_ity_data(i);
    if sum(data_index)>0
    tmp_index(data_index==0) = [];
    tmp_lat(data_index==0) = [];
    tmp_lon(data_index==0) = [];
    tmp_tilenum(data_index==0) = [];
    end

    dist_vec = sqrt((tmp_lat - lat_data(i)).^2 + (tmp_lon-lon_data(i)).^2);
    min_index = find(dist_vec==min(dist_vec));
    if length(min_index)>1
       tmp = randi([1 length(min_index)],1,1);
       min_index = min_index(tmp);
    end
    nearest_tile_cn(i) = tmp_tilenum(min_index);
    nearest_index_cn(i) = tmp_index(min_index);
    nearest_lat_cn(i) = tmp_lat(min_index);
    nearest_lon_cn(i) = tmp_lon(min_index);
end



mapped_rzexc_data = fz_rzexc_data(nearest_index);
mapped_catdef_data = fz_catdef_data(nearest_index); 
mapped_ghtcnt1_data = fz_ghtcnt1_data(nearest_index); 
mapped_htsnnn1_data = fz_htsnnn1_data(nearest_index); 
mapped_vgmax_data = fz_vgmax_data(nearest_index); 

mapped_tc_data = fz_tc_data(nearest_index_cn);
mapped_cnpft1_data = fz_cnpft1_data(nearest_index_cn); 
mapped_cnpft2_data = fz_cnpft2_data(nearest_index_cn);  
mapped_cnpft3_data = fz_cnpft3_data(nearest_index_cn); 

test_lat = unique(lat_data);
test_lon = unique(lon_data);

[lon_mat,lat_mat] = meshgrid(test_lon,test_lat);

mapped_rzexc_mat = NaN*ones(size(lon_mat));
gl_rzexc_mat = NaN*ones(size(lon_mat));
mapped_catdef_mat = NaN*ones(size(lon_mat));                                              
gl_catdef_mat = NaN*ones(size(lon_mat)); 
mapped_ghtcnt1_mat = NaN*ones(size(lon_mat));                                              
gl_ghtcnt1_mat = NaN*ones(size(lon_mat)); 
mapped_tc_mat = NaN*ones(size(lon_mat)); 
gl_tc_mat = NaN*ones(size(lon_mat)); 
mapped_htsnnn1_mat = NaN*ones(size(lon_mat));                                                 
gl_htsnnn1_mat = NaN*ones(size(lon_mat)); 
mapped_cnpft1_mat = NaN*ones(size(lon_mat));                                                 
gl_cnpft1_mat = NaN*ones(size(lon_mat));
mapped_cnpft2_mat = NaN*ones(size(lon_mat)); 
gl_cnpft2_mat = NaN*ones(size(lon_mat));
mapped_cnpft3_mat = NaN*ones(size(lon_mat)); 
gl_cnpft3_mat = NaN*ones(size(lon_mat));
mapped_vgmax_mat = NaN*ones(size(lon_mat));                                              
gl_vgmax_mat = NaN*ones(size(lon_mat)); 

for i = 1:length(lon_data)

    lat_index = test_lat == lat_data(i);
    lon_index = test_lon == lon_data(i);

    mapped_rzexc_mat(lat_index,lon_index) = mapped_rzexc_data(i);
    gl_rzexc_mat(lat_index,lon_index) = gl_rzexc_data(i);

    mapped_catdef_mat(lat_index,lon_index) = mapped_catdef_data(i);                        
    gl_catdef_mat(lat_index,lon_index) = gl_catdef_data(i);    

    mapped_ghtcnt1_mat(lat_index,lon_index) = mapped_ghtcnt1_data(i);                        
    gl_ghtcnt1_mat(lat_index,lon_index) = gl_ghtcnt1_data(i); 

    mapped_tc_mat(lat_index,lon_index) = mapped_tc_data(i);                        
    gl_tc_mat(lat_index,lon_index) = gl_tc_data(i);     

    mapped_htsnnn1_mat(lat_index,lon_index) = mapped_htsnnn1_data(i);                              
    gl_htsnnn1_mat(lat_index,lon_index) = gl_htsnnn1_data(i);

    mapped_cnpft1_mat(lat_index,lon_index) = mapped_cnpft1_data(i);                                  
    gl_cnpft1_mat(lat_index,lon_index) = gl_cnpft1_data(i); 

    mapped_cnpft2_mat(lat_index,lon_index) = mapped_cnpft2_data(i);  
    gl_cnpft2_mat(lat_index,lon_index) = gl_cnpft2_data(i); 

    mapped_cnpft3_mat(lat_index,lon_index) = mapped_cnpft3_data(i);  
    gl_cnpft3_mat(lat_index,lon_index) = gl_cnpft3_data(i);   

    mapped_vgmax_mat(lat_index,lon_index) = mapped_vgmax_data(i);                       
    gl_vgmax_mat(lat_index,lon_index) = gl_vgmax_data(i); 
end

rzexc_diff_mat = mapped_rzexc_mat - gl_rzexc_mat;
catdef_diff_mat = mapped_catdef_mat - gl_catdef_mat; 
ghtcnt1_diff_mat = mapped_ghtcnt1_mat - gl_ghtcnt1_mat; 
tc_diff_mat = mapped_tc_mat - gl_tc_mat; 
htsnnn1_diff_mat = mapped_htsnnn1_mat - gl_htsnnn1_mat;
cnpft1_diff_mat = mapped_cnpft1_mat - gl_cnpft1_mat;
cnpft2_diff_mat = mapped_cnpft2_mat - gl_cnpft2_mat;
cnpft3_diff_mat = mapped_cnpft3_mat - gl_cnpft3_mat;
vgmax_diff_mat = mapped_vgmax_mat - gl_vgmax_mat;  

load coastlines;
[latcells, loncells] = polysplit(coastlat, coastlon);


latlim = [-60 90];
lonlim = [-180 180];

cmap_mat = cmocean('balance',12);
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
s1 = quantile(rzexc_diff_mat(:),0.05);
s2 = quantile(rzexc_diff_mat(:),0.95);
caxis_min = -1*max(abs(s1),abs(s2));
caxis_max = max(abs(s1),abs(s2));
figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off; gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...
 'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...
    'MLabelParallel','south','FontSize',12)
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,rzexc_diff_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);                                                                  
tightmap                                                                                 
title('total diff mean RZEXC mapped-36')                                     
h = colorbar;                                                                            
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);                  
set(gcf,'paperunits','centimeters')                                                      
set(gcf,'paperposition',[1 1 17 15])                                                     
saveas(gcf,['/home/jkolassa/Projects/Eunjee_rst_comp/RZEXC_rst_regrid'],'png') 

cmap_mat = cmocean('balance',12);                                                        
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                         
cmap = cmap_mat;                                                                         
s1 = quantile(catdef_diff_mat(:),0.05);                                                   
s2 = quantile(catdef_diff_mat(:),0.95);                                                   
caxis_min = -1*max(abs(s1),abs(s2));                                                     
caxis_max = max(abs(s1),abs(s2));                                                        
figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                  
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);                                        
framem off; gridm on;                                                                    
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                   
 'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                        
    'MLabelParallel','south','FontSize',12)                                              
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])                                     
pcolorm(lat_mat,lon_mat,catdef_diff_mat)                                                  
caxis([caxis_min caxis_max]);                                                            
set(gcf,'Colormap',cmap);                                                                
plotm(coastlat, coastlon,'k')                                                            
hold on
set(gca,'XLim',lonlim);                                                                  
set(gca,'YLim',latlim);                                                                  
tightmap                                                                                 
title('total diff mean CATDEF mapped-36')                                                 
h = colorbar;                                                                            
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);                  
set(gcf,'paperunits','centimeters')                                                      
set(gcf,'paperposition',[1 1 17 15])                                                     
saveas(gcf,['/home/jkolassa/Projects/Eunjee_rst_comp/CATDEF_rst_regrid'],'png') 

cmap_mat = cmocean('balance',12);                                                        
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                         
cmap = cmap_mat;                                                                         
s1 = quantile(ghtcnt1_diff_mat(:),0.05);                                                   
s2 = quantile(ghtcnt1_diff_mat(:),0.95);                                                   
caxis_min = -1*max(abs(s1),abs(s2));                                                     
caxis_max = max(abs(s1),abs(s2));                                                        
figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                  
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);                                        
framem off; gridm on;                                                                    
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                   
 'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                        
    'MLabelParallel','south','FontSize',12)                                              
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])                                     
pcolorm(lat_mat,lon_mat,ghtcnt1_diff_mat)                                                  
caxis([caxis_min caxis_max]);                                                            
set(gcf,'Colormap',cmap);                                                                
plotm(coastlat, coastlon,'k')                                                            
hold on
set(gca,'XLim',lonlim);                                                                  
set(gca,'YLim',latlim);                                                                  
tightmap                                                                                 
title('total diff mean GHTCNT1 mapped-36')                                                 
h = colorbar;                                                                            
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);                  
set(gcf,'paperunits','centimeters')                                                      
set(gcf,'paperposition',[1 1 17 15])                                                     
saveas(gcf,['/home/jkolassa/Projects/Eunjee_rst_comp/GHTCNT1_rst_regrid'],'png')

cmap_mat = cmocean('balance',12);                                                        
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                         
cmap = cmap_mat;                                                                         
s1 = quantile(tc_diff_mat(:),0.05);                                                   
s2 = quantile(tc_diff_mat(:),0.95);                                                   
caxis_min = -1*max(abs(s1),abs(s2));                                                     
caxis_max = max(abs(s1),abs(s2));                                                        
figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                  
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);                                        
framem off; gridm on;                                                                    
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                   
 'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                        
    'MLabelParallel','south','FontSize',12)                                              
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])                                     
pcolorm(lat_mat,lon_mat,tc_diff_mat)                                                  
caxis([caxis_min caxis_max]);                                                            
set(gcf,'Colormap',cmap);                                                                
plotm(coastlat, coastlon,'k')                                                            
hold on
set(gca,'XLim',lonlim);                                                                  
set(gca,'YLim',latlim);                                                                  
tightmap                                                                                 
title('total diff mean TC mapped-36')                                                 
h = colorbar;                                                                            
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);                  
set(gcf,'paperunits','centimeters')                                                      
set(gcf,'paperposition',[1 1 17 15])                                                     
saveas(gcf,['/home/jkolassa/Projects/Eunjee_rst_comp/TC_rst_regrid'],'png')

cmap_mat = cmocean('balance',12);                                                        
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                         
cmap = cmap_mat;                                                                         
s1 = quantile(htsnnn1_diff_mat(:),0.05);                                                      
s2 = quantile(htsnnn1_diff_mat(:),0.95);                                                      
caxis_min = -1*max(abs(s1),abs(s2));                                                     
caxis_max = max(abs(s1),abs(s2));                                                        
figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                  
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);                                        
framem off; gridm on;                                                                    
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                   
 'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                        
    'MLabelParallel','south','FontSize',12)                                              
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])                                     
pcolorm(lat_mat,lon_mat,htsnnn1_diff_mat)                                                     
caxis([caxis_min caxis_max]);                                                            
set(gcf,'Colormap',cmap);                                                                
plotm(coastlat, coastlon,'k')                                                            
hold on
set(gca,'XLim',lonlim);                                                                  
set(gca,'YLim',latlim);                                                                  
tightmap                                                                                 
title('total diff mean HTSNNN1 mapped-36')                                                    
h = colorbar;                                                                            
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);                  
set(gcf,'paperunits','centimeters')                                                      
set(gcf,'paperposition',[1 1 17 15])                                                     
saveas(gcf,['/home/jkolassa/Projects/Eunjee_rst_comp/HTSNNN1_rst_regrid'],'png') 

cmap_mat = cmocean('balance',12);                                                        
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                         
cmap = cmap_mat;                                                                         
s1 = quantile(cnpft1_diff_mat(:),0.05);                                                               
s2 = quantile(cnpft1_diff_mat(:),0.95);                                                               
caxis_min = -1*max(abs(s1),abs(s2));                                                     
caxis_max = max(abs(s1),abs(s2));                                                        
figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                  
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);                                        
framem off; gridm on;                                                                    
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                   
 'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                        
    'MLabelParallel','south','FontSize',12)                                              
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])                                     
pcolorm(lat_mat,lon_mat,cnpft1_diff_mat)                                                              
caxis([caxis_min caxis_max]);                                                            
set(gcf,'Colormap',cmap);                                                                
plotm(coastlat, coastlon,'k')                                                            
hold on
set(gca,'XLim',lonlim);                                                                  
set(gca,'YLim',latlim);                                                                  
tightmap                                                                                 
title('total diff mean CNPFT1 mapped-36')                                                             
h = colorbar;                                                                            
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);                  
set(gcf,'paperunits','centimeters')                                                      
set(gcf,'paperposition',[1 1 17 15])                                                     
saveas(gcf,['/home/jkolassa/Projects/Eunjee_rst_comp/CNPFT1_rst_regrid'],'png')

cmap_mat = cmocean('balance',12);                                                        
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                         
cmap = cmap_mat;
s1 = quantile(cnpft2_diff_mat(:),0.05); 
s2 = quantile(cnpft2_diff_mat(:),0.95); 
caxis_min = -1*max(abs(s1),abs(s2));
caxis_max = max(abs(s1),abs(s2));
figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                  
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off; gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                   
 'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                        
    'MLabelParallel','south','FontSize',12)                                              
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,cnpft2_diff_mat) 
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('total diff mean CNPFT1 mapped-36')                                                                  
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);                  
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/home/jkolassa/Projects/Eunjee_rst_comp/CNPFT2_rst_regrid'],'png')

cmap_mat = cmocean('balance',12);                                                        
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                         
cmap = cmap_mat;
s1 = quantile(cnpft3_diff_mat(:),0.05); 
s2 = quantile(cnpft3_diff_mat(:),0.95); 
caxis_min = -1*max(abs(s1),abs(s2));
caxis_max = max(abs(s1),abs(s2));
figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],...                  
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
framem off; gridm on;
setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                   
 'ParallelLabel','on','MLabelLocation',60,'PLabelLocation',20,...                        
    'MLabelParallel','south','FontSize',12)                                              
geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
pcolorm(lat_mat,lon_mat,cnpft3_diff_mat) 
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('total diff mean CNPFT3 mapped-36')                                                                  
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);                  
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/home/jkolassa/Projects/Eunjee_rst_comp/CNPFT3_rst_regrid'],'png')                              
