clc
clear all
close all

test_file = '/discover/nobackup/jkolassa/test_dir/GEOSldas_CN40/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y1981/M01/GEOSldas_CN40.tavg24_1d_lnd_Nt.19810101_1200z.nc4';
ref_file = '/discover/nobackup/jkolassa/CN40_CN45_test_runs/GEOSldas_CN40/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y1981/M01/GEOSldas_CN40.tavg24_1d_lnd_Nt.19810101_1200z.nc4';

test_sfmc = double(ncread(test_file,'SFMC'));
ref_sfmc = double(ncread(ref_file,'SFMC'));

diff_sfmc = test_sfmc-ref_sfmc;

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

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

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

diff_sfmc_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);

    diff_sfmc_mat(lat_index,lon_index) = diff_sfmc(i);
end

disp('plot data')

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(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('mean SFMC dfference (exp-ref) [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/SFMC_diff'],'png')
