clc 
clear all
close all

orig_file = '/discover/nobackup/elee15/GEOSldas_4_0/restarts_link/SMAP_EASEv2_M09/e0004s_transientCO2/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y2001/M04/e0004s_transientCO2.catchcnclm40_internal_rst.20010401_0000';

mod_file = '/discover/nobackup/elee15/GEOSldas_4_0/restarts_link/rst_2013sm_snow/SMAP_EASEv2_M09/e0004s_transientCO2/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y2001/M04/e0004s_transientCO2.catchcnclm40_internal_rst.20010401_0000';

orig_srfexc = double(ncread(orig_file,'SRFEXC'));
mod_srfexc = double(ncread(mod_file,'SRFEXC'));

diff_srfexc = mod_srfexc - orig_srfexc;

figure 
plot(orig_srfexc,'-b','LineWidth',2)
hold on
plot(mod_srfexc,'-r','LineWidth',2)
hold off
legend('original','modified')
saveas(gcf,'eunjee_rst_test_SRFEXC.png')

test_file = '/gpfsm/dnb31/jkolassa/CN40_CN45_test_runs/GEOSldas_CN40/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y1981/M01/GEOSldas_CN40.tavg24_1d_lnd_Nt.19810131_1200z.nc4';

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

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

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

orig_srfexc_mat = NaN*ones(size(lon_mat));
mod_srfexc_mat = NaN*ones(size(lon_mat));

for i = 1:length(lon_vec)

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

    orig_srfexc_mat(lat_index,lon_index) = orig_srfexc(i);
    mod_srfexc_mat(lat_index,lon_index) = mod_srfexc(i);
end 

% plot data

disp('plot data')


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;
caxis_min = quantile(orig_srfexc(:),0.05);
caxis_max = quantile(orig_srfexc(:),0.95);

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,orig_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 [-] original')
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['eunjee_rst_SRFEXC_orig'],'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,mod_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 [-] modified')
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['eunjee_rst_SRFEXC_mod'],'png')

