clc
clear all
close all

% make selections
storm_name = 'Iota';
start_date = '20201103';
stop_date = '20201119';
lat_range = [0 33];      % [southern_border northern_border]
lon_range = [-110 -70];  % [western_border eastern_border]
period_length = 'one_month_prior';

smap_stream = 'Vv5014';
ol_stream = 'OL5030';

clim_file = ['/home/jkolassa/Projects/SMAP_TC/' storm_name '_clim_' period_length '.mat'];
clim_data = load(clim_file);


% convert selections
smap_l4_path = ['/css/smapl4/public/L4_Products/L4_SM/' smap_stream '/gph/'];
smap_ol_path = ['/css/smapl4/public/Model-Only/SoilMoisture/' ol_stream '/gph/'];

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

track_file = ['/home/jkolassa/Projects/SMAP_TC/' storm_name '_track.ascii'];
track_data = dlmread(track_file);
track_datenum = datenum(start_year, track_data(:,1), track_data(:,2), track_data(:,3), 0, 0);

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_evap_diff = VideoWriter(['/home/jkolassa/Projects/SMAP_TC/' storm_name '_evap_anom_prior']);
v_evap_diff.FrameRate = 1;
open(v_evap_diff);


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');
            mat_index = 0;
            for t = 1:length(t_vec)

                mat_index = mat_index + 1;

                smap_l4_file = [smap_l4_path 'Y' y_str '/M' m_str '/D' d_str ...
                           '/SMAP_L4_SM_gph_' y_str m_str d_str 'T' t_vec{t} '00_' smap_stream '_001.h5'];
                smap_ol_file = [smap_ol_path 'Y' y_str '/M' m_str '/D' d_str ...
                           '/SMAP_L4_SM_gph_' y_str m_str d_str 'T' t_vec{t} '00_' ol_stream '_001.h5'];

                if (mat_index==1)
                   lat_data = double(h5read(smap_l4_file,'/cell_lat'));
                   lon_data = double(h5read(smap_l4_file,'/cell_lon'));
                   lat_index = (lat_data(1,:)>=lat_range(1)) + (lat_data(1,:)<=lat_range(2));
                   lat_index = lat_index == 2;
                   lon_index = (lon_data(:,1)>=lon_range(1)) + (lon_data(:,1)<=lon_range(2));
                   lon_index = lon_index == 2;
                   domain_lats = lat_data(lon_index,lat_index);
                   domain_lons = lon_data(lon_index,lat_index);
                   domain_evap_l4 = NaN*ones(size(domain_lats,1),size(domain_lats,2),length(t_vec));
                   domain_evap_ol = NaN*ones(size(domain_lats,1),size(domain_lats,2),length(t_vec));
                   evap_diff_mat = NaN*ones(size(domain_lats,1),size(domain_lats,2));
                   lh_diff_mat = NaN*ones(size(domain_lats,1),size(domain_lats,2));
                end
                sm_l4_evap_data = double(h5read(smap_l4_file,'/Geophysical_Data/land_evapotranspiration_flux'));
                sm_ol_evap_data = double(h5read(smap_ol_file,'/Geophysical_Data/land_evapotranspiration_flux'));
                % select domain and convert all to mm/day
                domain_evap_l4(:,:,mat_index) = sm_l4_evap_data(lon_index,lat_index).*86400-clim_data.clim_et_l4;
                domain_evap_ol(:,:,mat_index) = sm_ol_evap_data(lon_index,lat_index).*86400-clim_data.clim_et_ol;
                end % end of t_vec loop

                % compute daily total
                mean_evap_l4 = nanmean(domain_evap_l4,3);
             %   mean_lh_l4 = nanmean(domain_lh_l4,3);

                % find nearest track point

                smap_datenum = datenum(y,m,d);
                datenum_dist = sqrt((track_datenum-smap_datenum).^2);
                date_index = find(datenum_dist==min(datenum_dist));
% plot domain-wide surface SM increments

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

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

                 figure
                 ax = axesm('miller','MapLatLimit',lat_range,'MapLonLimit',lon_range,...
                 'FLatLimit',lat_range,'FLonLimit',lon_range);
                 framem off; 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(domain_lats,domain_lons,mean_evap_l4)
                  caxis([caxis_min caxis_max]);
                  set(gcf,'Colormap',cmap);
                  plotm(coastlat, coastlon,'k')
                  hold on
                  plotm(track_data(1:date_index,5),track_data(1:date_index,4),'-r','LineWidth',2)
                  hold on
                  %geoshow(states,'facealpha',0);
                  set(gca,'XLim',lonlim);
                  set(gca,'YLim',latlim);
                  tightmap
                  date_time_str = datestr(datetime(smap_datenum,'ConvertFrom','datenum'));
                  title(['land evaporation anomaly [mm] (DA-clim) ' date_time_str])
                  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/SMAP_TC/' storm_name '_evap_anom_prior_' y_str m_str d_str],'png')

                  frame = getframe(gcf);
                  writeVideo(v_evap_diff,frame);
        end

    end
end

close(v_evap_diff);

