clc
clear all
close all

% plot elevation used in FP and L4SM for region of TC Idai (Mozambique, Madagascar)

latrange = [-27 -9];
lonrange = [31  51];

elevrange = [0 2200];

g0 = 9.8;

% FP constants file

fp_fpath = '/gpfsm/dnb04/projects/p14/intermediate/f5271_fp/diag/Y2021/M05';
fp_fname = 'f5271_fp.const_2d_asm_Nx.20210520_0300z.nc4';

% FP TOPO file

fp_topo_path = '/gpfsm/dnb44/ltakacs/bcs/Icarus/NCAR_TOPO_GMTED/TOPO_CF0720x6C/smoothed/';
fp_topo_file = [fp_topo_path 'gmted_DYN_ave_720x4320.nc4'];

til_file = '/discover/nobackup/jkolassa/CF0720x6C_CF0720x6C-Pfafstetter.til';
til_data = dlmread(til_file);
fp_topo_lon = til_data(:,3);
fp_topo_lat = til_data(:,4);
fp_ig = til_data(:,5);
fp_jg = til_data(:,6);

load('fp_topo_geoloc.mat')
%fp_topo_lat_new = NaN*ones(max(fp_ig),max(fp_jg));
%fp_topo_lon_new = NaN*ones(max(fp_ig),max(fp_jg));
%
%for i = 1:max(fp_ig)
%    i
%    for j = 1:max(fp_jg)
%
%        data_index = (fp_ig==i) + (fp_jg==j);
%        data_index = (data_index==2);
%        fp_topo_lat_new(i,j) = nanmean(fp_topo_lat(data_index));
%        fp_topo_lon_new(i,j) = nanmean(fp_topo_lon(data_index));
%     end
%end

% L4 constants file

l4_fpath = '/css/smapl4/public/L4_Products/L4_SM/Vv5030/lmc/';
l4_fname = 'SMAP_L4_SM_lmc_00000000T000000_Vv5030_001.h5';

% ----------------------------------------------------------------------

% read FP data

fp_phis = ncread([fp_fpath,'/',fp_fname],'PHIS');
fp_lon  = ncread([fp_fpath,'/',fp_fname],'lon');
fp_lat  = ncread([fp_fpath,'/',fp_fname],'lat');

%dlon = 1/12;
%fp_topo_lon = (-180+dlon/2):dlon:(180-dlon/2);
fp_topo = ncread(fp_topo_file,'z');
fp_topo = double(fp_topo);

fp_latind = find( latrange(1)<fp_lat & fp_lat<latrange(2) );
fp_lonind = find( lonrange(1)<fp_lon & fp_lon<lonrange(2) );
fp_topo_lonind = (fp_topo_lon_new>lonrange(1)) + (fp_topo_lon_new<lonrange(2) );
fp_topo_lonind = fp_topo_lonind == 2;
fp_topo_latind = (fp_topo_lat_new>latrange(1)) + (fp_topo_lat_new<latrange(2) );
fp_topo_latind = fp_topo_latind==2;
fp_topo_plotind = fp_topo_latind + fp_topo_lonind;
fp_topo_plotind = fp_topo_plotind==2;

% read L4 data

l4_elev = h5read([l4_fpath,'/',l4_fname],'/Land-Model-Constants_Data/cell_elevation');

l4_elev(l4_elev==-9999) = NaN;

l4_lon  = h5read([l4_fpath,'/',l4_fname],'/cell_lon');
l4_lat  = h5read([l4_fpath,'/',l4_fname],'/cell_lat');

l4_lon  =           l4_lon(:,1);
l4_lat  = transpose(l4_lat(1,:));

l4_latind = find( latrange(1)<l4_lat & l4_lat<latrange(2) );
l4_lonind = find( lonrange(1)<l4_lon & l4_lon<lonrange(2) );

% plot

figure
axesm('robinson','MapLatLimit',latrange,'MapLonLimit',lonrange);
pcolorm(fp_lat(fp_latind),fp_lon(fp_lonind),transpose(fp_phis(fp_lonind,fp_latind))/g0);
shading('flat');
geoshow('landareas.shp','EdgeColor','w','FaceColor','none')
load oceanlo;
hm=displaym(oceanmask);
set(hm,'FaceColor',[1 1 1 ])
caxis(elevrange)
colormap(parula(11))
colorbar
title('FP 5.27.1 PHIS/g0 (m)')
drawnow
print('-dpng','-r300','TC_Idai_FP_elev_FP5271_qtrdeg.png');

plot_lat = fp_topo_lat_new(fp_topo_plotind);
plot_lon = fp_topo_lon_new(fp_topo_plotind);
plot_topo = fp_topo(fp_topo_plotind);

%figure
%axesm('robinson','MapLatLimit',latrange,'MapLonLimit',lonrange);
%pcolorm(plot_lat,plot_lon,plot_topo);
%shading('flat');
%geoshow('landareas.shp','EdgeColor','w','FaceColor','none')
%load oceanlo;
%hm=displaym(oceanmask);
%set(hm,'FaceColor',[1 1 1 ])
%caxis(elevrange)
%colormap(parula(11))
%colorbar
%title('TOPO DYN AVE (m)')
%drawnow
%print('-dpng','-r300','TC_Idai_FP_elev_topo_dyn_ave.png');

figure
axesm('robinson','MapLatLimit',[-90 90],'MapLonLimit',[-180 180]);
pcolorm(fp_topo_lat_new,fp_topo_lon_new,fp_topo(:,:));
shading('flat');
geoshow('landareas.shp','EdgeColor','w','FaceColor','none')
load oceanlo;
hm=displaym(oceanmask);
set(hm,'FaceColor',[1 1 1 ])
caxis(elevrange)
colormap(parula(11))
colorbar
title('TOPO DYN AVE (m)')
drawnow
print('-dpng','-r300','TC_Idai_FP_elev_topo_dyn_ave_global.png');

figure
axesm('robinson','MapLatLimit',latrange,'MapLonLimit',lonrange);
pcolorm(double(l4_lat(l4_latind)),double(l4_lon(l4_lonind)),double(transpose(l4_elev(l4_lonind,l4_latind))));
shading('flat');
geoshow('landareas.shp','EdgeColor','w','FaceColor','none')
caxis(elevrange)
colormap(parula(11))
colorbar
title('L4SM cell\_elevation (m)')
drawnow
print('-dpng','-r300','TC_Idai_FP_elev_L4SM_9km.png');
