clc 
clear all
close all

s2s_forcing_path = '/discover/nobackup/projects/gmao/geos_carb/fzeng/FORECASTS/GEOS5/CLIM/GEOS5v2/BCSD_Final/6-hourly/2007/jan01/ens1/';
s2s_forcing_file = [s2s_forcing_path 'GEOS5.200701.nc4'];

tlml_data = ncread(s2s_forcing_file,'TLML');
mean_tlml_data = nanmean(tlml_data,3);
swgdn_data = ncread(s2s_forcing_file,'SWGDN');
mean_swgdn_data = nanmean(swgdn_data,3);
wind_data = ncread(s2s_forcing_file,'SPEED');
mean_wind_data = nanmean(wind_data,3);
lat_data = ncread(s2s_forcing_file,'latitude');
lon_data = ncread(s2s_forcing_file,'longitude');

exp_forcing_path = '/discover/nobackup/jkolassa/Eunjee_test/IDS_fcst_CTRL_M09_2007_ens1_zthFilt/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2007/M01/';
exp_forcing_file = [exp_forcing_path 'IDS_fcst_CTRL_M09_2007_ens1_zthFilt.tavg24_1d_lfs_Nt.monthly.200701.nc4'];

tlml_lfs_data = ncread(exp_forcing_file,'Tair');
%mean_tlml_lfs_data = nanmean(tlml_lfs_data,3);
swgdn_lfs_data = ncread(exp_forcing_file,'SWdown');
%mean_swgdn_lfs_data = nanmean(swgdn_lfs_data,3);
wind_lfs_data = ncread(exp_forcing_file,'Wind');
%mean_wind_lfs_data = nanmean(wind_lfs_data,3);
rainf_lfs_data = ncread(exp_forcing_file,'Rainf');
rainfc_lfs_data = ncread(exp_forcing_file,'Rainf_C');
qair_lfs_data = ncread(exp_forcing_file,'Qair');
psurf_lfs_data = ncread(exp_forcing_file,'Psurf');
snowf_lfs_data = ncread(exp_forcing_file,'Snowf');
lwdown_lfs_data = ncread(exp_forcing_file,'LWdown');
swnet_lfs_data = ncread(exp_forcing_file,'SWnet');

exp_data_file = [exp_forcing_path 'IDS_fcst_CTRL_M09_2007_ens1_zthFilt.tavg24_1d_lnd_Nt.monthly.200701.nc4'];
gpp_data = ncread(exp_data_file,'CNGPP');
gpp_data = gpp_data*1000*86400;

exp_tile_coord_file = '/discover/nobackup/jkolassa/Eunjee_test/IDS_fcst_CTRL_M09_2007_ens1_zthFilt/output/SMAP_EASEv2_M09_GLOBAL/rc_out/IDS_fcst_CTRL_M09_2007_ens1_zthFilt.ldas_tilecoord.bin';
exp_tile_coord = read_tilecoord_le(exp_tile_coord_file);

tile_coord_file = '/discover/nobackup/projects/gmao/geos_carb/fzeng/Catchment/SMAP_EASEv2_M09/e0004s_transientCO2_BChindcast_2007Jan01_ens1/output/SMAP_EASEv2_M09_GLOBAL/rc_out/e0004s_transientCO2_BChindcast_2007Jan01.ldas_tilecoord.bin';
tile_coord = read_tilecoord( tile_coord_file);
N_tile = tile_coord.N_tile;
N_out_fields = 80;


out_file = '/discover/nobackup/projects/gmao/geos_carb/fzeng/Catchment/SMAP_EASEv2_M09/e0004s_transientCO2_BChindcast_2007Jan01_ens1/output/SMAP_EASEv2_M09_GLOBAL/cat/ens_avg/Y2007/M01/e0004s_transientCO2_BChindcast_2007Jan01.ens_avg.ldas_tile_monthly_out.200701.bin'
out_data = read_tile_data(out_file, N_tile, N_out_fields);

swdown_vec = out_data(8,:);
wind_vec = out_data(9,:);
tair_vec = out_data(1,:);
qair_vec = out_data(2,:);
psurf_vec = out_data(3,:);
rainfc_vec = out_data(4,:);
rainf_vec = out_data(5,:);
snowf_vec = out_data(6,:);
lwdown_vec = out_data(7,:);
gpp_vec = out_data(62,:);
gpp_vec = gpp_vec*86400;
gpp_diff_vec = NaN*ones(size(gpp_vec));
swdown_diff_vec = NaN*ones(size(gpp_vec));
ld2gl_map = NaN*ones(size(gpp_vec));

load('ld2gl_map_index.mat')
ld2gl_map = ld2gl_map';

tmp_gpp = gpp_data(ld2gl_map);
tmp_swgdn = swgdn_lfs_data(ld2gl_map);
tmp_wind = wind_lfs_data(ld2gl_map);
tmp_tair = tlml_lfs_data(ld2gl_map);
tmp_rainf = rainf_lfs_data(ld2gl_map)*86400;
tmp_rainfc = rainfc_lfs_data(ld2gl_map)*86400; 
tmp_snowf = snowf_lfs_data(ld2gl_map)*86400;   % Fanwei's data are in mm/d; GEOSldas outputd kg m^-2 s^-1
tmp_qair = qair_lfs_data(ld2gl_map);
tmp_lwdown = lwdown_lfs_data(ld2gl_map);
tmp_psurf = psurf_lfs_data(ld2gl_map);

gpp_diff_vec = tmp_gpp - gpp_vec';
swdown_diff_vec = tmp_swgdn - swdown_vec';
wind_diff_vec = tmp_wind - wind_vec';
tair_diff_vec = tmp_tair - tair_vec';
rainf_diff_vec = tmp_rainf - rainf_vec';
rainfc_diff_vec = tmp_rainfc - rainfc_vec';
snowf_diff_vec = tmp_snowf - snowf_vec';
qair_diff_vec = tmp_qair - qair_vec';
lwdown_diff_vec = tmp_lwdown - lwdown_vec';
psurf_diff_vec = tmp_psurf - psurf_vec';   

%for i = 1:length(gpp_vec)
%    if (mod(i,10000)==0)
%        i
%    end
%    point_index = find(exp_tile_coord.tile_id==tile_coord.tile_id(i));
%    gpp_diff_vec(i) = gpp_vec(i) - gpp_data(point_index);
%    swdown_diff_vec(i) = swdown_vec(i) - swgdn_lfs_data(point_index);
%    ld2gl_map(i) = point_index;
%end
%
%save('ld2gl_map_index.mat','ld2gl_map');

test_lat = unique(tile_coord.com_lat);
test_lon = unique(tile_coord.com_lon);

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

gpp_diff_mat = NaN*ones(size(lon_mat));
swdown_diff_mat = NaN*ones(size(lon_mat));
wind_diff_mat = NaN*ones(size(lon_mat)); 
tair_diff_mat = NaN*ones(size(lon_mat)); 
rainf_diff_mat = NaN*ones(size(lon_mat)); 
rainfc_diff_mat = NaN*ones(size(lon_mat)); 
snowf_diff_mat = NaN*ones(size(lon_mat)); 
qair_diff_mat = NaN*ones(size(lon_mat)); 
lwdown_diff_mat = NaN*ones(size(lon_mat)); 
psurf_diff_mat = NaN*ones(size(lon_mat)); 

for i = 1:length(tile_coord.com_lon)

    lat_index = test_lat == tile_coord.com_lat(i);
    lon_index = test_lon == tile_coord.com_lon(i);

     gpp_diff_mat(lat_index,lon_index) = gpp_diff_vec(i);
     swdown_diff_mat(lat_index,lon_index) = swdown_diff_vec(i);
     wind_diff_mat(lat_index,lon_index) = wind_diff_vec(i); 
     tair_diff_mat(lat_index,lon_index) = tair_diff_vec(i); 
     rainf_diff_mat(lat_index,lon_index) = rainf_diff_vec(i); 
     rainfc_diff_mat(lat_index,lon_index) = rainfc_diff_vec(i); 
     snowf_diff_mat(lat_index,lon_index) = snowf_diff_vec(i); 
     qair_diff_mat(lat_index,lon_index) = qair_diff_vec(i); 
     lwdown_diff_mat(lat_index,lon_index) = lwdown_diff_vec(i); 
     psurf_diff_mat(lat_index,lon_index) = psurf_diff_vec(i); 
end

load coastlines;                                                                                               
[latcells, loncells] = polysplit(coastlat, coastlon);                                                          
                                                                                        
latlim = [-60 90];                                                                                             
lonlim = [-180 180];                                                                                           

cmap_mat = cmocean('balance',20);                                                                              
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                               
cmap = cmap_mat;                                                                                               
s1 = quantile(gpp_diff_mat(:),0.95);                                                                     
s2 = quantile(gpp_diff_mat(:),0.05);                                                                     
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,gpp_diff_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 Jan [gC/m^2/day] (GEOSldas - LDASsa)')                                                  
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_test/ZF_test/GPP_diff_Jan_2007'],'png')

cmap_mat = cmocean('balance',20);                                                                                            
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;                                                                                                   
s1 = quantile(swdown_diff_mat(:),0.95);                                                                           
s2 = quantile(swdown_diff_mat(:),0.05);                                                                           
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,swdown_diff_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 SWdown Jan [W/m^2] (GEOSldas - LDASsa)')                                                       
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_test/ZF_test/SWdown_diff_Jan_2007'],'png')  

cmap_mat = cmocean('balance',20);                                                                                            
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;                                                                                                   
s1 = quantile(wind_diff_mat(:),0.95);                                                                           
s2 = quantile(wind_diff_mat(:),0.05);                                                                           
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,wind_diff_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 Wind [m/s] Jan (GEOSldas - LDASsa)')                                                         
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_test/ZF_test/Wind_diff_Jan_2007'],'png')  
 
cmap_mat = cmocean('balance',20);                                                                                            
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;                                                                                                   
s1 = quantile(rainf_diff_mat(:),0.95);                                                                           
s2 = quantile(rainf_diff_mat(:),0.05);                                                                           
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,rainf_diff_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 Rainf Jan [mm/day] (GEOSldas - LDASsa)')                                                         
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_test/ZF_test/Rainf_diff_Jan_2007'],'png')  
                   
cmap_mat = cmocean('balance',20);                                                                                            
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;                                                                                                   
s1 = quantile(rainfc_diff_mat(:),0.95);                                                                           
s2 = quantile(rainfc_diff_mat(:),0.05);                                                                           
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,rainfc_diff_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 RainfC Jan [mm/day] (GEOSldas - LDASsa)')                                                         
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_test/ZF_test/RainfC_diff_Jan_2007'],'png')  
         
cmap_mat = cmocean('balance',20);                                                                                            
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;                                                                                                   
s1 = quantile(snowf_diff_mat(:),0.95);                                                                           
s2 = quantile(snowf_diff_mat(:),0.05);                                                                           
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,snowf_diff_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 Snowf Jan [mm/day] (GEOSldas - LDASsa)')                                                         
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_test/ZF_test/Snowf_diff_Jan_2007'],'png')  

cmap_mat = cmocean('balance',20);                                                                                            
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;                                                                                                   
s1 = quantile(qair_diff_mat(:),0.95);                                                                           
s2 = quantile(qair_diff_mat(:),0.05);                                                                           
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,qair_diff_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 Qair Jan [kg/kg] (GEOSldas - LDASsa)')                                                         
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_test/ZF_test/Qair_diff_Jan_2007'],'png')  

cmap_mat = cmocean('balance',20);                                                                                            
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;                                                                                                   
s1 = quantile(lwdown_diff_mat(:),0.95);                                                                           
s2 = quantile(lwdown_diff_mat(:),0.05);                                                                           
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,lwdown_diff_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 LWdown [W/m^2] Jan (GEOSldas - LDASsa)')                                                         
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_test/ZF_test/LWdown_diff_Jan_2007'],'png')  

cmap_mat = cmocean('balance',20);
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
s1 = quantile(tair_diff_mat(:),0.95);
s2 = quantile(tair_diff_mat(:),0.05);
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,tair_diff_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 Tair [K] Jan (GEOSldas - LDASsa)')
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_test/ZF_test/Tair_diff_Jan_2007'],'png')

cmap_mat = cmocean('balance',20);
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
s1 = quantile(psurf_diff_mat(:),0.95);
s2 = quantile(psurf_diff_mat(:),0.05);
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,psurf_diff_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 Psurf [Pa] Jan  (GEOSldas - LDASsa)')
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_test/ZF_test/Psurf_diff_Jan_2007'],'png')
