clc
clear all
close all

old_file = '/discover/nobackup/projects/gmao/geos_carb/elee15/IDS/GEOSldas40_sims/CTRL/IDS_fcst_CTRL_M09_2001_ens1/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2001/M01/IDS_fcst_CTRL_M09_2001_ens1.tavg24_1d_lfs_Nt.20010101_1200z.nc4';
new_file = '/discover/nobackup/jkolassa/staging_dir/IDS_fcst_CTRL_M09_2001_ens1_test/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2001/M01/IDS_fcst_CTRL_M09_2001_ens1_test.tavg24_1d_lfs_Nt.20010101_1200z.nc4';
exp_2_file = '/discover/nobackup/projects/gmao/geos_carb/elee15/IDS/GEOSldas40_sims/CTRL/IDS_fcst_CTRL_M09_2001_ens1_RcZero//output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2001/M01/IDS_fcst_CTRL_M09_2001_ens1_RcZero.tavg24_1d_lfslnd_fcst_Nt.20010101_1200z.nc4';
old_data_file = '/discover/nobackup/projects/gmao/geos_carb/elee15/IDS/GEOSldas40_sims/CTRL/IDS_fcst_CTRL_M09_2001_ens1/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2001/M01/IDS_fcst_CTRL_M09_2001_ens1.tavg24_1d_lfslnd_fcst_Nt.20010101_1200z.nc4';
new_data_file = '/discover/nobackup/jkolassa/staging_dir/IDS_fcst_CTRL_M09_2001_ens1_test/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2001/M01/IDS_fcst_CTRL_M09_2001_ens1_test.tavg24_1d_lnd_Nt.20010101_1200z.nc4';
exp_2_lfs_file = '/discover/nobackup/projects/gmao/geos_carb/elee15/IDS/GEOSldas40_sims/CTRL/IDS_fcst_CTRL_M09_2001_ens1_RcZero//output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2001/M01/IDS_fcst_CTRL_M09_2001_ens1_RcZero.tavg24_1d_lfs_Nt.20010101_1200z.nc4';


sfmc_old = double(ncread(old_data_file,'SFMC'));
sfmc_new = double(ncread(new_data_file,'SFMC'));
sfmc_exp_2 = double(ncread(exp_2_file,'SFMC'));

rzmc_old = double(ncread(old_data_file,'RZMC'));
rzmc_new = double(ncread(new_data_file,'RZMC'));

runoff_old = double(ncread(old_data_file,'RUNOFF'))*86400;
runoff_new = double(ncread(new_data_file,'RUNOFF'))*86400;

evland_old = double(ncread(old_data_file,'EVLAND'))*86400;
evland_new = double(ncread(new_data_file,'EVLAND'))*86400;

rainf_old = double(ncread(old_file,'Rainf'))*86400;
snowf_old = double(ncread(old_file,'Snowf'))*86400;
rainfc_old = double(ncread(old_file,'Rainf_C'))*86400;

rainf_new = double(ncread(new_file,'Rainf'))*86400;
snowf_new = double(ncread(new_file,'Snowf'))*86400;
rainfc_new = double(ncread(new_file,'Rainf_C'))*86400;

rainf_exp_2 = double(ncread(exp_2_lfs_file,'Rainf'))*86400;
rainfc_exp_2 = double(ncread(exp_2_lfs_file,'Rainf_C'))*86400;
snowf_exp_2 = double(ncread(exp_2_lfs_file,'Snowf'))*86400;

gpp_old = double(ncread(old_data_file,'CNGPP'))*86400*1000;
gpp_new = double(ncread(new_data_file,'CNGPP'))*86400*1000;
gpp_exp_2 = double(ncread(exp_2_file,'CNGPP'))*86400*1000;

tot_precip_old = rainf_old + snowf_old;
tot_precip_new = rainf_new + snowf_new;

diff_tot_precip = tot_precip_new - tot_precip_old;
diff_rainf = rainf_new - rainf_old;
diff_snowf = snowf_new - snowf_old;
diff_rainfc = rainfc_new - rainfc_old;
diff_rainfc_exp_2_old = rainfc_exp_2 - rainfc_old;
diff_rainfc_new_exp_2 = rainfc_new - rainfc_exp_2;

diff_gpp_new_old = gpp_new-gpp_old;
diff_gpp_new_exp_2 = gpp_new-gpp_exp_2;
diff_gpp_exp_2_old = gpp_exp_2-gpp_old;


diff_sfmc = sfmc_new - sfmc_old;
diff_runoff = runoff_new - runoff_old;
diff_evland = evland_new - evland_old;
diff_rzmc = rzmc_new - rzmc_old;

lat_vec = double(ncread(old_file,'lat'));
lon_vec = double(ncread(old_file,'lon'));

test_lat = unique(lat_vec);
test_lon = unique(lon_vec);

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

diff_tot_precip_mat = NaN*ones(size(lon_mat));
diff_rainf_mat = NaN*ones(size(lon_mat));
diff_snowf_mat = NaN*ones(size(lon_mat));
diff_sfmc_mat = NaN*ones(size(lon_mat));
diff_runoff_mat = NaN*ones(size(lon_mat));
diff_evland_mat = NaN*ones(size(lon_mat));
diff_rzmc_mat = NaN*ones(size(lon_mat));
diff_gpp_new_old_mat = NaN*ones(size(lon_mat));
diff_gpp_new_exp_2_mat = NaN*ones(size(lon_mat));
diff_gpp_exp_2_old_mat = NaN*ones(size(lon_mat));
diff_rainfc_mat = NaN*ones(size(lon_mat));  
diff_rainfc_exp_2_old_mat = NaN*ones(size(lon_mat));  
diff_rainfc_new_exp_2_mat = NaN*ones(size(lon_mat));  

for i = 1:length(lon_vec)

    lat_index = test_lat == lat_vec(i);
    lon_index = test_lon == lon_vec(i);

    diff_tot_precip_mat(lat_index,lon_index) = diff_tot_precip(i);
    diff_rainf_mat(lat_index,lon_index) = diff_rainf(i);
    diff_snowf_mat(lat_index,lon_index) = diff_snowf(i);
    diff_sfmc_mat(lat_index,lon_index) = diff_sfmc(i);
    diff_runoff_mat(lat_index,lon_index) = diff_runoff(i);
    diff_evland_mat(lat_index,lon_index) = diff_evland(i);
    diff_rzmc_mat(lat_index,lon_index) = diff_rzmc(i);
    diff_gpp_new_old_mat(lat_index,lon_index) = diff_gpp_new_old(i);
    diff_gpp_new_exp_2_mat(lat_index,lon_index) = diff_gpp_new_exp_2(i);
    diff_gpp_exp_2_old_mat(lat_index,lon_index) = diff_gpp_exp_2_old(i);
    diff_rainfc_mat(lat_index,lon_index) = diff_rainfc(i); 
    diff_rainfc_exp_2_old_mat(lat_index,lon_index) = diff_rainfc_exp_2_old(i);
    diff_rainfc_new_exp_2_mat(lat_index,lon_index) = diff_rainfc_new_exp_2(i);
end

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;
caxis_min = -5%min(diff_tot_precip);
caxis_max = 5%max(diff_tot_precip);

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,diff_tot_precip_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('total precipitation (Rainf + Snowf) difference [mm/day] (new - old) Jan-01')
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/lfs_tot_precip_diff'],'png')

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,diff_rainf_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('total liquid precipitation (Rainf) difference [mm/day] (new - old) Jan-01')
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/lfs_rainf_diff'],'png')

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,diff_snowf_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('total solid precipitation (Snowf) difference [mm/day] (new - old) Jan-01')
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/lfs_snowf_diff'],'png')

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,diff_rainfc_mat)                                                                                     
caxis([caxis_min caxis_max]);                                                                                                    
set(gcf,'Colormap',cmap);                                                                                                        
plotm(coastlat, coastlon,'k')                                                                                                    
hold on                                                                                                                          
%geoshow(states,'facealpha',0);                                                                                                  
set(gca,'XLim',lonlim);                                                                                                          
set(gca,'YLim',latlim);                                                                                                          
tightmap                                                                                                                         
title('convective precipitation (Rainf_C) difference [mm/day] (new - old) Jan-01')                                              
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/diff_rainfc_new_old'],'png')

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,diff_rainfc_new_exp_2_mat)    
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                                                                                  
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('convective precipitation (Rainf_C) difference [mm/day] (new - exp_2) Jan-01')        
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/diff_rainfc_new_exp_2'],'png')    

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,diff_rainfc_exp_2_old_mat)    
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                                                                                  
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('convective precipitation (Rainf_C) difference [mm/day] (exp_2 - old) Jan-01')        
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/diff_rainfc_exp_2_old'],'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(diff_sfmc,0.05);
s2 = quantile(diff_sfmc,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,diff_sfmc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('SFMC difference [m^3 m^{-3}] (new - old) Jan-01')
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/diff_sfmc'],'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(diff_rzmc,0.05);
s2 = quantile(diff_rzmc,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,diff_rzmc_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('RZMC difference [m^3 m^{-3}] (new - old) Jan-01')
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/diff_rzmc'],'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(diff_runoff,0.05);
s2 = quantile(diff_runoff,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,diff_runoff_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('runoff difference [mm/day] (new - old) Jan-01')
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/diff_runoff'],'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(diff_evland,0.05);
s2 = quantile(diff_evland,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,diff_evland_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('EVLAND difference [mm/day] (new - old) Jan-01')
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/diff_evland'],'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(diff_gpp_new_old(:),0.05);
s2 = quantile(diff_gpp_new_old(:),0.95);
caxis_min = -1;%-1*max(abs(s1),abs(s2));
caxis_max = 1;%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,diff_gpp_new_old_mat)                                                                            
caxis([caxis_min caxis_max]);                                                                                           
set(gcf,'Colormap',cmap);                                                                                               
plotm(coastlat, coastlon,'k')                                                                                           
hold on                                                                                                                 
%geoshow(states,'facealpha',0);                                                                                         
set(gca,'XLim',lonlim);                                                                                                 
set(gca,'YLim',latlim);                                                                                                 
tightmap                                                                                                                
title('\Delta GPP [gC/m^2/day] (new - old) Jan-01')                                     
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/diff_gpp_new_old'],'png') 

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,diff_gpp_new_exp_2_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                                                                         
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('\Delta GPP [gC/m^2/day] (new - exp_2) Jan-01')                                     
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/diff_gpp_new_exp_2'],'png')

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,diff_gpp_exp_2_old_mat)
caxis([caxis_min caxis_max]);
set(gcf,'Colormap',cmap);
plotm(coastlat, coastlon,'k')
hold on
%geoshow(states,'facealpha',0);                                                                                         
set(gca,'XLim',lonlim);
set(gca,'YLim',latlim);
tightmap
title('\Delta GPP [gC/m^2/day] (exp_2 - old) Jan-01')                                     
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/diff_gpp_exp_2_old'],'png')
