clc
clear all
close all

lnd_file = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/GEOSldas_1981_present/GEOSldas_CNCLM40_9km/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y1981/M05/GEOSldas_CNCLM40_9km.tavg24_1d_lnd_Nt.monthly.198105.nc4';
ref_lnd_file = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/202001_to_present/GEOSldas_e0004s_transientCO2_05/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2021/M11/GEOSldas_e0004s_transientCO2_05.tavg24_1d_lnd_Nt.monthly.202111.nc4';


rst_file = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/GEOSldas_1981_present/GEOSldas_CNCLM40_9km/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y1981/M01/GEOSldas_CNCLM40_9km.catchcnclm40_internal_rst.19810101_0000';
tsurf_rst = double(ncread(rst_file,'TSURF'));
srfexc_rst = double(ncread(rst_file,'SRFEXC'));
rst_info = ncinfo(rst_file);

rst_file_2 = '/discover/nobackup/projects/geoscm/fzeng/Catchment-CN40_9km/202001_to_present/GEOSldas_e0004s_transientCO2_05/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y2021/M08/GEOSldas_e0004s_transientCO2_05.catchcnclm40_internal_rst.20210801_0000';
%tsurf_rst_2 = double(ncread(rst_file_2,'TSURF'));
srfexc_rst_2 = double(ncread(rst_file_2,'SRFEXC'));
rst_2_info = ncinfo(rst_file_2);

test_rst = '/discover/nobackup/jkolassa/development_dir/v17.10.0/GEOSldas_CNCLM40_9km/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y1981/M01/GEOSldas_CNCLM40_9km.catchcnclm40_internal_rst.19810101_0000';
test_rst_info = ncinfo(test_rst);
srfexc_rst_test = double(ncread(test_rst,'SRFEXC'));
test_tsurf = ncread(test_rst,'TSURF');

ref_rst = '/discover/nobackup/jkolassa/CN40_CN45_test_runs/GEOSldas_CN40/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y1981/M01/GEOSldas_CN40.catchcnclm40_internal_rst.19810101_0000';
ref_rst_info = ncinfo(ref_rst);
ref_tsurf = ncread(ref_rst,'TSURF');
srfexc_rst_ref = double(ncread(ref_rst,'SRFEXC'));

lnd_info = ncinfo(lnd_file);

lat_data = double(ncread(lnd_file,'lat'));
lon_data = double(ncread(lnd_file,'lon'));
scuns_data = double(ncread(lnd_file,'SCUNS'));
ref_scuns_data = double(ncread(ref_lnd_file,'SCUNS'));
lai_data = double(ncread(lnd_file,'CNLAI'));
sfmc_data = double(ncread(lnd_file,'SFMC'));

%ccn_4_file = '/gpfsm/dnb44/ltakacs/bcs/Heracles-NL/SMAP_EASEv2_M09/clsm/CLM_veg_typs_fracs';
%ccn_4_data = dlmread(ccn_4_file);
%ccn_4_tile_id = ccn_4_data(:,1);
%ccn_4_veg_type = ccn_4_data(:,3:6);
%ccn_4_veg_fracs = ccn_4_data(:,7:10);

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

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

scuns_mat = NaN*ones(size(lon_mat));
ref_scuns_mat = NaN*ones(size(lon_mat));
lai_mat = NaN*ones(size(lon_mat));
sfmc_mat = NaN*ones(size(lon_mat));
tsurf_mat = NaN*ones(size(lon_mat));
tsurf_ref_mat = NaN*ones(size(lon_mat));
tsurf_test_mat = NaN*ones(size(lon_mat));
srfexc_mat = NaN*ones(size(lon_mat));
srfexc_ref_mat = NaN*ones(size(lon_mat));
srfexc_test_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);

    scuns_mat(lat_index,lon_index) = scuns_data(i);
    ref_scuns_mat(lat_index,lon_index) = ref_scuns_data(i);
    lai_mat(lat_index,lon_index) = lai_data(i);
    sfmc_mat(lat_index,lon_index) = sfmc_data(i);
    tsurf_mat(lat_index,lon_index) = tsurf_rst(i);
    tsurf_ref_mat(lat_index,lon_index) = ref_tsurf(i);
    tsurf_test_mat(lat_index,lon_index) = test_tsurf(i);
    srfexc_mat(lat_index,lon_index) = srfexc_rst(i);
    srfexc_ref_mat(lat_index,lon_index) = srfexc_rst_ref(i);
    srfexc_test_mat(lat_index,lon_index) = srfexc_rst_test(i);
    srfexc_2_mat(lat_index,lon_index) = srfexc_rst_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(lai_mat(:),0.05);
s2 = quantile(lai_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,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('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/NT_comp/SFMC_Fanwei_total'],'png')
saveas(gcf,['/home/jkolassa/Projects/LAI'],'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(scuns_mat(:),0.05);
s2 = quantile(scuns_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,scuns_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('SCUNS [m s^-1]')
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/SCUNS'],'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(scuns_mat(:),0.05);
s2 = quantile(scuns_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,ref_scuns_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('SCUNS [m s^-1]')
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/ref_SCUNS'],'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(sfmc_mat(:),0.05);
s2 = quantile(sfmc_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,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('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/NT_comp/SFMC_Fanwei_total'],'png')
saveas(gcf,['/home/jkolassa/Projects/SFMC'],'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(tsurf_mat(:),0.05);
s2 = quantile(tsurf_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,tsurf_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('TSURF [K]')
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/TSURF'],'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(tsurf_ref_mat(:),0.05);
s2 = quantile(tsurf_ref_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,tsurf_ref_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('TSURF [-]')
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/TSURF_ref'],'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(tsurf_test_mat(:),0.05);
s2 = quantile(tsurf_test_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,tsurf_test_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('TSURF [-]')
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/TSURF_test'],'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_test_mat(:),0.05);
s2 = quantile(srfexc_test_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_test_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_test'],'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_mat(:),0.05);
s2 = quantile(srfexc_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_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'],'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_ref_mat(:),0.05);
s2 = quantile(srfexc_ref_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_ref_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_ref'],'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'],'png')
