clc
clear all
close all

exp_id = 'AGCM_CN45';

agcm_path = ['/discover/nobackup/jkolassa/CN40_CN45_test_runs/' exp_id '/scratch/'];

t_vec = {'0130','0430','0730','1030','1330','1630','1930','2230'};

if strcmp(exp_id,'AGCM_CN40')

   ref_file = '/discover/nobackup/jkolassa/test_dir/CN40_ref/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2016/M08/CN40_ref.tavg24_1d_lnd_Nt.monthly.201608.nc4';

   ref_lat = double(ncread(ref_file,'lat'));
   ref_lon = double(ncread(ref_file,'lon'));

   test_lat = unique(ref_lat);
   test_lon = unique(ref_lon);
   [ref_lon_mat,ref_lat_mat] = meshgrid(test_lon,test_lat);

   ref_lai_vec = double(ncread(ref_file,'CNTLAI'));
   ref_nee_vec = double(ncread(ref_file,'CNNEE'))*1000;
   ref_wet1_vec = double(ncread(ref_file,'GWETTOP'));
   ref_wet2_vec = double(ncread(ref_file,'GWETROOT'));
   ref_wat10cm_vec = double(ncread(ref_file,'WAT10CM'));
   ref_evpveg_vec = double(ncread(ref_file,'EVPTRNS'))*0.0864*0.408;

   ref_lai_data = NaN*ones(size(ref_lon_mat));
   ref_nee_data = NaN*ones(size(ref_lon_mat));
   ref_wet1_data = NaN*ones(size(ref_lon_mat));
   ref_wet2_data = NaN*ones(size(ref_lon_mat));
   ref_wat10cm_data = NaN*ones(size(ref_lon_mat));
   ref_evpveg_data = NaN*ones(size(ref_lon_mat));

   for i = 1:length(ref_lon)                                                              
                                                                                        
    lat_index = test_lat == ref_lat(i);                                                
    lon_index = test_lon == ref_lon(i);

    ref_lai_data(lat_index,lon_index) = ref_lai_vec(i);
    ref_nee_data(lat_index,lon_index) = ref_nee_vec(i);
    ref_wet1_data(lat_index,lon_index) = ref_wet1_vec(i);
    ref_wet2_data(lat_index,lon_index) = ref_wet2_vec(i);
    ref_wat10cm_data(lat_index,lon_index) = ref_wat10cm_vec(i);
    ref_evpveg_data(lat_index,lon_index) = ref_evpveg_vec(i); 

   end

elseif strcmp(exp_id,'AGCM_CN45')

   ref_lat = (-90+0.25):0.5:(90-0.25);
   ref_lon = (-180+0.25):0.5:(180-0.25);

   [ref_lon_mat,ref_lat_mat] = meshgrid(ref_lon,ref_lat);

   ref_file = '/discover/nobackup/elee15/GEOSldas_4_5/sims/GEOSldas45_M36_rst_ldas06_16/output/SMAP_EASEv2_M36/cat/ens0000/Y2016/M07/monthly_0.5deg_lnd_201607.mat';
   ref_data = load(ref_file);
   ref_lai_data = ref_data.regrid_mat(:,:,52);
   ref_lai_data = ref_lai_data';
   ref_nee_data = ref_data.regrid_mat(:,:,61)*1000;
   ref_nee_data = ref_nee_data';
   ref_wet1_data = ref_data.regrid_mat(:,:,3);
   ref_wet1_data = ref_wet1_data';
   ref_wet2_data = ref_data.regrid_mat(:,:,2);
   ref_wet2_data = ref_wet2_data';
   ref_wat10cm_data = ref_data.regrid_mat(:,:,75);
   ref_wat10cm_data = ref_wat10cm_data';
   ref_evpveg_data = ref_data.regrid_mat(:,:,25)*0.0864*0.408;
   ref_evpveg_data = ref_evpveg_data';
end

mat_index = 0;
disp('read data')
for d = 17:30

    d_str = num2str(d,'%02i');

    for t = 1:length(t_vec)

        mat_index = mat_index + 1;

        agcm_file = [agcm_path exp_id '.geosgcm_surf.201607' d_str '_' t_vec{t} 'z.nc4'];

        if (mat_index==1)
           lat_data = double(ncread(agcm_file,'lat'));
           lon_data = double(ncread(agcm_file,'lon'));
        end

        wet1_data(:,:,mat_index) = double(ncread(agcm_file,'WET1'));
        wet2_data(:,:,mat_index) = double(ncread(agcm_file,'WET2'));
        cnlai_data(:,:,mat_index) = double(ncread(agcm_file,'CNTLAI'));
        cnnee_data(:,:,mat_index) = double(ncread(agcm_file,'CNNEE'))*1000;
        wat10cm_data(:,:,mat_index) = double(ncread(agcm_file,'WAT10CM'));
        tsurf_data(:,:,mat_index) = double(ncread(agcm_file,'TSURF'));
        evpveg_data_tmp = double(ncread(agcm_file,'EVPVEG'));
        evpveg_data(:,:,mat_index) = evpveg_data_tmp*0.0864*0.408;
    end
end

disp('compute averages')

wet1_plot_data = nanmean(wet1_data,3);
wet2_plot_data = nanmean(wet2_data,3);
cnlai_plot_data = nanmean(cnlai_data,3);
cnnee_plot_data = nanmean(cnnee_data,3);
wat10cm_plot_data = nanmean(wat10cm_data,3);
tsurf_plot_data = nanmean(tsurf_data,3);
evpveg_plot_data = nanmean(evpveg_data,3);

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 = min(wet1_plot_data(:));
caxis_max = max(wet1_plot_data(:));

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_data,lon_data,wet1_plot_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('surface wetness [-]')
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/' exp_id '_wet1'],'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(ref_lat_mat,ref_lon_mat,ref_wet1_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('surface wetness [-]')
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/' exp_id '_wet1_ref'],'png')

cmap_mat = cmocean('rain',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
caxis_min = min(wet2_plot_data(:));
caxis_max = max(wet2_plot_data(:));

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_data,lon_data,wet2_plot_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('rootzone wetness [-]')
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/' exp_id '_wet2'],'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(ref_lat_mat,ref_lon_mat,ref_wet2_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('rootzone wetness [-]')
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/' exp_id '_wet2_ref'],'png')

cmap_mat = cmocean('rain',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
caxis_min = min(wat10cm_plot_data(:));
caxis_max = max(wat10cm_plot_data(:));

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_data,lon_data,wat10cm_plot_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('top 10 cm soil moisture [kg m^{-2}]')
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/' exp_id '_wat10cm'],'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(ref_lat_mat,ref_lon_mat,ref_wat10cm_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('top 10 cm soil moisture [kg m^{-2}]')
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/' exp_id '_wat10cm_ref'],'png')

cmap_mat = cmocean('speed',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
caxis_min = min(cnlai_plot_data(:));
caxis_max = 10%max(cnlai_plot_data(:));

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_data,lon_data,cnlai_plot_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('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/' exp_id '_cnlai'],'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(ref_lat_mat,ref_lon_mat,ref_lai_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('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/' exp_id '_cnlai_ref'],'png')

cmap_mat = cmocean('rain',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
caxis_min = min(cnnee_plot_data(:));
caxis_max = max(cnnee_plot_data(:));

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_data,lon_data,cnnee_plot_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('NEE [kg m^{-2} 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/' exp_id '_cnnee'],'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(ref_lon_mat,ref_lon_mat,ref_nee_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('NEE [kg m^{-2} 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/' exp_id '_cnnee_ref'],'png')

cmap_mat = cmocean('rain',12);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap = cmap_mat;
caxis_min = min(evpveg_plot_data(:));
caxis_max = max(evpveg_plot_data(:));

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_data,lon_data,evpveg_plot_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('transpiration [mm/day]')
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/' exp_id '_evpveg'],'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(ref_lat_mat,ref_lon_mat,ref_evpveg_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('transpiration [mm/day]')
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/' exp_id '_evpveg_ref'],'png')
