clc
clear all
close all

s2s_forcing_path = '/discover/nobackup/projects/gmao/geos_carb/fzeng/FORECASTS/GEOS5/CLIM/GEOS5v2/BCSD_Final/6-hourly/2007/jan01/ens1/';
s2s_forcing_file = [s2s_forcing_path 'GEOS5.200701.nc4'];

tlml_data = ncread(s2s_forcing_file,'TLML');
mean_tlml_data = nanmean(tlml_data,3);
swgdn_data = ncread(s2s_forcing_file,'SWGDN');
mean_swgdn_data = nanmean(swgdn_data,3);
wind_data = ncread(s2s_forcing_file,'SPEED');
mean_wind_data = nanmean(wind_data,3);
lat_data = ncread(s2s_forcing_file,'latitude');
lat_data = double(lat_data);
lon_data = ncread(s2s_forcing_file,'longitude');
lon_data = double(lon_data);

exp_forcing_path = '/discover/nobackup/jkolassa/Eunjee_test/IDS_fcst_CTRL_M09_2007_ens1_timeFix/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/Y2007/M01/';
exp_forcing_file = [exp_forcing_path 'IDS_fcst_CTRL_M09_2007_ens1_timeFix.tavg24_1d_lfs_Nt.monthly.200701.nc4'];

tlml_lfs_data = ncread(exp_forcing_file,'Tair');
%mean_tlml_lfs_data = nanmean(tlml_lfs_data,3);
swgdn_lfs_data = ncread(exp_forcing_file,'SWdown');
%mean_swgdn_lfs_data = nanmean(swgdn_lfs_data,3);
wind_lfs_data = ncread(exp_forcing_file,'Wind');
%mean_wind_lfs_data = nanmean(wind_lfs_data,3);
rainf_lfs_data = ncread(exp_forcing_file,'Rainf');
rainfc_lfs_data = ncread(exp_forcing_file,'Rainf_C');
qair_lfs_data = ncread(exp_forcing_file,'Qair');
psurf_lfs_data = ncread(exp_forcing_file,'Psurf');
snowf_lfs_data = ncread(exp_forcing_file,'Snowf');
lwdown_lfs_data = ncread(exp_forcing_file,'LWdown');
swnet_lfs_data = ncread(exp_forcing_file,'SWnet');

exp_data_file = [exp_forcing_path 'IDS_fcst_CTRL_M09_2007_ens1_timeFix.tavg24_1d_lnd_Nt.monthly.200701.nc4'];
gpp_data = ncread(exp_data_file,'CNGPP');
gpp_data = gpp_data*1000*86400;

exp_tile_coord_file = '/discover/nobackup/jkolassa/Eunjee_test/IDS_fcst_CTRL_M09_2007_ens1_timeFix/output/SMAP_EASEv2_M09_GLOBAL/rc_out/IDS_fcst_CTRL_M09_2007_ens1_timeFix.ldas_tilecoord.bin';
exp_tile_coord = read_tilecoord_le(exp_tile_coord_file);

% mapping data to tile space 

dlon = 0.625;
dlat = 0.5;
mapped_mean_swgdn = NaN*ones(exp_tile_coord.N_tile,1);
for i = 1:exp_tile_coord.N_tile

    lon_dist = sqrt((lon_data-exp_tile_coord.com_lon(i)).^2);
    i_ind  = find(lon_dist==min(lon_dist));  
    lat_dist = sqrt((lat_data-exp_tile_coord.com_lat(i)).^2);                     
    j_ind = find(lat_dist==min(lat_dist)); 

    if (i_ind>length(lon_data))  
        i_ind=1                                                                                           
    end 

    tmp = mean_swgdn_data(i_ind,j_ind);
    mapped_mean_swgdn(i) = nanmean(tmp(:));
end  
%for i = 1:length(gpp_vec)
%    if (mod(i,10000)==0)
%        i
%    end
%    point_index = find(exp_tile_coord.tile_id==tile_coord.tile_id(i));
%    gpp_diff_vec(i) = gpp_vec(i) - gpp_data(point_index);
%    swdown_diff_vec(i) = swdown_vec(i) - swgdn_lfs_data(point_index);
%    ld2gl_map(i) = point_index;
%end
%
%save('ld2gl_map_index.mat','ld2gl_map');

test_lat = unique(exp_tile_coord.com_lat);
test_lon = unique(exp_tile_coord.com_lon);

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


exp_swdown_mat = NaN*ones(size(lon_mat));
forc_swgdn_mat = NaN*ones(size(lon_mat));

for i = 1:length(exp_tile_coord.com_lon)

    lat_index = test_lat == exp_tile_coord.com_lat(i);
    lon_index = test_lon == exp_tile_coord.com_lon(i);

     exp_swdown_mat(lat_index,lon_index) = swgdn_lfs_data(i);
     forc_swgdn_mat(lat_index,lon_index) = mapped_mean_swgdn(i); 
end

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

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

cmap_mat = cmocean('solar',20);
%cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
cmap_mat = flipud(cmap_mat)
cmap = cmap_mat;
s1 = quantile(exp_swdown_mat(:),0.95);
s2 = quantile(exp_swdown_mat(:),0.05);
caxis_min = s2;
caxis_max = s1;

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,exp_swdown_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('SWdown Jan [W/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/Eunjee_test/SWdown_GEOSldas_Jan_2007'],'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(double(lat_data),double(lon_data),mean_swgdn_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('SWdown Jan [W/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/Eunjee_test/SWdown_forcing_Jan_2007'],'png')

swdown_diff_mat = forc_swgdn_mat - exp_swdown_mat;

cmap_mat = cmocean('balance',20);                                                                              
cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);                                            
cmap_mat = flipud(cmap_mat)                                                                                  
cmap = cmap_mat;                                                                                             
s1 = quantile(swdown_diff_mat(:),0.95);                                                                       
s2 = quantile(swdown_diff_mat(:),0.05);                                                                       
caxis_min = s2;                                                                                              
caxis_max = s1; 

x = 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,swdown_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('\Delta SWdown Jan [W/m^2] (Hindcast - GEOSldas) / (in - out)')
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/Eunjee_test/SWdown_diff_forc_GEOSldas'],'png')
