clc
clear all
close all

pft2_file = '/discover/nobackup/jkolassa/restarts/GEOSldas_CNCLM45_NT/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CNCLM45_NT.catchcnclm51_internal_rst.19800101_0000';
pft4_file = '/discover/nobackup/jkolassa/restarts/GEOSldas_CNCLM45_NT/output/SMAP_EASEv2_M36/rs/ens0000/Y1980/M01/GEOSldas_CNCLM45_NT.catchcnclm51_internal_rst.19800101_0000.4PFT';
latlon_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_vec = double(ncread(latlon_file,'lat'));
lon_vec = double(ncread(latlon_file,'lon'));

tile_id = double(ncread(pft2_file,'TILE_ID'));

cnpft_2 = double(ncread(pft2_file,'CNPFT'));
cnpft_4 = double(ncread(pft4_file,'CNPFT'));

ityp_2 = double(ncread(pft2_file,'ITY'));
ityp_4 = double(ncread(pft4_file,'ITY'));

fveg_2 = double(ncread(pft2_file,'FVG'));
fveg_4 = double(ncread(pft4_file,'FVG'));

var_index = 69;

elai_2 = NaN*ones(length(lat_vec),3,2);
elai_4 = NaN*ones(length(lat_vec),3,4);

for nz = 1:3
    for nv = 1:2

        new_index = (var_index-1)*6 + (nv-1)*3 + nz;
        nv_old = (nv-1)*2 + [1 2];
        old_index = (var_index-1)*12 + (nv_old-1)*3 + nz;

        elai_2(:,nz,nv) = cnpft_2(:,new_index);
        elai_4(:,nz,nv_old) = cnpft_4(:,old_index);
 
  end 
end

test = nanmean(elai_4(:,1,1:2),3);
diff_vec = test - elai_2(:,1,1);

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

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

elai_2_plot_mat_1 = NaN*ones(size(lon_mat,1),size(lon_mat,2));
elai_4_plot_mat_1 = NaN*ones(size(lon_mat,1),size(lon_mat,2));
elai_4_plot_mat_2 = NaN*ones(size(lon_mat,1),size(lon_mat,2));

for i = 1:length(lon_vec)

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

    elai_2_plot_mat_1(lat_index,lon_index,:) = elai_2(i,1,1);
    elai_4_plot_mat_1(lat_index,lon_index,:) = elai_4(i,1,1);
    elai_4_plot_mat_2(lat_index,lon_index,:) = elai_4(i,1,2);
end

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

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

cmap_mat = cmocean('algae',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                                   
cmap = cmap_mat;
s1 = quantile(elai_2_plot_mat_1(:),0.01);
s2 = quantile(elai_2_plot_mat_1(:),0.99);
caxis_min = s1;
caxis_max = s2;

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,elai_2_plot_mat_1)
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(['ELAI [-] '])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['test_elai2_1'],'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,elai_4_plot_mat_1)
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(['ELAI [-] '])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['test_elai4_1'],'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,elai_4_plot_mat_2)
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(['ELAI [-] '])
h = colorbar;
set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,['test_elai4_2'],'png')
