clc
clear all
close all

% make selections
storm_name = 'Idai';
start_date = '20190223';
stop_date = '20190323';
lat_range = [-24 -10];      % [southern_border northern_border]
lon_range = [25 52];  % [western_border eastern_border]
save_path = '/home/jkolassa/Projects/SMAP_TC/';
fp_stream_name = 'f521_fp';

% convert selections
smap_l4_path = ['/discover/nobackup/jkolassa/SMAP_TC_data/'];

start_year = str2double(start_date(1:4));
start_month = str2double(start_date(5:6));
start_day = str2double(start_date(7:8));
stop_year = str2double(stop_date(1:4));
stop_month = str2double(stop_date(5:6));
stop_day = str2double(stop_date(7:8));

% read track data
if strcmp(storm_name,'Eta') | strcmp(storm_name,'Iota')
track_file_1 = ['/home/jkolassa/Projects/SMAP_TC/' storm_name '_track.ascii'];
track_data_1 = dlmread(track_file_1);
track_datenum_1 = datenum(start_year, track_data_1(:,1), track_data_1(:,2), track_data_1(:,3), 0, 0); 

track_file_2 = ['/home/jkolassa/Projects/SMAP_TC/Iota_track.ascii'];
track_data_2 = dlmread(track_file_2);
track_datenum_2 = datenum(start_year, track_data_2(:,1), track_data_2(:,2), track_data_2(:,3), 0, 0);
else
track_file_1 = ['/home/jkolassa/Projects/SMAP_TC/' storm_name '_track.ascii'];
track_data_1 = dlmread(track_file_1);
track_datenum_1 = datenum(start_year, track_data_1(:,1), track_data_1(:,2), track_data_1(:,3), 0, 0);
end

t_vec = {'0130' '0430' '0730' '1030' '1330' '1630' '1930' '2230'};
month_length = [31 28 31 30 31 30 31 31 30 31 30 31];

mat_index = 0;

v_sm = VideoWriter([save_path storm_name '_L4_SFMC_anom_daily']);
v_sm.FrameRate = 4;
open(v_sm);

v_rz = VideoWriter([save_path storm_name '_L4_RZMC_anom_daily']);
v_rz.FrameRate = 4;
open(v_rz);

v_et = VideoWriter([save_path storm_name '_L4_ET_anom_daily']);
v_et.FrameRate = 4;
open(v_et);

for y = start_year:stop_year

    y_str = num2str(y);

    m_first = 1;
    m_last = 12;
    if (y==start_year)
        m_first = start_month;
    end

    if (y==stop_year)
        m_last = stop_month;
    end

    for m = m_first:m_last

        m_str = num2str(m,'%02i');

        day_first = 1;
        day_last = month_length(m);
        if (y==start_year)&(m==start_month)
           day_first = start_day;
        end
        if (y==stop_year)&(m==stop_month)
                    day_last = stop_day;
        end

        for d = day_first:day_last

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


            for t = 1:length(t_vec)

                mat_index = mat_index + 1;
                smap_datenum = datenum(y,m,d,str2double(t_vec{t}(1:2)),str2double(t_vec{t}(3:4)),0);   

                if (strcmp(storm_name,'Eta')) & (smap_datenum>datenum(2020,11,12,0,0,0))
                    storm_name = 'Iota';
                end
                % read SMAP L4 data
                smap_l4_file = [smap_l4_path storm_name '_L4_SM_anom_' y_str m_str d_str '_' t_vec{t} 'z.nc'];
                lat_data = ncread(smap_l4_file,'lat');
                max_lat = max(lat_data(:));
                min_lat = min(lat_data(:));
                lon_data = ncread(smap_l4_file,'lon');
                max_lon = max(lon_data(:));
                min_lon = min(lon_data(:));
                sfmc_anom_data = ncread(smap_l4_file,'sfmc_anomaly');
                rzmc_anom_data = ncread(smap_l4_file,'rzmc_anomaly');
                et_anom_data = ncread(smap_l4_file,'et_anomaly');

                lat_range = [min_lat max_lat];
                lon_range = [min_lon max_lon];

                % find nearest track point
 
                if strcmp(storm_name,'Eta') | strcmp(storm_name,'Iota')
                   datenum_dist_1 = sqrt((track_datenum_1-smap_datenum).^2);
                date_index_1 = find(datenum_dist_1==min(datenum_dist_1));
                   datenum_dist_2 = sqrt((track_datenum_2-smap_datenum).^2);
                date_index_2 = find(datenum_dist_2==min(datenum_dist_2));
                else
                datenum_dist_1 = sqrt((track_datenum_1-smap_datenum).^2);
                date_index_1 = find(datenum_dist_1==min(datenum_dist_1));
                end

                smap_int_min = smap_datenum - (1/24*1.5);
                smap_int_max = smap_datenum + (1/24*1.5);
                 % read atmospheric data from FP

                for hh = 0:23
                    hh_str = num2str(hh,'%02i');
                    inst_datenum = datenum(y,m,d,hh,0,0);
                    if (inst_datenum >= smap_int_min) & (inst_datenum<smap_int_max)
                fp_file = ['/archive/u/dao_ops/GEOS-5.25/GEOSadas-5_25/' fp_stream_name '/ana/Y' y_str '/M' m_str '/' fp_stream_name '.ana.inst3d_met_p.' y_str m_str d_str '_' hh_str 'z.nc4']
                fp_file = ['/archive/u/dao_ops/GEOS-5.21/GEOSadas-5_21/' fp_stream_name '/ana/Y' y_str '/M' m_str '/' fp_stream_name '.ana.inst3d_met_p.' y_str m_str d_str '_' hh_str 'z.nc4']
                fp_lon = ncread(fp_file,'lon');
                fp_lat = ncread(fp_file,'lat');
                fp_lev = ncread(fp_file,'lev');
                fp_u_data = ncread(fp_file,'U');
                fp_v_data = ncread(fp_file,'V');
                lat_index = (fp_lat< min_lat) + (fp_lat > max_lat);
                lat_index = lat_index == 0;
                lon_index = (fp_lon < min_lon) + (fp_lon > max_lon);
                lon_index = lon_index == 0;
                lev_index = fp_lev == 900;
                fp_domain_lats = fp_lat(lat_index);
                fp_domain_lons = fp_lon(lon_index);
                u_900 = fp_u_data(lon_index,lat_index,lev_index);
                v_900 = fp_v_data(lon_index,lat_index,lev_index);


                [X,Y] = meshgrid(fp_domain_lons,fp_domain_lats);
                plot_u = u_900';
                plot_v = v_900';                 

                load coastlines;
                [latcells, loncells] = polysplit(coastlat, coastlon);
                 latlim = lat_range;
                 lonlim = lon_range;

                 cmap_mat = cmocean('balance',20);
                 cmap_mat = flipud(cmap_mat);
                 cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
                 cmap = cmap_mat;
             %    s1 = quantile(sfmc_anom_data(:),0.99);
              %   s2 = quantile(sfmc_anom_data(:),0.01);
              %   caxis_min = -1*max(abs(s1),abs(s2));
              %   caxis_max = max(abs(s1),abs(s2));
                 caxis_min = -0.15;
                 caxis_max = 0.15;

                 max_thin = min(size(X,1),size(X,2));
                 thin_index = 1:1:max_thin;


                 figure
                 ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
                 'FLatLimit',lat_range,'FLonLimit',lon_range);
                 framem
                 gridm on;
                 setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...
                 'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...
                 'MLabelParallel','south','FontSize',12)
                  geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
                  pcolorm(lat_data,lon_data,sfmc_anom_data)
                  caxis([caxis_min caxis_max]);
                  set(gcf,'Colormap',cmap);
                  plotm(coastlat, coastlon,'k')
                  hold on
                  plotm(track_data_1(1:date_index_1,5),track_data_1(1:date_index_1,4),'-r','LineWidth',2)
                  hold on
                  if strcmp(storm_name,'Iota')
                     plotm(track_data_2(1:date_index_2,5),track_data_2(1:date_index_2,4),'-b','LineWidth',2)                     
                     hold on 
                  end
                 % h=streamslice(ax,X,Y,u_900,v_900);
                  quiverm(Y(thin_index,thin_index),X(thin_index,thin_index),plot_v(thin_index,thin_index),...
                          plot_u(thin_index,thin_index),'m');
                  %trimcart(h)
                  hold off
                  %geoshow(states,'facealpha',0);
                  set(gca,'XLim',lonlim);
                  set(gca,'YLim',latlim);
                  tightmap
                  title(['L4 SFMC anomaly [m^{3} m^{-3}] ' y_str m_str d_str '-'  hh_str 'z'])
                  h = colorbar;
                  set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
                  set(gcf,'paperunits','centimeters')
                  set(gcf,'paperposition',[1 1 17 15])
 %                 saveas(gcf,[save_path storm_name '_L4_SFMC_anom_flow_daily_' y_str m_str d_str t_vec{t}],'png')

                  frame = getframe(gcf);
                  writeVideo(v_sm,frame);

                  cmap_mat = cmocean('balance',20);
                 cmap_mat = flipud(cmap_mat);
                 cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
                 cmap = cmap_mat;
             %    s1 = quantile(sfmc_anom_data(:),0.99);
              %   s2 = quantile(sfmc_anom_data(:),0.01);
              %   caxis_min = -1*max(abs(s1),abs(s2));
              %   caxis_max = max(abs(s1),abs(s2));
                 caxis_min = -0.1;
                 caxis_max = 0.1;

                  figure
                 ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
                 'FLatLimit',lat_range,'FLonLimit',lon_range);
                 framem
                 gridm on;
                 setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...
                 'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...
                 'MLabelParallel','south','FontSize',12)
                  geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
                  pcolorm(lat_data,lon_data,rzmc_anom_data)
                  caxis([caxis_min caxis_max]);
                  set(gcf,'Colormap',cmap);
                  plotm(coastlat, coastlon,'k')
                  hold on
                  plotm(track_data_1(1:date_index_1,5),track_data_1(1:date_index_1,4),'-r','LineWidth',2)                     
                  hold on                                                                                                     
                  if strcmp(storm_name,'Iota')                                                                                
                     plotm(track_data_2(1:date_index_2,5),track_data_2(1:date_index_2,4),'-b','LineWidth',2)                     
                     hold on                                                                                                  
                  end 
                 % h=streamslice(ax,X,Y,u_900,v_900);
                  quiverm(Y(thin_index,thin_index),X(thin_index,thin_index),plot_v(thin_index,thin_index),...
                          plot_u(thin_index,thin_index),'m');
                  %trimcart(h)
                  hold off
                  %geoshow(states,'facealpha',0);
                  set(gca,'XLim',lonlim);
                  set(gca,'YLim',latlim);
                  tightmap
                  title(['L4 RZMC anomaly [m^{3} m^{-3}] ' y_str m_str d_str '-'  hh_str 'z'])
                  h = colorbar;
                  set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
                  set(gcf,'paperunits','centimeters')
                  set(gcf,'paperposition',[1 1 17 15])
%                  saveas(gcf,[save_path storm_name '_L4_RZMC_anom_flow_daily_' y_str m_str d_str t_vec{t}],'png')

                  frame = getframe(gcf);
                  writeVideo(v_rz,frame);

                  cmap_mat = cmocean('balance',20);
                 cmap_mat = flipud(cmap_mat);
                 cmap_mat(length(cmap_mat)/2:length(cmap_mat)/2+1,:) = ones(2,3);
                 cmap = cmap_mat;
                 caxis_min = -4;
                 caxis_max = 4;

                   figure
                 ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
                 'FLatLimit',lat_range,'FLonLimit',lon_range);
                 framem
                 gridm on;
                 setm(ax,'MLineLocation',5,'PLineLocation',5,'MeridianLabel','on',...
                 'ParallelLabel','on','MLabelLocation',5,'PLabelLocation',5,...
                 'MLabelParallel','south','FontSize',12)
                  geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
                  pcolorm(lat_data,lon_data,et_anom_data)
                  caxis([caxis_min caxis_max]);
                  set(gcf,'Colormap',cmap);
                  plotm(coastlat, coastlon,'k')
                  hold on
                  plotm(track_data_1(1:date_index_1,5),track_data_1(1:date_index_1,4),'-r','LineWidth',2)                     
                  hold on                                                                                                     
                  if strcmp(storm_name,'Iota')                                                                                
                     plotm(track_data_2(1:date_index_2,5),track_data_2(1:date_index_2,4),'-b','LineWidth',2)                     
                     hold on                                                                                                  
                  end 
                 % h=streamslice(ax,X,Y,u_900,v_900);
                  quiverm(Y(thin_index,thin_index),X(thin_index,thin_index),plot_v(thin_index,thin_index),...
                          plot_u(thin_index,thin_index),'m');
                  %trimcart(h)
                  hold off
                  %geoshow(states,'facealpha',0);
                  set(gca,'XLim',lonlim);
                  set(gca,'YLim',latlim);
                  tightmap
                  title(['L4 ET anomaly [mm/day] ' y_str m_str d_str '-' hh_str 'z'])
                  h = colorbar;
                  set(h,'Location','eastoutside','Orientation','vertical','FontSize',14);
                  set(gcf,'paperunits','centimeters')
                  set(gcf,'paperposition',[1 1 17 15])
%                  saveas(gcf,[save_path storm_name '_L4_ET_anom_flow_daily_' y_str m_str d_str t_vec{t}],'png')

                  frame = getframe(gcf);
                  writeVideo(v_et,frame);

                end % if
               end % hh
            end %% t_vec
         end % d
    end % m
end % y

close(v_sm);
close(v_rz);
close(v_et);
