clc
clear all
close all

cn40_exp_dir = '/discover/nobackup/jkolassa/test_dir/GEOSldas_CN40/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y1981/M01/';
cn45_exp_dir = '/discover/nobackup/jkolassa/test_dir/GEOSldas_CN45/output/SMAP_EASEv2_M09/cat/ens0000/Y1980/M01/';

cn40_ref_dir = '/discover/nobackup/jkolassa/CN40_CN45_test_runs/GEOSldas_CN40/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y1981/M01/';
cn45_ref_dir = '/discover/nobackup/jkolassa/CN40_CN45_test_runs/GEOSldas_CN45/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/';

cn40_exp_file = [cn40_exp_dir 'GEOSldas_CN40.tavg24_1d_lnd_Nt.19810131_1200z.nc4'];
cn40_ref_file = [cn40_ref_dir 'GEOSldas_CN40.tavg24_1d_lnd_Nt.19810131_1200z.nc4'];

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

cn40_exp_sfmc_vec = double(ncread(cn40_exp_file,'SFMC'));
cn40_ref_sfmc_vec = double(ncread(cn40_ref_file,'SFMC'));

cn40_exp_lai_vec = double(ncread(cn40_exp_file,'CNLAI'));
cn40_ref_lai_vec = double(ncread(cn40_ref_file,'CNLAI'));

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

cn40_exp_sfmc_mat = NaN*ones(size(lon_mat));
cn40_ref_sfmc_mat = NaN*ones(size(lon_mat));
cn40_exp_lai_mat = NaN*ones(size(lon_mat));
cn40_ref_lai_mat = NaN*ones(size(lon_mat));

for i = 1:length(lat_vec)

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

    cn40_exp_sfmc_mat(lat_index,lon_index) = cn40_exp_sfmc_vec(i);
    cn40_ref_sfmc_mat(lat_index,lon_index) = cn40_ref_sfmc_vec(i);
    cn40_exp_lai_mat(lat_index,lon_index) = cn40_exp_lai_vec(i);
    cn40_ref_lai_mat(lat_index,lon_index) = cn40_ref_lai_vec(i);
end

cn40_diff_sfmc_mat = cn40_exp_sfmc_mat - cn40_ref_sfmc_mat;
cn40_diff_lai_mat = cn40_exp_lai_mat - cn40_ref_lai_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(cn40_diff_sfmc_mat(:),0.99);
s2 = quantile(cn40_diff_sfmc_mat(:),0.01);
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,cn40_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('\Delta SFMC [m^3 m^{-3}]')
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/CN40_news_bcs_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(cn40_diff_lai_mat(:),0.99);
s2 = quantile(cn40_diff_lai_mat(:),0.01);
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,cn40_diff_lai_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 LAI [-]')
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/CN40_new_bcs_LAI'],'png')
