clc
clear all
close all

NT_file = '/discover/nobackup/jkolassa/output/SMAP_EASEv2_M36/cat/ens0000/Y2015/M03/CURRENT.tavg24_1d_lnd_Nt.20150301_1200z.nc4';
Jana_file = '/discover/nobackup/jkolassa/GEOSgcm_GC_287_fixes_new_tag/GEOSldas_CNCLM40_NT_rst/output/SMAP_EASEv2_M36/cat/ens0000/Y2015/M03/GEOSldas_CNCLM40_NT_rst.tavg24_1d_lnd_Nt.20150301_1200z.nc4';
Fanwei_file = '/gpfsm/dnb05/projects/p144/fzeng/NOBACKUP/Catchment/SMAP_EASEv2_M09/e0004s_transientCO2_05/output/SMAP_EASEv2_M09_GLOBAL/cat/ens_avg/Y2015/M03/e0004s_transientCO2.ens_avg.ldas_tile_monthly_out.201503.bin';

NT_info = ncinfo(NT_file);
Jana_info = ncinfo(Jana_file);

lat_data = double(ncread(NT_file,'lat'));
lon_data =double( ncread(NT_file,'lon'));
lat_unique = unique(lat_data);
lon_unique = unique(lon_data);
[lon_mat,lat_mat] = meshgrid(lon_unique,lat_unique);


collectionID = 1;
nrOutVar = 44;
[data_tag, data_units] = get_data_tag(collectionID, nrOutVar, nrOutVar);
Fanwei_tile_file = '/gpfsm/dnb05/projects/p144/fzeng/NOBACKUP/Catchment/SMAP_EASEv2_M09/e0004s_transientCO2_05/output/SMAP_EASEv2_M09_GLOBAL/rc_out/e0004s_transientCO2.ldas_tilecoord.bin';
[tile_coord ] = read_tilecoord( Fanwei_tile_file, 0 );
ld_points = length(tile_coord.com_lat);
ld_lats = tile_coord.com_lat(1:ld_points);
ld_lons = tile_coord.com_lon(1:ld_points);
ld_tilenum = tile_coord.tile_id(1:ld_points);
ld_area = tile_coord.area(1:ld_points);

load coastlines;
[latcells, loncells] = polysplit(coastlat, coastlon);
latlim = [-60 90];
lonlim = [-180 180];

for i = 5:length(NT_info.Variables)

    var_name = NT_info.Variables(i).Name
    NT_vec = double(ncread(NT_file,var_name));
    Jana_vec = double(ncread(Jana_file,var_name));
  
    NT_data = NaN*ones(size(lon_mat));
    Jana_data = NaN*ones(size(lon_mat));

    for j = 1:length(lon_data)
 
        lat_index = lat_unique == lat_data(j);
        lon_index = lon_unique == lon_data(j);

        NT_data(lat_index,lon_index) = NT_vec(j);
        Jana_data(lat_index,lon_index) = NT_vec(j);
     end   

    figure 
    subplot(3,1,1)
    cmap_mat = cmocean('thermal',12);                                         
    %cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
    cmap = cmap_mat;          
    caxis_min = quantile(NT_data(:),0.05);                                                             
    caxis_max = quantile(NT_data(:),0.95); 

    ax1 = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],... 
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
    framem off; gridm on;
    setm(ax1,'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,NT_data)
    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([var_name ' NT'])
    h1 = colorbar(ax1);
    set(h1,'Location','eastoutside','Orientation','vertical','FontSize',14);
    freezeColors
    cbfreeze(h1,cmap)

    subplot(3,1,2)
    ax2 = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],... 
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
    framem off; gridm on; 
    setm(ax2,'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,Jana_data)
    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([var_name ' GEOSldas new tag'])
    h2 = colorbar(ax2);
    set(h2,'Location','eastoutside','Orientation','vertical','FontSize',14);
    freezeColors
    cbfreeze(h2,cmap)

    diff_mat = Jana_data-NT_data;
    cmap_mat = cmocean('balance',12);
    %cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
    cmap = cmap_mat;
    s1 = quantile(NT_data(:),0.05); 
    s2 = quantile(NT_data(:),0.95); 
    caxis_min = -1*max(abs(s1),abs(s2));
    caxis_max = max(abs(s1),abs(s2));

    subplot(3,1,3)
    ax3 = axesm('miller','MapLatLimit',[-60 90],'MapLonLimit',[-180 180],... 
    'FLatLimit',[-60 90],'FLonLimit',[-180 180]);
    framem off; gridm on; 
    setm(ax3,'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_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([var_name 'GEOSldas_{new} - NT'])
    h3 = colorbar(ax3);
    set(h3,'Location','eastoutside','Orientation','vertical','FontSize',14);
    freezeColors
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,['/home/jkolassa/Projects/NT_comp/' var_name],'png')  
    

end

