clc
clear all
close all

nt_in_rst = '/discover/nobackup/jkolassa/NT_CNCLM45/GEOSldas_CNCLM45_NT_rst/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CNCLM45_NT_rst.catchcnclm45_internal_rst.19800101_0000';
ref_in_rst = '/discover/nobackup/jkolassa/CN40_CN45_test_runs/GEOSldas_CN45/input/restart/catchcnclm45_internal_rst';
rst_gen_in_rst = '/discover/nobackup/jkolassa/NT_CNCLM45/GEOSldas_CNCLM45_NT_rst/input/restart/catchcnclm45_internal_rst';
matt_nt_in_rst = '/discover/nobackup/mathomp4/SystemTests/runs/LDAS_GLOBALCNCLM45/model/BASELINE/input/restart/catchcnclm45_internal_rst';

nt_in_info = ncinfo(nt_in_rst);
ref_in_info = ncinfo(ref_in_rst);
rst_gen_in_info = ncinfo(rst_gen_in_rst);
matt_nt_in_info = ncinfo(matt_nt_in_rst);

lnd_file = '/discover/nobackup/jkolassa/CN40_CN45_test_runs/GEOSldas_CN45/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/GEOSldas_CN45.tavg24_1d_lnd_Nt.monthly.198001.nc4';

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

file_1 = '/discover/nobackup/jkolassa/CN40_CN45_test_runs/GEOSldas_CN45/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CN45.catchcnclm45_internal_rst.19800101_0000';
file_2 = '/discover/nobackup/jkolassa/NT_CNCLM45/GEOSldas_CNCLM45_NT_rst/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CNCLM45_NT_rst.catchcnclm45_internal_rst.19800101_0000';
file_3 = '/discover/nobackup/mathomp4/SystemTests/runs/LDAS_GLOBALCNCLM45/model/CURRENT/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/CURRENT.tavg24_1d_lnd_Nt.19800101_1200z.nc4';
%file_3 = '/discover/nobackup/jkolassa/development_dir/GEOSldas_CNCLM45_NT_rst_v17_10_0/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CNCLM45_NT_rst_v17_10_0.catchcnclm45_internal_rst.19800101_0000';
%file_4 = '/discover/nobackup/jkolassa/development_dir/undo_fix/GEOSldas_CNCLM45_NT_rst_undo/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CNCLM45_NT_rst_undo.catchcnclm45_internal_rst.19800101_0000';
%file_5 = '/discover/nobackup/jkolassa/development_dir/develop/GEOSldas_CN45_develop_LT_bcs/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CN45_develop_LT_bcs.catchcnclm45_internal_rst.19800101_0000';
%file_6 = '/discover/nobackup/jkolassa/development_dir/develop/GEOSldas_CN45_develop_BO_bcs/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CN45_develop_BO_bcs.catchcnclm45_internal_rst.19800101_0000';
%file_7 = '/discover/nobackup/jkolassa/development_dir/new_fix/GEOSldas_CNCLM45_NT_rst/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CNCLM45_NT_rst.catchcnclm45_internal_rst.19800101_0000';
orig_file = '/discover/nobackup/jkolassa/clm4.5_restart/clm4.5_DE720/output/global/rs/ens0000/Y2017/M01/clm4.5_DE720.ens0000.catchcn_ldas_rst.20170101_0000z';
test_file_1 = '/discover/nobackup/jkolassa/e0004s_transientCO2_05_rst/e0004s_transientCO2/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y1981/M01/e0004s_transientCO2.ens0000.catchcn_ldas_rst.19810101_0000z';
test_file_2 = '/discover/nobackup/jkolassa/e0004s_transientCO2_05_rst/e0004s_transientCO2/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y1981/M01/e0004s_transientCO2.ens0000.catchcn40_ldas_rst.19810101_0000z';

file_1_info = ncinfo(file_1);
file_2_info = ncinfo(file_2);

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

tile_id_1 = double(ncread(file_1,'TILE_ID'));
tile_id_orig = double(ncread(orig_file,'TILE_ID'));
tile_id_test_1 = double(ncread(test_file_1,'TILE_ID'));
tile_id_test_2 = double(ncread(test_file_2,'TILE_ID'));

cnpft_1 = double(ncread(file_1,'CNPFT'));
cnpft_2 = double(ncread(file_2,'CNPFT'));
%cnpft_3 = ncread(file_3,'CNPFT');
%cnpft_4 = ncread(file_4,'CNPFT');
%cnpft_5 = ncread(file_5,'CNPFT');
%cnpft_6 = ncread(file_6,'CNPFT');
%cnpft_7 = ncread(file_7,'CNPFT');

cncol_1 = double(ncread(file_1,'CNCOL'));
cncol_2 = double(ncread(file_2,'CNCOL'));
%cncol_3 = ncread(file_3,'CNCOL');
%cncol_4 = ncread(file_4,'CNCOL');
%cncol_5 = ncread(file_5,'CNCOL');
%cncol_6 = ncread(file_6,'CNCOL');
%cncol_7 = ncread(file_7,'CNCOL');

index = 0;
elai_vec_1 = NaN*ones(length(srfexc_1)*3*4,1);
elai_vec_2 = NaN*ones(length(srfexc_1)*3*4,1);
%elai_vec_3 = NaN*ones(length(srfexc_1)*3*4,1);
%elai_vec_4 = NaN*ones(length(srfexc_1)*3*4,1);
%elai_vec_5 = NaN*ones(length(srfexc_1)*3*4,1);
%elai_vec_6 = NaN*ones(length(srfexc_1)*3*4,1);
%elai_vec_7 = NaN*ones(length(srfexc_1)*3*4,1);

for i = 1:length(srfexc_1)

    for nz = 1:3
        for nv = 1:4
            index = index + 1;
            elai_vec_1(index) = cnpft_1(i,nz*nv*(71-1));
            elai_vec_2(index) = cnpft_2(i,nz*nv*(71-1));
 %           elai_vec_3(index) = cnpft_3(i,nz*nv*(71-1));
 %           elai_vec_4(index) = cnpft_4(i,nz*nv*(71-1));
 %           elai_vec_5(index) = cnpft_5(i,nz*nv*(71-1));
 %           elai_vec_6(index) = cnpft_6(i,nz*nv*(71-1));
 %           elai_vec_7(index) = cnpft_7(i,nz*nv*(71-1));
        end
    end
end

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

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

elai_1_mat = NaN*ones(size(lon_mat));
srfexc_1_mat = NaN*ones(size(lon_mat));
%srfexc_3_mat = NaN*ones(size(lon_mat));
%srfexc_4_mat = NaN*ones(size(lon_mat));
%srfexc_5_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);

    elai_1_mat(lat_index,lon_index) = elai_vec_1(i);
    srfexc_1_mat(lat_index,lon_index) = srfexc_1(i);
    %srfexc_3_mat(lat_index,lon_index) = srfexc_3(i);
    %srfexc_4_mat(lat_index,lon_index) = srfexc_4(i);
    %srfexc_5_mat(lat_index,lon_index) = srfexc_5(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;
caxis_min = -0.5;
caxis_max = 0.5;

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_CNCLM45_1'],'png')

science_file_1 = '/gpfsm/dnb31/jkolassa/NT_CNCLM45/GEOSldas_CNCLM45_NT_rst/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/GEOSldas_CNCLM45_NT_rst.tavg24_1d_lnd_Nt.monthly.198001.nc4';
science_file_2 = '/discover/nobackup/jkolassa/development_dir/v17.10.0/GEOSldas_CN45_v17_10/output/SMAP_EASEv2_M36/cat/ens0000/Y1980/M01/GEOSldas_CN45_v17_10.tavg24_1d_lnd_Nt.monthly.198001.nc4';

%science_file_info_1 = ncinfo(science_file_1);
%science_file_info_2 = ncinfo(science_file_2);

%cnclm40_rst_develop = '/discover/nobackup/jkolassa/development_dir/develop/GEOSldas_CN40_develop/output/SMAP_EASEv2_M09_GLOBAL/rs/ens0000/Y1981/M01/GEOSldas_CN40_develop.catchcnclm40_internal_rst.19810101_0000';
%cnclm40_rst_v17 = '
