clc
clear all 
close all

lnd_file = '/discover/nobackup/jkolassa/CN_restart_debug/GEOSldas_CNCLM40_9km/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y1981/M01/GEOSldas_CNCLM40_9km.tavg24_1d_lnd_Nt.monthly.198101.nc4';
lat_data = double(ncread(lnd_file,'lat'));
lon_data = double(ncread(lnd_file,'lon'));

file_1 = '/discover/nobackup/jkolassa/CN_restart_debug/GEOSldas_CNCLM40_9km/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y1981/M01/GEOSldas_CNCLM40_9km.catchcnclm40_internal_rst.19810101_0000';
file_2 = '/discover/nobackup/jkolassa/CN_restart_debug/GEOSldas_CCN_comp_M09/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y1981/M01/GEOSldas_CCN_comp_M09.catchcnclm40_internal_rst.19810101_0000';

srfexc_1 = double(ncread(file_1,'SRFEXC'));
srfexc_2 = double(ncread(file_2,'SRFEXC'));

tile_id_1 = double(ncread(file_1,'TILE_ID'));
tile_id_2 = double(ncread(file_2,'TILE_ID'));

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

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

srfexc_1_mat = NaN*ones(size(lon_mat));
srfexc_2_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);

    srfexc_1_mat(lat_index,lon_index) = srfexc_1(i);
    srfexc_2_mat(lat_index,lon_index) = srfexc_2(i);

end

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


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

cmap_mat = cmocean('rain',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;
s1 = quantile(srfexc_1_mat(:),0.05);
s2 = quantile(srfexc_1_mat(:),0.95);
caxis_min = s1;
caxis_max = 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,srfexc_1_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('SRFEXC')
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/NT_comp/SFMC_Fanwei_total'],'png')
saveas(gcf,['/home/jkolassa/Projects/SRFEXC_1_debug'],'png')

cmap_mat = cmocean('rain',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;
s1 = quantile(srfexc_2_mat(:),0.05);
s2 = quantile(srfexc_2_mat(:),0.95);
caxis_min = s1;
caxis_max = 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,srfexc_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('SRFEXC')
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/NT_comp/SFMC_Fanwei_total'],'png')
saveas(gcf,['/home/jkolassa/Projects/SRFEXC_2_debug'],'png')


