clc
clear all
close all

read_path = '/discover/nobackup/jkolassa/test_dir/GEOSldas_CN51/output/SMAP_EASEv2_M36/cat/ens0000/';
clm45_read_path = '/discover/nobackup/projects/geoscm/fzeng/GEOSldas45_M36_rst_ldas06_15/output/SMAP_EASEv2_M36/cat/ens0000/';

start_year = 1980;
stop_year = 1988;
start_month = 1;
stop_month = 4;

leapYears = 1980:4:2022;
monthLength_nly = [31 28 31 30 31 30 31 31 30 31 30 31];
monthLength_ly = [31 29 31 30 31 30 31 31 30 31 30 31];

% read in PFT data

rst_file = '/discover/nobackup/jkolassa/test_dir/GEOSldas_CN51/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CN51.catchcnclm51_internal_rst.19800101_0000';
tile_id = ncread(rst_file,'TILE_ID');
ityp = ncread(rst_file,'ITY');
fvg = ncread (rst_file,'FVG');
    
for i = 1:length(ityp)
    [max_value,max_index] = max(fvg(i,:)); 
    if (max_value>=0.7)
       dom_pft(i) = ityp(i,max_index);
    else
       dom_pft(i) = NaN;
    end 
end 

pft_name = {...
           'Needleleaf evergreen temperate tree',...
           'Needleleaf evergreen boreal tree',...
           'Needleleaf deciduous boreal tree',...
           'Broadleaf evergreen tropical tree',...
           'Broadleaf evergreen temperate tree',...
           'Broadleaf deciduous tropical tree',...
           'Broadleaf deciduous temperate tree',...
           'Broadleaf deciduous boreal tree',...
           'Broadleaf evergreen temperate shrub',...
           'Broadleaf deciduous temperate shrub',...
           'Broadleaf deciduous boreal shrub',...
           'Arctic C3 grass',...
           'Non-arctic C3 grass',...
           'C4 grass',...
           'C3 crop'...
           };


% read in model output

mat_index = 0;

for y = start_year:stop_year

    y_str = num2str(y,'%04i');

    % define month length based on leap year status
    if sum(y==leapYears)>0
        monthLength = monthLength_ly;
    else
        monthLength = monthLength_nly;
    end

    % loop over months
    for m = 1:12

        mat_index = mat_index + 1;

        m_str = num2str(m,'%02i');

        read_file = [read_path 'Y' y_str '/M' m_str '/GEOSldas_CN51.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];
        clm45_read_file = [clm45_read_path 'Y' y_str '/M' m_str '/GEOSldas45_M36_rst_ldas06_15.tavg24_1d_lnd_Nt.monthly.' y_str m_str '.nc4'];


        if (mat_index==1)
           lat_vec = double(ncread(read_file,'lat'));
           lon_vec = double(ncread(read_file,'lon'));

           cntlai_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           cngpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);

           clm45_cntlai_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
           clm45_cngpp_mat = NaN*ones(length(lat_vec),(stop_year-start_year+1)*12);
        end

        datenum_vec(mat_index,1) = datenum(y,m,1);
        year_vec(mat_index) = y;
        month_vec(mat_index) = m;

        cntlai_mat(:,mat_index) = double(ncread(read_file,'CNTLAI'));
        cngpp_mat(:,mat_index) = double(ncread(read_file,'CNGPP'));

        clm45_cntlai_mat(:,mat_index) = double(ncread(clm45_read_file,'CNTLAI'));
        clm45_cngpp_mat(:,mat_index) = double(ncread(clm45_read_file,'CNGPP'));

         
    end
end 

cntlai_mat(:,mat_index+1:end) = [];
cngpp_mat(:,mat_index+1:end) = [];
clm45_cntlai_mat(:,mat_index+1:end) = [];
clm45_cngpp_mat(:,mat_index+1:end) = [];



%nan_index = (cntlai_mat==0) + (lat_vec>66);
%nan_index = nan_index == 2;
cntlai_mat(cntlai_mat==0) = NaN;
cngpp_mat(cngpp_mat==0) = NaN;
clm45_cntlai_mat(cntlai_mat==0) = NaN;
clm45_cngpp_mat(cngpp_mat==0) = NaN;

% compute and plot annual global averages

ann_index = 0;

for y = start_year:stop_year

    ann_index = ann_index + 1;
    year_index = (year_vec==y);

    tmp = cngpp_mat(:,year_index);
    mean_global_ann_cngpp(ann_index) = nanmean(tmp(:));
    mean_local_ann_cngpp(:,ann_index) = nanmean(tmp,2);

    tmp = clm45_cngpp_mat(:,year_index);
    mean_global_ann_clm45_cngpp(ann_index) = nanmean(tmp(:));
    mean_local_ann_clm45_cngpp(:,ann_index) = nanmean(tmp,2);

    for p = 1:15

       pft_index = find(dom_pft==p);

       pft_index_cell{p} = pft_index;

       tmp = cngpp_mat(pft_index,year_index);
       mean_global_ann_cngpp_pft(p,ann_index) = nanmean(tmp(:));

       tmp = clm45_cngpp_mat(pft_index,year_index);
       mean_global_ann_clm45_cngpp_pft(p,ann_index) = nanmean(tmp(:));
    end
end 

figure
plot(start_year:stop_year,mean_global_ann_cngpp,'LineWidth',2);
hold on
plot(start_year:stop_year,mean_global_ann_clm45_cngpp,'LineWidth',2);
hold off
grid on
xlabel('year')
ylabel('GPP [kg m^{-2} s^{-1}]')
title('annual average GPP [kg m^{-2} s^{-1}]')
legend('CN51','CN45','Location','bestoutside')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/GPP_ann_avg'],'png')

% plot PFT-specific annual averages

for p = 1:15
     if length(pft_index_cell{p})>1
       figure
       plot(start_year:stop_year,mean_global_ann_cngpp_pft(p,:),'LineWidth',2);
       hold on
       plot(start_year:stop_year,mean_global_ann_clm45_cngpp_pft(p,:),...
       'LineStyle','--','LineWidth',2);
       hold off
       grid on
       xlabel('year')
       ylabel('GPP [kgC m^{-2} s^{-1}]')
       legend('CN51','CN45')
       title([pft_name{p} ' (' num2str(length(pft_index_cell{p})) ')'])
       set(gcf,'paperunits','centimeters')
       set(gcf,'paperposition',[1 1 17 15])
       saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/GPP_ann_avg_pft' num2str(p)],'png')
     end
end

% make global map plots

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

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

cntlai_plot_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2),size(cntlai_mat,2));
cngpp_plot_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2),size(cngpp_mat,2));
clm45_cntlai_plot_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2),size(clm45_cntlai_mat,2));
clm45_cngpp_plot_mat = NaN*ones(size(lon_mat,1),size(lon_mat,2),size(clm45_cngpp_mat,2));

for i = 1:length(lon_vec)

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

    cntlai_plot_mat(lat_index,lon_index,:) = cntlai_mat(i,:);
    cngpp_plot_mat(lat_index,lon_index,:) = cngpp_mat(i,:);
    clm45_cntlai_plot_mat(lat_index,lon_index,:) = clm45_cntlai_mat(i,:);
    clm45_cngpp_plot_mat(lat_index,lon_index,:) = clm45_cngpp_mat(i,:);

end

[mod_aug, B]        = geotiffread('BlueMarbleNG-TB_2004-08-01_rgb_3600x1800.TIFF');

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

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

cmap_mat = cmocean('algae',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;
s1 = quantile(cntlai_plot_mat(:),0.01);                                                                             
s2 = quantile(cntlai_plot_mat(:),0.99);                                                                             
caxis_min = s1;
caxis_max = s2;

caxis_min = s1;
caxis_max = s2;

for i = 1:size(cntlai_plot_mat,3)

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,cntlai_plot_mat(:,:,i)) 
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 LAI [-] ' datestr(datenum_vec(i),'yyyy-mm')])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CNTLAI_' datestr(datenum_vec(i),'yyyy-mm')],'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,clm45_cntlai_plot_mat(:,:,i))
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 LAI [-] ' datestr(datenum_vec(i),'yyyy-mm')])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CLM45_CNTLAI_' datestr(datenum_vec(i),'yyyy-mm')],'png')


%ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 179],...                                             
%    'FLatLimit',[-60 90],'FLonLimit',[-180 179]);
%%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])
%geoshow(mod_aug, B, 'facealpha', 1)
%pcolorm(lat_mat,lon_mat,cntlai_plot_mat(:,:,i))
%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 LAI [-] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
%h = colorbar;
%set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
%set(gcf,'paperunits','centimeters')
%set(gcf,'paperposition',[1 1 17 15])
%saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CNTLAI_MODIS_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')
%
%figure
%ax = axesm('robinson','MapLatLimit',[-60 90],'MapLonLimit',[-180 179],...                                             
%    'FLatLimit',[-60 90],'FLonLimit',[-180 179]);
%%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])
%geoshow(mod_aug, B, 'facealpha', 0.7)
%pcolorm(lat_mat,lon_mat,cntlai_plot_mat(:,:,i))
%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 LAI [-] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
%h = colorbar;
%set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
%set(gcf,'paperunits','centimeters')
%set(gcf,'paperposition',[1 1 17 15])
%saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CNTLAI_MODIS_robinson_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')

end

cntlai_avg_plot_mat = nanmean(cntlai_plot_mat,3);

figure
ax = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 179],...                                             
    'FLatLimit',[-60 90],'FLonLimit',[-180 179]);
%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])
geoshow(mod_aug, B, 'facealpha', 1)
pcolorm(lat_mat,lon_mat,cntlai_avg_plot_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(['annual average LAI [-] ' datestr(datenum_vec(i),'yyyy')])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CNTLAI_MODIS_ann_avg_' datestr(datenum_vec(i),'yyyy')],'png')

figure
ax = axesm('robinson','MapLatLimit',[-60 90],'MapLonLimit',[-180 179],...                                             
    'FLatLimit',[-60 90],'FLonLimit',[-180 179]);
%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])
geoshow(mod_aug, B, 'facealpha', 1)
pcolorm(lat_mat,lon_mat,cntlai_avg_plot_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(['annual average LAI [-] ' datestr(datenum_vec(i),'yyyy')])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CNTLAI_MODIS_ann_avg_robinson_' datestr(datenum_vec(i),'yyyy')],'png')



cmap_mat = cmocean('algae',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;
s1 = quantile(cngpp_plot_mat(:),0.01);
s2 = quantile(cngpp_plot_mat(:),0.99);
caxis_min = s1;
caxis_max = s2;

caxis_min = s1;
caxis_max = s2;

for i = 1:size(cngpp_plot_mat,3)

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,cngpp_plot_mat(:,:,i))
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(['GPP [kgC m^{-2} s^{-1}] ' datestr(datenum_vec(i),'yyyy-mm')])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CNGPP_' datestr(datenum_vec(i),'yyyy-mm')],'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,clm45_cngpp_plot_mat(:,:,i))
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(['GPP [kgC m^{-2} s^{-1}] ' datestr(datenum_vec(i),'yyyy-mm')])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CLM45_CNGPP_' datestr(datenum_vec(i),'yyyy-mm')],'png')

end

cngpp_diff_plot_mat = clm45_cngpp_plot_mat - cngpp_plot_mat;

cmap_mat = cmocean('balance',12);
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;
s1 = quantile(cngpp_diff_plot_mat(:),0.01);
s2 = quantile(cngpp_diff_plot_mat(:),0.99);
caxis_min = -1*(max(abs(s1),abs(s2)));
caxis_max = (max(abs(s1),abs(s2)));

for i = 1:size(cngpp_diff_plot_mat,3)

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,cngpp_diff_plot_mat(:,:,i))
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 (CN45 - CN51) [kgC m^{-2} s^{-1}] ' datestr(datenum_vec(i),'yyyy-mm')])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CNGPP_diff_' datestr(datenum_vec(i),'yyyy-mm')],'png')

end

cntlai_diff_plot_mat = clm45_cntlai_plot_mat - cntlai_plot_mat;

cmap_mat = cmocean('balance',12);
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
s1 = quantile(cntlai_diff_plot_mat(:),0.01);
s2 = quantile(cntlai_diff_plot_mat(:),0.99);
caxis_min = -1*(max(abs(s1),abs(s2)));
caxis_max = (max(abs(s1),abs(s2)));

for i = 1:size(cntlai_diff_plot_mat,3)

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,cntlai_diff_plot_mat(:,:,i))
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 LAI (CN45 - CN51) [kgC m^{-2} s^{-1}] ' datestr(datenum_vec(i),'yyyy-mm')])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/CNTLAI_diff_' datestr(datenum_vec(i),'yyyy-mm')],'png')

end


figure
ax = axesm('miller','MapLatLimit',[-60 85],'MapLonLimit',[-180 179],... 
    'FLatLimit',[-60 85],'FLonLimit',[-180 179]);
% setm(ax,'MLineLocation',20,'PLineLocation',20,'MeridianLabel','on',...                                                      
%     'ParallelLabel','on','MLabelLocation',20,'PLabelLocation',20,...                                                      
%     'MLabelParallel','south','FontSize',14)

geoshow(mod_aug, B, 'facealpha', 0.7)
saveas(gcf,['/discover/nobackup/jkolassa/CN51_plots/global/test_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')


