clc
clear all 
close all

%-----------------------------------------
% EXPERIMENT DEFINITION
%-----------------------------------------

storm_name = 'Idai';
save_path = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/' storm_name '/'];

control_dir = '/discover/nobackup/projects/smaptc/jkolassa/archive/';
control_name = 'LADAS_control_4D_new';
exp_dir = '/discover/nobackup/projects/smaptc/jkolassa/';
exp_name = 'Idai_SMAP_4D';

lat_range = [-24 -10];
lon_range = [24 52];

start_year = 2019;
start_month = 2;
start_day = 4;

end_year = 2019;
end_month = 3;
end_day = 22;

land_mask = 'on'; % when making summary flux plots, use entire domain or mask to only use land pixels expected to influence storm

%---------------------------------------------
% DATA READ PREPARATION
%---------------------------------------------

leapYears = 1980:4:2022;
monthLength_nly = [31 28 31 30 31 30 31 31 30 31 30 31];
monthLength_ly = [31 29 31 30 31 30 31 31 30 31 30 31];
mat_index = 0;
mask_index  = 0;

%---------------------------------------------
% DATA READ
%---------------------------------------------

% read track data

track_file = ['/discover/nobackup/jkolassa/SMAP_TC_project/track_files/' 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);

% loop over years
for y = start_year:end_year

    y_str = num2str(y,'%04i');
 
    % define month length based on leap year status
    if sum(y==leapYears)>0
        monthLength = monthLength_ly;
    else
        monthLength = monthLength_nly;
    end

    % loop over months
    for m = start_month:end_month

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

        % define start and end day for current month
        day_first = 1;
        day_last = monthLength(m);

        if (m == start_month)
           day_first = start_day;
        end

        if (m == end_month)
           day_last = end_day;
        end

        % loop over days
        for d = day_first:day_last

            mat_index = mat_index + 1;

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

            % loop over time steps
            for t = 0:23

                mask_index = mask_index + 1;

                t_str = num2str(t,'%02i');

                % define diagnostic file names
                control_diag_file = [control_dir control_name '/diag/Y' y_str '/M' m_str '/' control_name ...
                                    '.tavg1_2d_lnd_Nx.' y_str m_str d_str '_' t_str '30z.nc4'];

                exp_diag_file = [exp_dir exp_name '/diag_old/Y' y_str '/M' m_str '/' exp_name ... 
                                '.tavg1_2d_lnd_Nx.' y_str m_str d_str '_' t_str '30z.nc4'];

                control_ana_file = [control_dir control_name '/ana/Y' y_str '/M' m_str '/' control_name ...
                                    '.ana.inst3d_met_p.' y_str m_str d_str '_' t_str '00z.nc4'];

                exp_ana_file = [exp_dir exp_name '/ana/Y' y_str '/M' m_str '/' exp_name ...
                                '.ana.inst3d_met_p.' y_str m_str d_str '_' t_str '00z.nc4'];

                control_ocn_file = [control_dir control_name '/diag/Y' y_str '/M' m_str '/' control_name ...
                                    '.tavg1_2d_ocn_Nx.' y_str m_str d_str '_' t_str '30z.nc4'];

                exp_ocn_file = [exp_dir exp_name '/diag_old/Y' y_str '/M' m_str '/' exp_name ...
                                '.tavg1_2d_ocn_Nx.' y_str m_str d_str '_' t_str '30z.nc4'];

                % on first time step read coordinate information and initialize matrices

                if (mat_index == 1) && (t == 0)

                   lat_vec = double(ncread(control_diag_file,'lat'));
                   lon_vec = double(ncread(control_diag_file,'lon'));

                   datenum_vec = NaN*ones(sum(monthLength(start_month:end_month)),1);

                   control_sfmc_tavg24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
                   control_rzmc_tavg24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
                   control_totevptrns_ttot24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
                   control_precip_ttot24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
                   control_runoff_ttot24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));

                   exp_sfmc_tavg24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
                   exp_rzmc_tavg24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
                   exp_totevptrns_ttot24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
                   exp_precip_ttot24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
                   exp_runoff_ttot24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));

                   control_land_mask_daily = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
                   exp_land_mask_daily = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));

                   control_sfmc_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   control_rzmc_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   control_evpintr_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   control_evpsoil_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   control_evptrns_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   control_totevptrns_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   control_precip_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   control_runoff_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);

                   exp_sfmc_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   exp_rzmc_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   exp_evpintr_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   exp_evpsoil_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   exp_evptrns_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   exp_totevptrns_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   exp_precip_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
                   exp_runoff_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);

                   exp_land_mask = NaN*ones(length(lon_vec),length(lat_vec),24);
                   control_land_mask = NaN*ones(length(lon_vec),length(lat_vec),24);

                end % if

                % read hourly data

                control_sfmc_hour_mat(:,:,t+1) = double(ncread(control_diag_file,'SFMC'));
                control_rzmc_hour_mat(:,:,t+1) = double(ncread(control_diag_file,'RZMC'));
                control_evpintr_hour_mat(:,:,t+1) = double(ncread(control_diag_file,'EVPINTR'));
                control_evpsoil_hour_mat(:,:,t+1) = double(ncread(control_diag_file,'EVPSOIL'));
                control_evptrns_hour_mat(:,:,t+1) = double(ncread(control_diag_file,'EVPTRNS'));
                control_precip_hour_mat(:,:,t+1) = double(ncread(control_diag_file,'PRECTOT'));
                control_runoff_hour_mat(:,:,t+1) = double(ncread(control_diag_file,'RUNOFF'));

                exp_sfmc_hour_mat(:,:,t+1) = double(ncread(exp_diag_file,'SFMC'));
                exp_rzmc_hour_mat(:,:,t+1) = double(ncread(exp_diag_file,'RZMC'));
                exp_evpintr_hour_mat(:,:,t+1) = double(ncread(exp_diag_file,'EVPINTR'));
                exp_evpsoil_hour_mat(:,:,t+1) = double(ncread(exp_diag_file,'EVPSOIL'));
                exp_evptrns_hour_mat(:,:,t+1) = double(ncread(exp_diag_file,'EVPTRNS'));
                exp_precip_hour_mat(:,:,t+1) = double(ncread(exp_diag_file,'PRECTOT'));
                exp_runoff_hour_mat(:,:,t+1) = double(ncread(exp_diag_file,'RUNOFF'));

                % read in wind data for filtering

                control_hPa = double(ncread(control_ana_file,'lev'));
                control_level_index = find(control_hPa==850);
                tmp_mat = double(ncread(control_ana_file,'U'));
                control_u_mat = tmp_mat(:,:,control_level_index);
                tmp_mat = double(ncread(control_ana_file,'V'));
                control_v_mat = tmp_mat(:,:,control_level_index);

                control_ana_lat = double(ncread(control_ana_file,'lat'));
                control_ana_lon = double(ncread(control_ana_file,'lon')); 

                exp_hPa = double(ncread(exp_ana_file,'lev'));
                exp_level_index = find(exp_hPa==850);
                tmp_mat = double(ncread(exp_ana_file,'U'));
                exp_u_mat = tmp_mat(:,:,exp_level_index);
                tmp_mat = double(ncread(exp_ana_file,'V'));
                exp_v_mat = tmp_mat(:,:,exp_level_index);

                exp_ana_lat = double(ncread(exp_ana_file,'lat'));
                exp_ana_lon = double(ncread(exp_ana_file,'lon'));

                [exp_ana_lat_mat,exp_ana_lon_mat] = meshgrid(exp_ana_lat,exp_ana_lon);
                control_u_interp_mat = interp2(control_ana_lat,control_ana_lon,control_u_mat,exp_ana_lat_mat,exp_ana_lon_mat);
                control_v_interp_mat = interp2(control_ana_lat,control_ana_lon,control_v_mat,exp_ana_lat_mat,exp_ana_lon_mat);

                % compute total wind speed from components

                control_ws_mat = sqrt((control_u_interp_mat).^2 + (control_v_interp_mat).^2);
                exp_ws_mat = sqrt((exp_u_mat).^2 + (exp_v_mat).^2);

                % read land fraction

                control_frland = double(ncread(control_ocn_file,'FRLAND'));
                exp_frland = double(ncread(exp_ocn_file,'FRLAND'));

                % compute land influence mask

               exp_land_mask(:,:,t+1) = (exp_ws_mat>=17) & (exp_frland>=0.5);
               control_land_mask(:,:,t+1) = (control_ws_mat>=17) & (control_frland>=0.5);

                % compute total evaporation

                control_totevptrns_hour_mat(:,:,t+1) = control_evptrns_hour_mat(:,:,t+1) + ...
                                                       control_evpsoil_hour_mat(:,:,t+1) + ...
                                                       control_evpintr_hour_mat(:,:,t+1);

                exp_totevptrns_hour_mat(:,:,t+1) = exp_evptrns_hour_mat(:,:,t+1) + ...
                                                   exp_evpsoil_hour_mat(:,:,t+1) + ...
                                                   exp_evpintr_hour_mat(:,:,t+1);


                 % filter out points that are not expected to influence the storm

                 if (strcmp(land_mask,'on'))
           
                    tmp_mat = control_totevptrns_hour_mat(:,:,t+1);
                    tmp_mat(control_land_mask(:,:,t+1)==0) = NaN;
                    control_totevptrns_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = exp_totevptrns_hour_mat(:,:,t+1);
                    tmp_mat(exp_land_mask(:,:,t+1)==0) = NaN;
                    exp_totevptrns_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = control_precip_hour_mat(:,:,t+1);
                    tmp_mat(control_land_mask(:,:,t+1)==0) = NaN; 
                    control_precip_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = exp_precip_hour_mat(:,:,t+1);
                    tmp_mat(exp_land_mask(:,:,t+1)==0) = NaN; 
                    exp_precip_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = control_runoff_hour_mat(:,:,t+1);
                    tmp_mat(control_land_mask(:,:,t+1)==0) = NaN; 
                    control_runoff_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = exp_runoff_hour_mat(:,:,t+1);
                    tmp_mat(exp_land_mask(:,:,t+1)==0) = NaN; 
                    exp_runoff_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = control_sfmc_hour_mat(:,:,t+1);
                    tmp_mat(control_land_mask(:,:,t+1)==0) = NaN; 
                    control_sfmc_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = exp_sfmc_hour_mat(:,:,t+1);
                    tmp_mat(exp_land_mask(:,:,t+1)==0) = NaN; 
                    exp_sfmc_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = control_rzmc_hour_mat(:,:,t+1);
                    tmp_mat(control_land_mask(:,:,t+1)==0) = NaN;
                    control_rzmc_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = exp_rzmc_hour_mat(:,:,t+1);
                    tmp_mat(exp_land_mask(:,:,t+1)==0) = NaN;
                    exp_rzmc_hour_mat(:,:,t+1) = tmp_mat;

                 end

%                if ((m==3) & (d==12) & (t==0))
%
%                    test_control_sfmc = control_sfmc_hour_mat(:,:,t+1);
%                    test_exp_sfmc = exp_sfmc_hour_mat(:,:,t+1);
%   
%                end
%
            end  %t

            % compute mm/hr temporary matrices for daily total calculations

            tmp_cntrl_totevptrns_mat = control_totevptrns_hour_mat*0.0036*0.408;
            tmp_exp_totevptrns_mat = exp_totevptrns_hour_mat*0.0036*0.408;
            tmp_cntrl_precip_mat = control_precip_hour_mat*3600;
            tmp_exp_precip_mat = exp_precip_hour_mat*3600;
            tmp_cntrl_runoff_mat = control_runoff_hour_mat*3600;
            tmp_exp_runoff_mat = exp_runoff_hour_mat*3600;

            % convert total evaporation from W m^-2 to mm/day

            control_totevptrns_hour_mat = control_totevptrns_hour_mat*0.0864*0.408;
            exp_totevptrns_hour_mat = exp_totevptrns_hour_mat*0.0864*0.408;

            % log date

            datenum_vec(mat_index,1) = datenum(y,m,d);

            % compute daily averages

            control_sfmc_tavg24_mat(:,:,mat_index) = nanmean(control_sfmc_hour_mat,3);
            control_rzmc_tavg24_mat(:,:,mat_index) = nanmean(control_rzmc_hour_mat,3);

            control_totevptrns_ttot24_mat(:,:,mat_index) = sum(tmp_cntrl_totevptrns_mat,3,'omitnan');
            nan_index = (sum(isnan(tmp_cntrl_totevptrns_mat))==size(tmp_cntrl_totevptrns_mat,3));
            tmp = control_totevptrns_ttot24_mat(:,:,mat_index);
            tmp(nan_index) = NaN;
            control_totevptrns_ttot24_mat(:,:,mat_index) = tmp;

            control_precip_ttot24_mat(:,:,mat_index) = sum(tmp_cntrl_precip_mat,3,'omitnan');
            nan_index = (sum(isnan(tmp_cntrl_precip_mat))==size(tmp_cntrl_precip_mat,3));
            tmp = control_precip_ttot24_mat(:,:,mat_index);
            tmp(nan_index) = NaN;
            control_precip_ttot24_mat(:,:,mat_index) = tmp;

            control_runoff_ttot24_mat(:,:,mat_index) = sum(tmp_cntrl_runoff_mat,3,'omitnan');
            nan_index = (sum(isnan(tmp_cntrl_runoff_mat))==size(tmp_cntrl_runoff_mat,3));
            tmp = control_runoff_ttot24_mat(:,:,mat_index);
            tmp(nan_index) = NaN;
            control_runoff_ttot24_mat(:,:,mat_index) = tmp;

            exp_sfmc_tavg24_mat(:,:,mat_index) = nanmean(exp_sfmc_hour_mat,3);
            exp_rzmc_tavg24_mat(:,:,mat_index) = nanmean(exp_rzmc_hour_mat,3);

            exp_totevptrns_ttot24_mat(:,:,mat_index) = sum(tmp_exp_totevptrns_mat,3,'omitnan');
            nan_index = (sum(isnan(tmp_exp_totevptrns_mat))==size(tmp_exp_totevptrns_mat,3));
            tmp = exp_totevptrns_ttot24_mat(:,:,mat_index);
            tmp(nan_index) = NaN;
            exp_totevptrns_ttot24_mat(:,:,mat_index) = tmp;

            exp_precip_ttot24_mat(:,:,mat_index) = sum(tmp_exp_precip_mat,3,'omitnan');
            nan_index = (sum(isnan(tmp_exp_precip_mat))==size(tmp_exp_precip_mat,3));
            tmp = exp_precip_ttot24_mat(:,:,mat_index);
            tmp(nan_index) = NaN;
            exp_precip_ttot24_mat(:,:,mat_index) = tmp;

            exp_runoff_ttot24_mat(:,:,mat_index) = sum(tmp_exp_runoff_mat,3,'omitnan');
            nan_index = (sum(isnan(tmp_exp_runoff_mat))==size(tmp_exp_runoff_mat,3));
            tmp = exp_runoff_ttot24_mat(:,:,mat_index);
            tmp(nan_index) = NaN;
            exp_runoff_ttot24_mat(:,:,mat_index) = tmp;

            control_land_mask_daily(:,:,mat_index) = sum(control_land_mask,3,'omitnan');
            exp_land_mask_daily(:,:,mat_index) = sum(exp_land_mask,3,'omitnan');

            % reset hourly matrices

            control_sfmc_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            control_rzmc_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            control_evpintr_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            control_evpsoil_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            control_evptrns_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            control_totevptrns_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);   
            control_precip_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);         
            control_runoff_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);

            exp_sfmc_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            exp_rzmc_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            exp_evpintr_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            exp_evpsoil_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            exp_evptrns_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            exp_totevptrns_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            exp_precip_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);
            exp_runoff_hour_mat = NaN*ones(length(lon_vec),length(lat_vec),24);

            control_land_mask = NaN*ones(length(lon_vec),length(lat_vec),24);
            exp_land_mask = NaN*ones(length(lon_vec),length(lat_vec),24);

        end %d
    end %m
end %y

% remove excess allocated slots

datenum_vec((mat_index+1):end) = [];
control_sfmc_tavg24_mat(:,:,(mat_index+1):end) = [];
exp_sfmc_tavg24_mat(:,:,(mat_index+1):end) = [];
control_rzmc_tavg24_mat(:,:,(mat_index+1):end) = [];
exp_rzmc_tavg24_mat(:,:,(mat_index+1):end) = [];
control_totevptrns_ttot24_mat(:,:,(mat_index+1):end) = [];
exp_totevptrns_ttot24_mat(:,:,(mat_index+1):end) = [];
control_precip_ttot24_mat(:,:,(mat_index+1):end) = [];
exp_precip_ttot24_mat(:,:,(mat_index+1):end) = [];
control_runoff_ttot24_mat(:,:,(mat_index+1):end) = [];
exp_runoff_ttot24_mat(:,:,(mat_index+1):end) = [];
control_land_mask_daily(:,:,(mat_index+1):end) = [];
exp_land_mask_daily(:,:,(mat_index+1):end) = [];

% compute differences between experiment and control

diff_sfmc_tavg24_mat = exp_sfmc_tavg24_mat - control_sfmc_tavg24_mat;
diff_rzmc_tavg24_mat = exp_rzmc_tavg24_mat - control_rzmc_tavg24_mat;
diff_totevptrns_ttot24_mat = exp_totevptrns_ttot24_mat - control_totevptrns_ttot24_mat;
diff_precip_ttot24_mat = exp_precip_ttot24_mat - control_precip_ttot24_mat;
diff_runoff_ttot24_mat = exp_runoff_ttot24_mat - control_runoff_ttot24_mat;


% create animation and plot of differences

% define plot domain
[lon_mat,lat_mat] = meshgrid(lon_vec,lat_vec);
load coastlines;
[latcells, loncells] = polysplit(coastlat, coastlon);
latlim = lat_range;
lonlim = lon_range;

% define colorbars and limits

lat_index = (lat_vec<lat_range(1)) + (lat_vec>lat_range(2));
lat_index = lat_index == 0;
lon_index = (lon_vec<lon_range(1)) + (lon_vec>lon_range(2));
lon_index = lon_index == 0;

cmap_land_mask = [ 1 1 1; 255/255 51/255 51/255];
caxis_min_land_mask = 0;
caxis_max_land_mask = 1;


v_cm = VideoWriter([save_path storm_name '_control_land_mask']);
v_cm.FrameRate = 4;
open(v_cm);

v_em = VideoWriter([save_path storm_name '_exp_land_mask']);
v_em.FrameRate = 4;
open(v_em);

for i = 1:length(datenum_vec)

    datenum_dist = sqrt((track_datenum-datenum_vec(i)).^2);
    date_index = find(datenum_dist==min(datenum_dist));

    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_vec,lon_vec,double(control_land_mask_daily(:,:,i))')
    caxis([caxis_min_land_mask caxis_max_land_mask]);
    set(gcf,'Colormap',cmap_land_mask);
    plotm(coastlat, coastlon,'k')
    hold on
    plotm(track_data(1:date_index,5),track_data(1:date_index,4),'-m','LineWidth',2)
    hold off
    set(gca,'XLim',lonlim);
    set(gca,'YLim',latlim);
    tightmap
    title(['control land mask ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    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 '_diag_SFMC_diff_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')

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

    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_vec,lon_vec,double(exp_land_mask_daily(:,:,i))')
    caxis([caxis_min_land_mask caxis_max_land_mask]);
    set(gcf,'Colormap',cmap_land_mask);
    plotm(coastlat, coastlon,'k')
    hold on
    plotm(track_data(1:date_index,5),track_data(1:date_index,4),'-m','LineWidth',2)
    hold off
    set(gca,'XLim',lonlim);
    set(gca,'YLim',latlim);
    tightmap
    title(['exp land mask ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    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 '_diag_SFMC_diff_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')

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

  % compute average values over storm domain

    tmp_mat = control_sfmc_tavg24_mat(lon_index,lat_index,i);
    control_mean_sfmc(i) = nanmean(tmp_mat(:));
    tmp_mat = control_rzmc_tavg24_mat(lon_index,lat_index,i);
    control_mean_rzmc(i) = nanmean(tmp_mat(:));
    tmp_mat = control_totevptrns_ttot24_mat(lon_index,lat_index,i);
    control_mean_totevptrns(i) = nanmean(tmp_mat(:));
    tmp_mat = control_precip_ttot24_mat(lon_index,lat_index,i);
    control_mean_precip(i) = nanmean(tmp_mat(:));
    tmp_mat = control_runoff_ttot24_mat(lon_index,lat_index,i);
    control_mean_runoff(i) = nanmean(tmp_mat(:));

    tmp_mat = exp_sfmc_tavg24_mat(lon_index,lat_index,i);
    exp_mean_sfmc(i) = nanmean(tmp_mat(:));
    tmp_mat = exp_rzmc_tavg24_mat(lon_index,lat_index,i);
    exp_mean_rzmc(i) = nanmean(tmp_mat(:));
    tmp_mat = exp_totevptrns_ttot24_mat(lon_index,lat_index,i);
    exp_mean_totevptrns(i) = nanmean(tmp_mat(:));
    tmp_mat = exp_precip_ttot24_mat(lon_index,lat_index,i);
    exp_mean_precip(i) = nanmean(tmp_mat(:));
    tmp_mat = exp_runoff_ttot24_mat(lon_index,lat_index,i);
    exp_mean_runoff(i) = nanmean(tmp_mat(:));

    tmp_mat = diff_sfmc_tavg24_mat(lon_index,lat_index,i);
    diff_mean_sfmc(i) = nanmean(tmp_mat(:));
    tmp_mat = diff_rzmc_tavg24_mat(lon_index,lat_index,i);
    diff_mean_rzmc(i) = nanmean(tmp_mat(:));
    tmp_mat = diff_totevptrns_ttot24_mat(lon_index,lat_index,i);
    diff_mean_totevptrns(i) = nanmean(tmp_mat(:));
    tmp_mat = diff_precip_ttot24_mat(lon_index,lat_index,i);
    diff_mean_precip(i) = nanmean(tmp_mat(:));
    tmp_mat = diff_runoff_ttot24_mat(lon_index,lat_index,i);
    diff_mean_runoff(i) = nanmean(tmp_mat(:));

end

close(v_cm)
close(v_em)

%close(v_sm);
%close(v_rz);
%close(v_evpintr);
%close(v_evpsoil);
%close(v_evptrns);
%close(v_totevptrns);

figure
yyaxis left
b = bar(datenum_vec,diff_mean_precip);
yyaxis right
p = plot(datenum_vec,diff_mean_sfmc);
title('Area average differences (EXP-CNTRL)')
datetick('x','dd-mmm')
xlabel('Day')
yyaxis left
ylabel('Daily avg precipitation [mm]')
yyaxis right
ylabel('Soil moisture [m^3 m^{-3}]')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,[save_path storm_name '_areaAvgDiff_landMask'],'png')

% convert soil moisture from m^3 m^-3 to mm 
% using 5cm depth for surface and 100 cm depth for root zone

control_mean_sfmc_mm = control_mean_sfmc * 50;
exp_mean_sfmc_mm = exp_mean_sfmc * 50;
control_mean_rzmc_mm = control_mean_rzmc * 1000; 
exp_mean_rzmc_mm = exp_mean_rzmc * 1000;

mean_precip_mat = [control_mean_precip; exp_mean_precip];
gray_mat = [0.4 0.4 0.4; 0.7 0.7 0.7];
blue_vec = [30 136 229]./250;
dark_blue_vec = [1 87 155]./250;
green_vec = [0 137 123]./250;
red_vec = [183 28 28]./250;


figure
yyaxis left
b = bar(datenum_vec,mean_precip_mat);
set(b(1),'FaceColor',gray_mat(1,:));
set(b(2),'FaceColor',gray_mat(2,:));
hold on
p5 = plot(datenum_vec,control_mean_totevptrns,'Color',green_vec,'LineStyle','-','LineWidth',2);
hold on
p6 = plot(datenum_vec,exp_mean_totevptrns,'Color',green_vec,'LineStyle','--','LineWidth',2);
hold on
p7 = plot(datenum_vec,control_mean_runoff,'Color',red_vec,'LineStyle','-','LineWidth',2,'Marker','none');
hold on
p8 = plot(datenum_vec,exp_mean_runoff,'Color',red_vec,'LineStyle','--','LineWidth',2,'Marker','none');
hold off
yyaxis right
p1 = plot(datenum_vec,control_mean_sfmc,'Color',blue_vec,'LineStyle','-','LineWidth',2);
hold on 
p2 = plot(datenum_vec,exp_mean_sfmc,'Color',blue_vec,'LineStyle','--','LineWidth',2);
hold on
p3 = plot(datenum_vec,control_mean_rzmc,'Color',dark_blue_vec,'LineStyle','-','LineWidth',2);
hold on
p4 = plot(datenum_vec,exp_mean_rzmc,'Color',dark_blue_vec,'LineStyle','--','LineWidth',2);
hold off
title('water budget timeseries')
datetick('x','dd-mmm')
xlabel('Day')
yyaxis left
ylabel('Daily avg water flux [mm]')
yyaxis right
ylabel('Daily avg soil moisture content [m^3 m^{-3}]')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,[save_path storm_name '_areaAvgBudget_landMask'],'png')

mean_totevptrns_mat = [control_mean_totevptrns; exp_mean_totevptrns];
mean_runoff_mat = [control_mean_runoff; exp_mean_runoff];
%mean_mc_mat = [control_mean_sfmc; exp_mean_sfmc; control_mean_rzmc; exp_mean_rzmc];
mean_mc_mat = [control_mean_sfmc; exp_mean_sfmc];

figure
tiledlayout(2,2)
ax1 = nexttile;
b = bar(datenum_vec,mean_precip_mat);
set(b(1),'FaceColor',gray_mat(1,:));
set(b(2),'FaceColor',gray_mat(2,:));
title('daily total precipitation [mm]')
datetick('x','dd-mmm')
xlabel('Day')
ylabel('precipitation [mm]')
ax2=nexttile;
e = bar(datenum_vec,mean_totevptrns_mat);
set(e(1),'FaceColor',green_vec);
set(e(2),'FaceColor',green_vec,'FaceAlpha',0.5);
title('daily total evapotranspiration [mm]')
datetick('x','dd-mmm')
xlabel('Day')
ylabel('evapotranspiration [mm]')
ax3=nexttile;
r = bar(datenum_vec,mean_runoff_mat);
set(r(1),'FaceColor',red_vec);
set(r(2),'FaceColor',red_vec,'FaceAlpha',0.5);
title('daily total runoff [mm]')
datetick('x','dd-mmm')
xlabel('Day')
ylabel('runoff [mm]')
ax4=nexttile;
s = bar(datenum_vec,mean_mc_mat);
set(s(1),'FaceColor',blue_vec);
set(s(2),'FaceColor',blue_vec,'FaceAlpha',0.5);
%set(s(3),'FaceColor',dark_blue_vec);
%set(s(4),'FaceColor',dark_blue_vec,'FaceAlpha',0.5);
title('daily avg soil moisture [m^3 m^{-3}]')
datetick('x','dd-mmm')
xlabel('Day')
ylabel('soil moisture [m^3 m^{-3}]')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 35 15])
saveas(gcf,[save_path storm_name '_areaAvgBudget_Bar_landMask'],'png')



mean_flux_diff_mat = [diff_mean_precip;...
                      diff_mean_totevptrns;...
                      diff_mean_runoff];

mean_water_content_diff_mat = [diff_mean_sfmc;...
                               diff_mean_rzmc];


figure 
tiledlayout(2,1)
ax1 = nexttile;
b = bar(datenum_vec,mean_flux_diff_mat);
set(b(1),'FaceColor',gray_mat(1,:));
set(b(2),'FaceColor',green_vec);
set(b(3),'FaceColor',red_vec);
title('Daily total water flux difference (EXP - CNTRL)')
datetick('x','dd-mmm')
xlabel('Day')
ylabel('\Delta water flux [mm]')
ax2 = nexttile;
p = bar(datenum_vec,mean_water_content_diff_mat);
set(p(1),'FaceColor',blue_vec);
set(p(2),'FaceColor',dark_blue_vec);
title('Daily avg. water content difference (EXP - CNTRL)')
datetick('x','dd-mmm')
xlabel('Day')
ylabel('\Delta soil moisture [m^3 m^{-3}]')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,[save_path storm_name '_areaAvgBudgetDiff_landMask'],'png')

diff_mean_sfmc_mm = diff_mean_sfmc * 50;
diff_mean_rzmc_mm = diff_mean_rzmc * 1000;

mean_water_content_diff_mat_mm = [diff_mean_sfmc_mm;...
                               diff_mean_rzmc_mm];

figure
tiledlayout(2,1)
ax1 = nexttile;
b = bar(datenum_vec,mean_flux_diff_mat);
set(b(1),'FaceColor',gray_mat(1,:));
set(b(2),'FaceColor',green_vec);
set(b(3),'FaceColor',red_vec);
title('Daily total water flux difference (EXP - CNTRL)')
datetick('x','dd-mmm')
xlabel('Day')
ylabel('\Delta water flux [mm]')
ax2 = nexttile;
p = bar(datenum_vec,mean_water_content_diff_mat_mm);
set(p(1),'FaceColor',blue_vec);
set(p(2),'FaceColor',dark_blue_vec);
title('Daily avg. water content difference (EXP - CNTRL)')
datetick('x','dd-mmm')
xlabel('Day')
ylabel('\Delta soil water content [mm]')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 25 15])
saveas(gcf,[save_path storm_name '_areaAvgBudgetDiff_allMM_landMask'],'png')

%test_diff_sfmc = test_exp_sfmc - test_control_sfmc;
%
%% define colorbars and limits
%cmap_mat = jet(19);
%cmap_mat = flipud(cmap_mat);
%cmap_mat(ceil(length(cmap_mat)/2),:) = 0.7*ones(1,3);
%cmap = cmap_mat;
%
%s1 = quantile(test_diff_sfmc(:),0.99);
%s2 = quantile(test_diff_sfmc(:),0.01);
%caxis_min_sfmc = -0.09; %-1*max(abs(s1),abs(s2));
%caxis_max_sfmc = 0.09; %max(abs(s1),abs(s2));
%
%c_length = caxis_max_sfmc/0.01;
%c_length = c_length*2;
%
%cmap_mat_sfmc = jet(c_length);
%cmap_mat_sfmc = flipud(cmap_mat_sfmc);
%cmap_mat_sfmc(length(cmap_mat_sfmc)/2:length(cmap_mat_sfmc)/2+1,:) = 0.85*ones(2,3);
%cmap_sfmc = cmap_mat_sfmc;
%
%
%    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_vec,lon_vec,test_diff_sfmc')
%    caxis([caxis_min_sfmc caxis_max_sfmc]);
%    set(gcf,'Colormap',cmap_sfmc);
%    plotm(coastlat, coastlon,'k')
%    set(gca,'XLim',lonlim);
%    set(gca,'YLim',latlim);
%    tightmap
%    title(['\Delta SFMC (EXP - CTRL) [m^{3} m^{-3}] 2019-03-12 0z'])
%    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 '_diag_SFMC_diff_20190312_0z'],'png')

