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/' storm_name '/'];
control_name = 'LADAS_control_4D_new';
exp_dir = ['/discover/nobackup/projects/smaptc/jkolassa/archive/' storm_name '/'];
exp_name = 'Idai_SMAP_4D';

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

start_year = 2019;
start_month = 3;
start_day = 4;

end_year = 2019;
end_month = 3;
end_day = 15;

land_mask = 'off'; % 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;

%---------------------------------------------
% 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

                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/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/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'));

                   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;


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

                   control_sfmc_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   control_rzmc_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
%                   control_evpintr_tavg24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
%                   control_evpsoil_tavg24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
%                   control_evptrns_tavg24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
                    control_totevptrns_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                    control_precip_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                    control_runoff_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                    control_totevptrns_ttot24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                    control_precip_ttot24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                    control_runoff_ttot24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));

                   exp_sfmc_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   exp_rzmc_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
%                   exp_evpintr_tavg24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
%                   exp_evpsoil_tavg24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
%                   exp_evptrns_tavg24_mat = NaN*ones(length(lon_vec),length(lat_vec),sum(monthLength(start_month:end_month)));
                    exp_totevptrns_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                    exp_precip_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                    exp_runoff_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                    exp_totevptrns_ttot24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                    exp_precip_ttot24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                    exp_runoff_ttot24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));

                   control_sfmc_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   control_rzmc_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   control_evpintr_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   control_evpsoil_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   control_evptrns_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   control_totevptrns_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   control_precip_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   control_runoff_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);

                   exp_sfmc_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   exp_rzmc_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   exp_evpintr_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   exp_evpsoil_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   exp_evptrns_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   exp_totevptrns_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   exp_precip_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   exp_runoff_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                end % if

                % read hourly data

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

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

                % 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')); 

                lat_index_c = (control_ana_lat<lat_range(1)) + (control_ana_lat>lat_range(2));
                lat_index_c = lat_index_c == 0;
                lon_index_c = (control_ana_lon<lon_range(1)) + (control_ana_lon>lon_range(2));
                lon_index_c = lon_index_c == 0;

                control_u_mat = control_u_mat(lon_index_c,lat_index_c);
                control_v_mat = control_v_mat(lon_index_c,lat_index_c);

                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'));

                lat_index_e = (exp_ana_lat<lat_range(1)) + (exp_ana_lat>lat_range(2));
                lat_index_e = lat_index_e == 0;
                lon_index_e = (exp_ana_lon<lon_range(1)) + (exp_ana_lon>lon_range(2));
                lon_index_e = lon_index_e == 0;

                exp_u_mat = exp_u_mat(lon_index_e,lat_index_e);
                exp_v_mat = exp_v_mat(lon_index_e,lat_index_e);

                [exp_ana_lat_mat,exp_ana_lon_mat] = meshgrid(exp_ana_lat(lat_index_e),exp_ana_lon(lon_index_e));
                control_u_interp_mat = interp2(control_ana_lat(lat_index_c),control_ana_lon(lon_index_c),control_u_mat,exp_ana_lat_mat,exp_ana_lon_mat);
                control_v_interp_mat = interp2(control_ana_lat(lat_index_c),control_ana_lon(lon_index_c),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 = (exp_ws_mat>=17) & (exp_frland>=0.5);
               %control_land_mask = (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==0) = NaN;
                    control_totevptrns_hour_mat(:,:,t+1) = tmp_mat;

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

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

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

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

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

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

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

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

                    tmp_mat = exp_rzmc_hour_mat(:,:,t+1);
                    tmp_mat(exp_land_mask==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;

            % convert precipitation and runoff from kg m^-2 s^-1 to mm/day

            control_precip_hour_mat(:,:,t+1) = control_precip_hour_mat(:,:,t+1) * 86400;
            control_runoff_hour_mat(:,:,t+1) = control_runoff_hour_mat(:,:,t+1) * 86400;

            exp_precip_hour_mat(:,:,t+1) = exp_precip_hour_mat(:,:,t+1) * 86400;
            exp_runoff_hour_mat(:,:,t+1) = exp_runoff_hour_mat(:,:,t+1) * 86400;

            % 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_evpintr_tavg24_mat(:,:,mat_index) = nanmean(control_evpintr_hour_mat,3);
%            control_evpsoil_tavg24_mat(:,:,mat_index) = nanmean(control_evpsoil_hour_mat,3);
%            control_evptrns_tavg24_mat(:,:,mat_index) = nanmean(control_evptrns_hour_mat,3);
            control_totevptrns_tavg24_mat(:,:,mat_index) = nanmean(control_totevptrns_hour_mat,3);
            control_precip_tavg24_mat(:,:,mat_index) = nanmean(control_precip_hour_mat,3);
            control_runoff_tavg24_mat(:,:,mat_index) = nanmean(control_runoff_hour_mat,3);
            control_totevptrns_ttot24_mat(:,:,mat_index) = sum(tmp_cntrl_totevptrns_mat,3,'omitnan');
            control_precip_ttot24_mat(:,:,mat_index) = sum(tmp_cntrl_precip_mat,3,'omitnan');
            control_runoff_ttot24_mat(:,:,mat_index) = sum(tmp_cntrl_runoff_mat,3,'omitnan');

            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_evpintr_tavg24_mat(:,:,mat_index) = nanmean(exp_evpintr_hour_mat,3);
%            exp_evpsoil_tavg24_mat(:,:,mat_index) = nanmean(exp_evpsoil_hour_mat,3);
%            exp_evptrns_tavg24_mat(:,:,mat_index) = nanmean(exp_evptrns_hour_mat,3);
            exp_totevptrns_tavg24_mat(:,:,mat_index) = nanmean(exp_totevptrns_hour_mat,3);
            exp_precip_tavg24_mat(:,:,mat_index) = nanmean(exp_precip_hour_mat,3);
            exp_runoff_tavg24_mat(:,:,mat_index) = nanmean(exp_runoff_hour_mat,3);
            exp_totevptrns_ttot24_mat(:,:,mat_index) = sum(tmp_exp_totevptrns_mat,3,'omitnan');
            nanmat = isnan(tmp_exp_totevptrns_mat);
            nan_index = (sum(nanmat,3)==size(tmp_exp_totevptrns_mat,3));
            tmp_mat = exp_totevptrns_ttot24_mat(:,:,mat_index);
            tmp_mat(nan_index) = NaN;
            exp_totevptrns_ttot24_mat(:,:,mat_index) = tmp_mat;
            exp_precip_ttot24_mat(:,:,mat_index) = sum(tmp_exp_precip_mat,3,'omitnan');
            exp_runoff_ttot24_mat(:,:,mat_index) = sum(tmp_exp_runoff_mat,3,'omitnan');

            % reset hourly matrices

            control_sfmc_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            control_rzmc_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            control_evpintr_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            control_evpsoil_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            control_evptrns_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            control_totevptrns_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);   
            control_precip_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);         
            control_runoff_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);

            exp_sfmc_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            exp_rzmc_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            exp_evpintr_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            exp_evpsoil_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            exp_evptrns_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            exp_totevptrns_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            exp_precip_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            exp_runoff_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),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_evpintr_tavg24_mat(:,:,(mat_index+1):end) = [];
%exp_evpintr_tavg24_mat(:,:,(mat_index+1):end) = [];
%control_evpsoil_tavg24_mat(:,:,(mat_index+1):end) = [];
%exp_evpsoil_tavg24_mat(:,:,(mat_index+1):end) = [];
%control_evptrns_tavg24_mat(:,:,(mat_index+1):end) = [];
%exp_evptrns_tavg24_mat(:,:,(mat_index+1):end) = [];
control_totevptrns_tavg24_mat(:,:,(mat_index+1):end) = [];
exp_totevptrns_tavg24_mat(:,:,(mat_index+1):end) = [];
control_precip_tavg24_mat(:,:,(mat_index+1):end) = [];
exp_precip_tavg24_mat(:,:,(mat_index+1):end) = [];
control_runoff_tavg24_mat(:,:,(mat_index+1):end) = [];
exp_runoff_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) = [];


% 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_evpintr_tavg24_mat = exp_evpintr_tavg24_mat - control_evpintr_tavg24_mat;
%diff_evpsoil_tavg24_mat = exp_evpsoil_tavg24_mat - control_evpsoil_tavg24_mat;
%diff_evptrns_tavg24_mat = exp_evptrns_tavg24_mat - control_evptrns_tavg24_mat;
diff_totevptrns_tavg24_mat = exp_totevptrns_tavg24_mat - control_totevptrns_tavg24_mat;
diff_precip_tavg24_mat = exp_precip_tavg24_mat - control_precip_tavg24_mat;
diff_runoff_tavg24_mat = exp_runoff_tavg24_mat - control_runoff_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;

diff_sfmc_tot_mat = nanmean(diff_sfmc_tavg24_mat,3);
diff_rzmc_tot_mat = nanmean(diff_rzmc_tavg24_mat,3);

% 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;

tmp_mat = diff_sfmc_tavg24_mat;

s1 = quantile(tmp_mat(:),0.99);
s2 = quantile(tmp_mat(:),0.01);
caxis_min_sfmc = -1*max(abs(s1),abs(s2));
caxis_min_sfmc = round(caxis_min_sfmc,2);
caxis_max_sfmc = max(abs(s1),abs(s2));
caxis_max_sfmc = round(caxis_max_sfmc,2);

c_length = caxis_max_sfmc/0.01;
c_length = c_length*2;

cmap_mat_sfmc = cmocean('balance',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;

tmp_mat = diff_rzmc_tavg24_mat;

s1 = quantile(tmp_mat(:),0.99);
s2 = quantile(tmp_mat(:),0.01);
caxis_min_rzmc = -1*max(abs(s1),abs(s2));
caxis_min_rzmc = round(caxis_min_rzmc,2);
caxis_max_rzmc = max(abs(s1),abs(s2));
caxis_max_rzmc = round(caxis_max_rzmc,2);

c_length = caxis_max_rzmc/0.005;
c_length = c_length*2;

cmap_mat_rzmc = cmocean('balance',c_length);
cmap_mat_rzmc = flipud(cmap_mat_rzmc);
cmap_mat_rzmc(length(cmap_mat_rzmc)/2:length(cmap_mat_rzmc)/2+1,:) = 0.85*ones(2,3);
cmap_rzmc = cmap_mat_rzmc;

%tmp_mat = diff_evpintr_tavg24_mat(lon_index,lat_index,:);
%
%s1 = quantile(tmp_mat(:),0.99);
%s2 = quantile(tmp_mat(:),0.01);
%caxis_min_evpintr = -1*max(abs(s1),abs(s2));
%caxis_min_evpintr = round(caxis_min_evpintr,1);
%caxis_max_evpintr = max(abs(s1),abs(s2));
%caxis_max_evpintr = round(caxis_max_evpintr,1);
%
%c_length = caxis_max_evpintr/0.1;
%c_length = c_length*2;
%
%cmap_mat_evpintr = cmocean('balance',c_length);
%cmap_mat_evpintr = flipud(cmap_mat_evpintr);
%cmap_mat_evpintr(length(cmap_mat_evpintr)/2:length(cmap_mat_evpintr)/2+1,:) = 0.85*ones(2,3);
%cmap_evpintr = cmap_mat_evpintr;
%
%tmp_mat = diff_evpsoil_tavg24_mat(lon_index,lat_index,:);
%
%s1 = quantile(tmp_mat(:),0.99);
%s2 = quantile(tmp_mat(:),0.01);
%caxis_min_evpsoil = -1*max(abs(s1),abs(s2));
%caxis_min_evpsoil = round(caxis_min_evpsoil,1);
%caxis_max_evpsoil = max(abs(s1),abs(s2));
%caxis_max_evpsoil = round(caxis_max_evpsoil,1);
%
%c_length = caxis_max_evpsoil/0.1;
%c_length = c_length*2;
%
%cmap_mat_evpsoil = cmocean('balance',c_length);
%cmap_mat_evpsoil = flipud(cmap_mat_evpsoil);
%cmap_mat_evpsoil(length(cmap_mat_evpsoil)/2:length(cmap_mat_evpsoil)/2+1,:) = 0.85*ones(2,3);
%cmap_evpsoil = cmap_mat_evpsoil;
%
%tmp_mat = diff_evptrns_tavg24_mat(lon_index,lat_index,:);
%
%s1 = quantile(tmp_mat(:),0.99);
%s2 = quantile(tmp_mat(:),0.01);
%caxis_min_evptrns = -1*max(abs(s1),abs(s2));
%caxis_min_evptrns = round(caxis_min_evptrns,1);
%caxis_max_evptrns = max(abs(s1),abs(s2));
%caxis_max_evptrns = round(caxis_max_evptrns,1);
%
%c_length = caxis_max_evptrns/0.1;
%c_length = c_length*2;
%
%cmap_mat_evptrns = cmocean('balance',c_length);
%cmap_mat_evptrns = flipud(cmap_mat_evptrns);
%cmap_mat_evptrns(length(cmap_mat_evptrns)/2:length(cmap_mat_evptrns)/2+1,:) = 0.85*ones(2,3);
%cmap_evptrns = cmap_mat_evptrns;

tmp_mat = diff_totevptrns_ttot24_mat;

s1 = quantile(tmp_mat(:),0.99);
s2 = quantile(tmp_mat(:),0.01);
caxis_min_totevptrns = -1*max(abs(s1),abs(s2));
caxis_min_totevptrns = round(caxis_min_totevptrns,1);
caxis_max_totevptrns = max(abs(s1),abs(s2));
caxis_max_totevptrns = round(caxis_max_totevptrns,1);

c_length = caxis_max_totevptrns/0.2;
c_length = c_length*2;

cmap_mat_totevptrns = cmocean('balance',c_length);
cmap_mat_totevptrns = flipud(cmap_mat_totevptrns);
cmap_mat_totevptrns(length(cmap_mat_totevptrns)/2:length(cmap_mat_totevptrns)/2+1,:) = 0.85*ones(2,3);
cmap_totevptrns = cmap_mat_totevptrns;

tmp_mat = diff_precip_ttot24_mat;
tmp_mat(tmp_mat==0) = NaN;

s1 = quantile(tmp_mat(:),0.999);
s2 = quantile(tmp_mat(:),0.001);
%s1 = max(tmp_mat(:));
%s2 = min(tmp_mat(:));
caxis_min_precip = -1*max(abs(s1),abs(s2));
caxis_min_precip = round(caxis_min_precip);
caxis_max_precip = max(abs(s1),abs(s2));
caxis_max_precip = round(caxis_max_precip);

c_length = caxis_max_precip/5;
c_length = c_length*2;

cmap_mat_precip = cmocean('balance',c_length);
cmap_mat_precip = flipud(cmap_mat_precip);
cmap_mat_precip(length(cmap_mat_precip)/2:length(cmap_mat_precip)/2+1,:) = 0.85*ones(2,3);
cmap_precip = cmap_mat_precip;

tmp_mat = control_sfmc_tavg24_mat;

s1 = quantile(tmp_mat(:),0.99);
s2 = quantile(tmp_mat(:),0.01);
caxis_min_sfmc_raw = -1*max(abs(s1),abs(s2));
caxis_min_sfmc_raw = round(caxis_min_sfmc_raw,1);
caxis_max_sfmc_raw = max(abs(s1),abs(s2));
caxis_max_sfmc_raw = round(caxis_max_sfmc_raw,1);

c_length = caxis_max_sfmc_raw/0.1;
c_length = c_length*2;

cmap_mat_sfmc_raw = cmocean('balance',c_length);
cmap_mat_sfmc_raw = flipud(cmap_mat_sfmc_raw);
cmap_mat_sfmc_raw(length(cmap_mat_sfmc_raw)/2:length(cmap_mat_sfmc_raw)/2+1,:) = 0.85*ones(2,3);
cmap_sfmc_raw = cmap_mat_sfmc_raw;

tmp_mat = control_rzmc_tavg24_mat;

s1 = quantile(tmp_mat(:),0.99);
s2 = quantile(tmp_mat(:),0.01);
caxis_min_rzmc_raw = -1*max(abs(s1),abs(s2));
caxis_min_rzmc_raw = round(caxis_min_rzmc_raw,1);
caxis_max_rzmc_raw = max(abs(s1),abs(s2));
caxis_max_rzmc_raw = round(caxis_max_rzmc_raw,1);

c_length = caxis_max_rzmc_raw/0.1;
c_length = c_length*2;

cmap_mat_rzmc_raw = cmocean('balance',c_length);
cmap_mat_rzmc_raw = flipud(cmap_mat_rzmc_raw);
cmap_mat_rzmc_raw(length(cmap_mat_rzmc_raw)/2:length(cmap_mat_rzmc_raw)/2+1,:) = 0.85*ones(2,3);
cmap_rzmc_raw = cmap_mat_rzmc_raw;



%v_sm = VideoWriter([save_path storm_name '_SFMC_diff_tavg24']);
%v_sm.FrameRate = 4;
%open(v_sm);
%
%v_rz = VideoWriter([save_path storm_name '_RZMC_diff_tavg24']);
%v_rz.FrameRate = 4;
%open(v_rz);
%
%v_evpintr = VideoWriter([save_path storm_name '_EVPINTR_diff_tavg24']);
%v_evpintr.FrameRate = 2;
%open(v_evpintr);
%
%v_evpsoil = VideoWriter([save_path storm_name '_EVPSOIL_diff_tavg24']);
%v_evpsoil.FrameRate = 2;
%open(v_evpsoil);
%
%v_evptrns = VideoWriter([save_path storm_name '_EVPTRNS_diff_tavg24']);
%v_evptrns.FrameRate = 2;
%open(v_evptrns);
%
%v_totevptrns = VideoWriter([save_path storm_name '_TOTEVPTRNS_diff_tavg24']);
%v_totevptrns.FrameRate = 2;
%open(v_totevptrns);

control_mean_sfmc = NaN*ones(1,length(datenum_vec));
control_mean_rzmc = NaN*ones(1,length(datenum_vec));
control_mean_totevptrns = NaN*ones(1,length(datenum_vec));
control_mean_precip = NaN*ones(1,length(datenum_vec));
control_mean_runoff = NaN*ones(1,length(datenum_vec));

exp_mean_sfmc = NaN*ones(1,length(datenum_vec));
exp_mean_rzmc = NaN*ones(1,length(datenum_vec));
exp_mean_totevptrns = NaN*ones(1,length(datenum_vec));
exp_mean_precip = NaN*ones(1,length(datenum_vec));
exp_mean_runoff = NaN*ones(1,length(datenum_vec));

diff_mean_sfmc = NaN*ones(1,length(datenum_vec));
diff_mean_rzmc = NaN*ones(1,length(datenum_vec));
diff_mean_totevptrns = NaN*ones(1,length(datenum_vec));
diff_mean_precip = NaN*ones(1,length(datenum_vec));
diff_mean_runoff = NaN*ones(1,length(datenum_vec));


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,diff_sfmc_tavg24_mat(:,:,i)')
%    caxis([caxis_min_sfmc caxis_max_sfmc]);
%    set(gcf,'Colormap',cmap_sfmc);
%    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(['\Delta SFMC (EXP - CTRL) [m^{3} m^{-3}] ' 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')

%    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,control_sfmc_tavg24_mat(:,:,i)')
%    caxis([caxis_min_sfmc_raw caxis_max_sfmc_raw]);
%    set(gcf,'Colormap',cmap_sfmc_raw);
%    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(['SFMC CTRL [m^{3} m^{-3}] ' 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_CNTRL_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')
%
%    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,exp_sfmc_tavg24_mat(:,:,i)')
%    caxis([caxis_min_sfmc_raw caxis_max_sfmc_raw]);
%    set(gcf,'Colormap',cmap_sfmc_raw);
%    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(['SFMC EXP [m^{3} m^{-3}] ' 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_EXP_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')


  %  frame = getframe(gcf);
  %  writeVideo(v_sm,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,diff_rzmc_tavg24_mat(:,:,i)')
%    caxis([caxis_min_rzmc caxis_max_rzmc]);
%    set(gcf,'Colormap',cmap_rzmc);
%    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(['\Delta RZMC (EXP - CTRL) [m^{3} m^{-3}] ' 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_RZMC_diff_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')

%    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,control_rzmc_tavg24_mat(:,:,i)')
%    caxis([caxis_min_rzmc_raw caxis_max_rzmc_raw]);
%    set(gcf,'Colormap',cmap_rzmc_raw);
%    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(['RZMC CTRL [m^{3} m^{-3}] ' 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_RZMC_CNTRL_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')
%
%    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,exp_rzmc_tavg24_mat(:,:,i)')
%    caxis([caxis_min_rzmc_raw caxis_max_rzmc_raw]);
%    set(gcf,'Colormap',cmap_rzmc_raw);
%    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(['RZMC EXP [m^{3} m^{-3}] ' 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_RZMC_EXP_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')


  %  frame = getframe(gcf);
  %  writeVideo(v_rz,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,diff_evpintr_tavg24_mat(:,:,i)')
%    caxis([caxis_min_evpintr caxis_max_evpintr]);
%    set(gcf,'Colormap',cmap_evpintr);
%    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(['\Delta EVPINTR (EXP - CTRL) [mm day^{-1}] ' 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_EVPINTR_diff_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')
%
%  %  frame = getframe(gcf);
%  %  writeVideo(v_evpintr,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,diff_evpsoil_tavg24_mat(:,:,i)')
%    caxis([caxis_min_evpsoil caxis_max_evpsoil]);
%    set(gcf,'Colormap',cmap_evpsoil);
%    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(['\Delta EVPSOIL (EXP - CTRL) [mm day^{-1}] ' 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_EVPSOIL_diff_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')
%
%  %  frame = getframe(gcf);
%  %  writeVideo(v_evpsoil,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,diff_evptrns_tavg24_mat(:,:,i)')
%    caxis([caxis_min_evptrns caxis_max_evptrns]);
%    set(gcf,'Colormap',cmap_evptrns);
%    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);
%    seggt(gca,'YLim',latlim);
%    tightmap
%    title(['\Delta EVPTRNS (EXP - CTRL) [mm day^{-1}] ' 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_EVPTRNS_diff_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')

 %   frame = getframe(gcf);
 %   writeVideo(v_evptrns,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,diff_totevptrns_ttot24_mat(:,:,i)')
%    caxis([caxis_min_totevptrns caxis_max_totevptrns]);
%    set(gcf,'Colormap',cmap_totevptrns);
%    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(['\Delta total daily evapotranspiration (EXP - CTRL) [mm] ' 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_TOTEVPTRNS_diff_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')
%
%  %  frame = getframe(gcf);
%  %  writeVideo(v_totevptrns,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,diff_precip_ttot24_mat(:,:,i)')
%    caxis([caxis_min_precip caxis_max_precip]);
%    set(gcf,'Colormap',cmap_precip);
%    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(['\Delta total daily precipitation (EXP - CTRL) [mm] ' 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_PRECIP_diff_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')



  % 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


avg_sfmc_diff_plot_mat = nanmean(diff_sfmc_tavg24_mat,3);
avg_rzmc_diff_plot_mat = nanmean(diff_rzmc_tavg24_mat,3);

tmp_mat = avg_sfmc_diff_plot_mat;

s1 = quantile(tmp_mat(:),0.99);
s2 = quantile(tmp_mat(:),0.01);
caxis_min_sfmc_avg = -1*max(abs(s1),abs(s2));
caxis_min_sfmc_avg = round(caxis_min_sfmc_avg,2);
caxis_max_sfmc_avg = max(abs(s1),abs(s2));
caxis_max_sfmc_avg = round(caxis_max_sfmc_avg,2);

c_length = caxis_max_sfmc_avg/0.01;
c_length = c_length*2;

cmap_mat_sfmc_avg = cmocean('balance',c_length);
cmap_mat_sfmc_avg = flipud(cmap_mat_sfmc_avg);
cmap_mat_sfmc_avg(length(cmap_mat_sfmc_avg)/2:length(cmap_mat_sfmc_avg)/2+1,:) = 0.85*ones(2,3);
cmap_sfmc_avg = cmap_mat_sfmc_avg;

tmp_mat = avg_rzmc_diff_plot_mat;

s1 = quantile(tmp_mat(:),0.99);
s2 = quantile(tmp_mat(:),0.01);
caxis_min_rzmc_avg = -1*max(abs(s1),abs(s2));
caxis_min_rzmc_avg = round(caxis_min_rzmc_avg,2);
caxis_max_rzmc_avg = max(abs(s1),abs(s2));
caxis_max_rzmc_avg = round(caxis_max_rzmc_avg,2);

c_length = caxis_max_rzmc_avg/0.01;
c_length = c_length*2;

cmap_mat_rzmc_avg = cmocean('balance',c_length);
cmap_mat_rzmc_avg = flipud(cmap_mat_rzmc_avg);
cmap_mat_rzmc_avg(length(cmap_mat_rzmc_avg)/2:length(cmap_mat_rzmc_avg)/2+1,:) = 0.85*ones(2,3);
cmap_rzmc_avg = cmap_mat_rzmc_avg;



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(lat_index),lon_vec(lon_index),avg_sfmc_diff_plot_mat')
caxis([caxis_min_sfmc_avg caxis_max_sfmc_avg]);
set(gcf,'Colormap',cmap_sfmc_avg);
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(['\Delta SFMC (SMAP\_ANA - CTRL\_ANA) [m^3 m^{-3}]'])
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_avg'],'png')

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(lat_index),lon_vec(lon_index),avg_rzmc_diff_plot_mat')
caxis([caxis_min_rzmc_avg caxis_max_rzmc_avg]);
set(gcf,'Colormap',cmap_rzmc_avg);
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(['\Delta RZMC (SMAP\_ANA - CTRL\_ANA) [m^3 m^{-3}]'])
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_RZMC_diff_avg'],'png')



%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];

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')

