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 = 2;
%start_day = 4;

start_year = 2019;
start_month = 3;
start_day = 1;

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

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

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

bt_time_str = '72h';   % string specifying how far back in time backtrajectories are computed that are used to create land_mask

%---------------------------------------------
% 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;
month_vec = ['Jan';'Feb';'Mar';'Apr';'May';,'Jun';'Jul';'Aug';'Sep';'Oct';'Nov';'Dec'];

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


% define latitude and longitude index for domain

coord_file = [control_dir control_name '/diag/Y' num2str(start_year,'%04i') '/M' num2str(start_month,'%02i') ...
             '/' control_name '.tavg1_2d_lnd_Nx.' num2str(start_year,'%04i') num2str(start_month,'%02i') ...
              num2str(start_day,'%02i') '_2130z.nc4'];

lat_vec = double(ncread(coord_file,'lat'));
lon_vec = double(ncread(coord_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;

coord_ana_file = [control_dir control_name '/ana/Y' num2str(start_year,'%04i') '/M' num2str(start_month,'%02i') '/' control_name ... 
                 '.ana.inst3d_met_p.' num2str(start_year,'%04i') num2str(start_month,'%02i') num2str(start_day,'%02i')...
                 '_2100z.nc4'];

lat_vec_coarse = double(ncread(coord_ana_file,'lat'));
lon_vec_coarse = double(ncread(coord_ana_file,'lon'));

lat_index_coarse = (lat_vec_coarse<lat_range(1)) + (lat_vec_coarse>lat_range(2));
lat_index_coarse = lat_index_coarse == 0;
lon_index_coarse = (lon_vec_coarse<lon_range(1)) + (lon_vec_coarse>lon_range(2));
lon_index_coarse = lon_index_coarse == 0;

% loop over years
for y = start_year:end_year

    disp('read data')

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

                land_mask_file = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/' ...
                                storm_name '/land_mask/land_mask_' d_str '-' month_vec(m,:) '-' y_str '_' ...
                                t_str  ':30:00.mat'];

                land_mask_file_control = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/' ...
                                storm_name '/land_mask/land_mask_interp_' bt_time_str '_control_' d_str '-' month_vec(m,:) '-' y_str '_' ...
                                '.mat'];

                land_mask_file_exp = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/backtrajectories/' ...
                                storm_name '/land_mask/land_mask_interp_' bt_time_str '_experiment_' d_str '-' month_vec(m,:) '-' y_str '_' ...
                                '.mat'];

                % on first time step read coordinate information and initialize matrices

                if (mat_index == 1) && (t == 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_prmc_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_totlh_ttot24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   control_totsh_ttot24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   control_KEint_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   control_QVint_tavg24_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_prmc_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_totlh_ttot24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   exp_totsh_ttot24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   exp_KEint_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   exp_QVint_tavg24_mat = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));

                   control_land_mask_heat_map_daily = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   exp_land_mask_heat_map_daily = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   control_land_mask_daily = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   exp_land_mask_daily = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   control_land_mask_time_daily = NaN*ones(sum(lon_index),sum(lat_index),sum(monthLength(start_month:end_month)));
                   exp_land_mask_time_daily = 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_prmc_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_totlh_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   control_totsh_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   control_KEint_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   control_QVint_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_prmc_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_totlh_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   exp_totsh_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   exp_KEint_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
                   exp_QVint_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);

                   exp_land_mask = NaN*ones(sum(lon_index),sum(lat_index),24);
                   control_land_mask = NaN*ones(sum(lon_index),sum(lat_index),24);

                end % if

                % read hourly data

                tmp_read_mat = double(ncread(control_diag_file,'SFMC'));
                control_sfmc_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(control_diag_file,'RZMC'));
                control_rzmc_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(control_diag_file,'PRMC'));
                control_prmc_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(control_diag_file,'EVPINTR'));
                control_evpintr_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(control_diag_file,'EVPSOIL'));
                control_evpsoil_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(control_diag_file,'EVPTRNS'));
                control_evptrns_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(control_diag_file,'LHLAND'));
                control_totlh_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(control_diag_file,'SHLAND'));
                control_totsh_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);

                tmp_read_mat = double(ncread(exp_diag_file,'SFMC'));
                exp_sfmc_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(exp_diag_file,'RZMC'));
                exp_rzmc_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(exp_diag_file,'PRMC'));
                exp_prmc_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(exp_diag_file,'EVPINTR'));
                exp_evpintr_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(exp_diag_file,'EVPSOIL'));
                exp_evpsoil_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(exp_diag_file,'EVPTRNS'));
                exp_evptrns_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(exp_diag_file,'LHLAND'));
                exp_totlh_hour_mat(:,:,t+1) = tmp_read_mat(lon_index,lat_index);
                tmp_read_mat = double(ncread(exp_diag_file,'SHLAND'));
                exp_totsh_hour_mat(:,:,t+1) = tmp_read_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>=150);
                tmp_read_mat = double(ncread(control_ana_file,'U'));
                control_u_mat = tmp_read_mat(lon_index_coarse,lat_index_coarse,control_level_index);
                tmp_read_mat = double(ncread(control_ana_file,'V'));
                control_v_mat = tmp_read_mat(lon_index_coarse,lat_index_coarse,control_level_index);
                tmp_read_mat = double(ncread(control_ana_file,'T'));
                control_temp_mat = tmp_read_mat(lon_index_coarse,lat_index_coarse,control_level_index);
                tmp_read_mat = double(ncread(control_ana_file,'QV'));
                control_qv_mat = tmp_read_mat(lon_index_coarse,lat_index_coarse,control_level_index);

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

                exp_hPa = double(ncread(exp_ana_file,'lev'));
                exp_level_index = find(exp_hPa>=600);
                tmp_read_mat = double(ncread(exp_ana_file,'U'));
                exp_u_mat = tmp_read_mat(lon_index,lat_index,control_level_index);
                tmp_read_mat = double(ncread(exp_ana_file,'V'));
                exp_v_mat = tmp_read_mat(lon_index,lat_index,control_level_index);
                tmp_read_mat = double(ncread(exp_ana_file,'T'));
                exp_temp_mat = tmp_read_mat(lon_index,lat_index,control_level_index);
                tmp_read_mat = double(ncread(exp_ana_file,'QV'));
                exp_qv_mat = tmp_read_mat(lon_index,lat_index,control_level_index);
                
                tmp_read_mat = double(ncread(exp_ana_file,'lat'));
                exp_ana_lat = tmp_read_mat(lat_index);
                tmp_read_mat = double(ncread(exp_ana_file,'lon'));
                exp_ana_lon = tmp_read_mat(lon_index);

                [exp_ana_lat_mat,exp_ana_lon_mat] = meshgrid(exp_ana_lat,exp_ana_lon);

                control_KE_lev = NaN*ones(sum(lon_index),sum(lat_index),length(exp_level_index));
                exp_KE_lev = NaN*ones(sum(lon_index),sum(lat_index),length(exp_level_index));
 
                control_QV_lev = NaN*ones(sum(lon_index),sum(lat_index),length(exp_level_index));
                exp_QV_lev = exp_qv_mat;

                for l = 1:length(exp_level_index)

                    control_u_interp_mat = interp2(control_ana_lat,control_ana_lon,control_u_mat(:,:,l),exp_ana_lat_mat,exp_ana_lon_mat);
                    control_v_interp_mat = interp2(control_ana_lat,control_ana_lon,control_v_mat(:,:,l),exp_ana_lat_mat,exp_ana_lon_mat);
                    control_temp_interp_mat = interp2(control_ana_lat,control_ana_lon,control_temp_mat(:,:,l),exp_ana_lat_mat,exp_ana_lon_mat);
                    control_QV_lev(:,:,l) = interp2(control_ana_lat,control_ana_lon,control_qv_mat(:,:,l),exp_ana_lat_mat,exp_ana_lon_mat);


                    control_rho_mat = 1000*control_hPa(l)./(287*control_temp_interp_mat); 
                    exp_rho_mat = 1000*exp_hPa(l)./(287*exp_temp_mat(:,:,l));

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

                    control_KE_tmp = 0.5*control_rho_mat.*(control_ws_mat.^2);
                    exp_KE_tmp = 0.5*exp_rho_mat.*(exp_ws_mat.^2);

                    % remove points outside of tropical cyclone (wind speed exceeds 17 m/s)

                    TC_index_control = control_ws_mat >= 17;
                    TC_index_exp = exp_ws_mat >= 17;

                    control_KE_tmp(TC_index_control==0) = NaN;
                    exp_KE_tmp(TC_index_exp==0) = NaN;

                    control_KE_lev(:,:,l) = control_KE_tmp;
                    exp_KE_lev(:,:,l) = exp_KE_tmp;

                    control_QV_tmp = control_QV_lev(:,:,l);
                    exp_QV_tmp = exp_QV_lev(:,:,l);

                    control_QV_tmp(TC_index_control==0) = NaN;
                    exp_QV_tmp(TC_index_exp==0) = NaN;

                    control_QV_lev(:,:,l) = control_QV_tmp;
                    exp_QV_lev(:,:,l) = exp_QV_tmp;

               
                end % l 

                control_KEint_hour_mat(:,:,t+1) = sum(control_KE_lev,3,'omitnan');
                nan_mat = isnan(control_KE_lev);
                nan_sum = sum(nan_mat,3);
                nan_index = (nan_sum == length(exp_level_index));
                control_tmp = control_KEint_hour_mat(:,:,t+1);
                control_tmp(nan_index) = NaN;
                control_KEint_hour_mat(:,:,t+1) = control_tmp; 

                exp_KEint_hour_mat(:,:,t+1) = sum(exp_KE_lev,3,'omitnan');
                nan_mat = isnan(exp_KE_lev);
                nan_sum = sum(nan_mat,3);
                nan_index = (nan_sum == length(exp_level_index));
                exp_tmp = exp_KEint_hour_mat(:,:,t+1);
                exp_tmp(nan_index) = NaN;
                exp_KEint_hour_mat(:,:,t+1) = exp_tmp;

                control_QVint_hour_mat(:,:,t+1) = sum(control_QV_lev,3,'omitnan');
                nan_mat = isnan(control_QV_lev);
                nan_sum = sum(nan_mat,3);
                nan_index = (nan_sum == length(exp_level_index));
                control_tmp = control_QVint_hour_mat(:,:,t+1);
                control_tmp(nan_index) = NaN;
                control_QVint_hour_mat(:,:,t+1) = control_tmp;

                exp_QVint_hour_mat(:,:,t+1) = sum(exp_QV_lev,3,'omitnan');
                nan_mat = isnan(exp_QV_lev);
                nan_sum = sum(nan_mat,3);
                nan_index = (nan_sum == length(exp_level_index));
                exp_tmp = exp_QVint_hour_mat(:,:,t+1);
                exp_tmp(nan_index) = NaN;
                exp_QVint_hour_mat(:,:,t+1) = exp_tmp;
  

                % read land mask
                if (exist(land_mask_file_control))
                land_mask_control = load(land_mask_file_control);
                land_mask_exp = load(land_mask_file_exp);
                end

                % 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_flag,'on')) & (exist(land_mask_file_control))
           
                    tmp_mat = control_totevptrns_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_control.land_mask_mat(lon_index,lat_index)==0) = NaN;
                    control_totevptrns_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = exp_totevptrns_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_exp.land_mask_mat(lon_index,lat_index)==0) = NaN;
                    exp_totevptrns_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = control_sfmc_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_control.land_mask_mat(lon_index,lat_index)==0) = NaN; 
                    control_sfmc_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = exp_sfmc_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_exp.land_mask_mat(lon_index,lat_index)==0) = NaN; 
                    exp_sfmc_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = control_rzmc_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_control.land_mask_mat(lon_index,lat_index)==0) = NaN;
                    control_rzmc_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = exp_rzmc_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_exp.land_mask_mat(lon_index,lat_index)==0) = NaN;
                    exp_rzmc_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = control_prmc_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_control.land_mask_mat(lon_index,lat_index)==0) = NaN;
                    control_prmc_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = exp_prmc_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_exp.land_mask_mat(lon_index,lat_index)==0) = NaN;
                    exp_prmc_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = control_totlh_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_control.land_mask_mat(lon_index,lat_index)==0) = NaN;
                    control_totlh_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = exp_totlh_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_exp.land_mask_mat(lon_index,lat_index)==0) = NaN;
                    exp_totlh_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = control_totsh_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_control.land_mask_mat(lon_index,lat_index)==0) = NaN;
                    control_totsh_hour_mat(:,:,t+1) = tmp_mat;

                    tmp_mat = exp_totsh_hour_mat(:,:,t+1);
                    tmp_mat(land_mask_exp.land_mask_mat(lon_index,lat_index)==0) = NaN;
                    exp_totsh_hour_mat(:,:,t+1) = tmp_mat;

                 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_totlh_mat = control_totlh_hour_mat*0.0036*0.408;
            tmp_exp_totlh_mat = exp_totlh_hour_mat*0.0036*0.408;
            tmp_cntrl_totsh_mat = control_totsh_hour_mat*0.0036*0.408;
            tmp_exp_totsh_mat = exp_totsh_hour_mat*0.0036*0.408;

            % 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;
            control_totlh_hour_mat = control_totlh_hour_mat*0.0864*0.408;
            exp_totlh_hour_mat = exp_totlh_hour_mat*0.0864*0.408;            
            control_totsh_hour_mat = control_totsh_hour_mat*0.0864*0.408;
            exp_totsh_hour_mat = exp_totsh_hour_mat*0.0864*0.408;

            % log date

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

            % compute daily averages
            if (exist(land_mask_file_control))
            control_land_mask_heat_map_daily(:,:,mat_index) = land_mask_control.land_mask_heat_mat(lon_index,lat_index);
            exp_land_mask_heat_map_daily(:,:,mat_index) = land_mask_exp.land_mask_heat_mat(lon_index,lat_index);
            control_land_mask_daily(:,:,mat_index) = land_mask_control.land_mask_mat(lon_index,lat_index);
            exp_land_mask_daily(:,:,mat_index) = land_mask_exp.land_mask_mat(lon_index,lat_index);
            control_land_mask_time_daily(:,:,mat_index) = land_mask_control.land_mask_time_mat(lon_index,lat_index);
            exp_land_mask_time_daily(:,:,mat_index) = land_mask_exp.land_mask_time_mat(lon_index,lat_index);
            end

            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_prmc_tavg24_mat(:,:,mat_index) = nanmean(control_prmc_hour_mat,3);
            control_KEint_tavg24_mat(:,:,mat_index) = nanmean(control_KEint_hour_mat,3);
            control_QVint_tavg24_mat(:,:,mat_index) = nanmean(control_QVint_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));
            nan_index = isnan(control_rzmc_tavg24_mat(:,:,mat_index));
            tmp = control_totevptrns_ttot24_mat(:,:,mat_index);
            tmp(nan_index) = NaN;
            control_totevptrns_ttot24_mat(:,:,mat_index) = tmp;

            control_totlh_ttot24_mat(:,:,mat_index) = sum(tmp_cntrl_totlh_mat,3,'omitnan');
            %nan_index = (sum(isnan(tmp_cntrl_totlh_mat))==size(tmp_cntrl_totlh_mat,3));
            nan_index = isnan(control_rzmc_tavg24_mat(:,:,mat_index));
            tmp = control_totlh_ttot24_mat(:,:,mat_index);
            tmp(nan_index) = NaN;
            control_totlh_ttot24_mat(:,:,mat_index) = tmp;

            control_totsh_ttot24_mat(:,:,mat_index) = sum(tmp_cntrl_totsh_mat,3,'omitnan');
           % nan_index = (sum(isnan(tmp_cntrl_totsh_mat))==size(tmp_cntrl_totsh_mat,3));
            nan_index = isnan(control_rzmc_tavg24_mat(:,:,mat_index));
            tmp = control_totsh_ttot24_mat(:,:,mat_index);
            tmp(nan_index) = NaN;
            control_totsh_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_prmc_tavg24_mat(:,:,mat_index) = nanmean(exp_prmc_hour_mat,3);
            exp_KEint_tavg24_mat(:,:,mat_index) = nanmean(exp_KEint_hour_mat,3);
            exp_QVint_tavg24_mat(:,:,mat_index) = nanmean(exp_QVint_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_totlh_ttot24_mat(:,:,mat_index) = sum(tmp_exp_totlh_mat,3,'omitnan');
            nan_index = (sum(isnan(tmp_exp_totlh_mat))==size(tmp_exp_totlh_mat,3));
            tmp = exp_totlh_ttot24_mat(:,:,mat_index);
            tmp(nan_index) = NaN;
            exp_totlh_ttot24_mat(:,:,mat_index) = tmp;

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

            % 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_prmc_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_totlh_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            control_totsh_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            control_KEint_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            control_QVint_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_prmc_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_totlh_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            exp_totsh_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            exp_KEint_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);
            exp_QVint_hour_mat = NaN*ones(sum(lon_index),sum(lat_index),24);

        end %d
    end %m
end %y

disp('done reading data')

% 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_prmc_tavg24_mat(:,:,(mat_index+1):end) = [];
exp_prmc_tavg24_mat(:,:,(mat_index+1):end) = [];
control_totevptrns_ttot24_mat(:,:,(mat_index+1):end) = [];
exp_totevptrns_ttot24_mat(:,:,(mat_index+1):end) = [];
control_KEint_tavg24_mat(:,:,(mat_index+1):end) = []; 
exp_KEint_tavg24_mat(:,:,(mat_index+1):end) = []; 
control_totlh_ttot24_mat(:,:,(mat_index+1):end) = [];
exp_totlh_ttot24_mat(:,:,(mat_index+1):end) = [];
control_totsh_ttot24_mat(:,:,(mat_index+1):end) = [];
exp_totsh_ttot24_mat(:,:,(mat_index+1):end) = [];
control_QVint_tavg24_mat(:,:,(mat_index+1):end) = [];
exp_QVint_tavg24_mat(:,:,(mat_index+1):end) = [];
control_land_mask_heat_map_daily(:,:,(mat_index+1):end) = [];
exp_land_mask_heat_map_daily(:,:,(mat_index+1):end) = [];
control_land_mask_daily(:,:,(mat_index+1):end) = [];
exp_land_mask_daily(:,:,(mat_index+1):end) = [];
control_land_mask_time_daily(:,:,(mat_index+1):end) = [];
exp_land_mask_time_daily(:,:,(mat_index+1):end) = [];


%control_totlh_ttot24_mat_corr = NaN*ones(size(control_totlh_ttot24_mat)); 
%exp_totlh_ttot24_mat_corr = NaN*ones(size(exp_totlh_ttot24_mat));
%control_totsh_ttot24_mat_corr = NaN*ones(size(control_totsh_ttot24_mat));
%exp_totsh_ttot24_mat_corr = NaN*ones(size(exp_totsh_ttot24_mat));
%control_totevptrns_ttot24_mat_corr = NaN*ones(size(control_totevptrns_ttot24_mat));
%exp_totevptrns_ttot24_mat_corr = NaN*ones(size(exp_totevptrns_ttot24_mat));
%control_sfmc_tavg24_mat_corr = NaN*ones(size(control_sfmc_tavg24_mat));
%exp_sfmc_tavg24_mat_corr = NaN*ones(size(exp_sfmc_tavg24_mat));
%control_rzmc_tavg24_mat_corr = NaN*ones(size(control_rzmc_tavg24_mat));
%exp_rzmc_tavg24_mat_corr = NaN*ones(size(exp_rzmc_tavg24_mat));
%control_prmc_tavg24_mat_corr = NaN*ones(size(control_prmc_tavg24_mat));
%exp_prmc_tavg24_mat_corr = NaN*ones(size(exp_prmc_tavg24_mat));
%
%% compute matrices mixing correct time steps from land maks
%
%for i = 1:length(datenum_vec)
%
%    if (i>3)
%    lm_tmp = control_land_mask_heat_map_daily(:,:,i);
%    lm_tmp_exp = exp_land_mask_heat_map_daily(:,:,i);
%    sum_control_heat_mask(i) = sum(lm_tmp(:),'omitnan');
%
%    control_lh_tmp = zeros(size(control_totlh_ttot24_mat(:,:,i)));
%    control_sh_tmp = zeros(size(control_totsh_ttot24_mat(:,:,i)));
%    control_totevptrns_tmp = zeros(size(control_totevptrns_ttot24_mat(:,:,i)));
%    control_sfmc_tmp = zeros(size(control_sfmc_tavg24_mat(:,:,i)));
%    control_rzmc_tmp = zeros(size(control_rzmc_tavg24_mat(:,:,i)));
%    control_prmc_tmp = zeros(size(control_prmc_tavg24_mat(:,:,i)));
%
%    exp_lh_tmp = zeros(size(exp_totlh_ttot24_mat(:,:,i)));
%    exp_sh_tmp = zeros(size(exp_totsh_ttot24_mat(:,:,i)));
%    exp_totevptrns_tmp = zeros(size(exp_totevptrns_ttot24_mat(:,:,i)));
%    exp_sfmc_tmp = zeros(size(exp_sfmc_tavg24_mat(:,:,i)));
%    exp_rzmc_tmp = zeros(size(exp_rzmc_tavg24_mat(:,:,i)));
%    exp_prmc_tmp = zeros(size(exp_prmc_tavg24_mat(:,:,i)));
%
%    control_lh_count = zeros(size(control_totlh_ttot24_mat(:,:,i)));
%    control_sh_count = zeros(size(control_totsh_ttot24_mat(:,:,i)));
%    control_totevptrns_count = zeros(size(control_totevptrns_ttot24_mat(:,:,i)));
%    control_sfmc_count = zeros(size(control_sfmc_tavg24_mat(:,:,i)));
%    control_rzmc_count = zeros(size(control_rzmc_tavg24_mat(:,:,i)));
%    control_prmc_count = zeros(size(control_prmc_tavg24_mat(:,:,i)));
%
%    exp_lh_count = zeros(size(exp_totlh_ttot24_mat(:,:,i)));
%    exp_sh_count = zeros(size(exp_totsh_ttot24_mat(:,:,i)));
%    exp_totevptrns_count = zeros(size(exp_totevptrns_ttot24_mat(:,:,i)));
%    exp_sfmc_count = zeros(size(exp_sfmc_tavg24_mat(:,:,i)));
%    exp_rzmc_count = zeros(size(exp_rzmc_tavg24_mat(:,:,i)));
%    exp_prmc_count = zeros(size(exp_prmc_tavg24_mat(:,:,i)));
%
%    for j = 0:3
%
%        control_data_index = (control_land_mask_time_daily(:,:,i)==j);
%        control_data_index = control_data_index;
%
%        tmp = control_totlh_ttot24_mat(:,:,i-j);
%        control_lh_tmp(control_data_index) = control_lh_tmp(control_data_index) + tmp(control_data_index); 
%        control_lh_count(control_data_index) = control_lh_count(control_data_index) + 1;
%       % control_lh_tmp(control_data_index) = control_lh_tmp(control_data_index).*lm_tmp(control_index);
%
%        tmp = control_totsh_ttot24_mat(:,:,i-j);
%        control_sh_tmp(control_data_index) = control_sh_tmp(control_data_index) + tmp(control_data_index); 
%        control_sh_count(control_data_index) = control_sh_count(control_data_index) + 1;
%        %control_sh_tmp(control_data_index) = control_sh_tmp(control_data_index).*lm_tmp(control_index);
%
%        tmp = control_totevptrns_ttot24_mat(:,:,i-j);
%        control_totevptrns_tmp(control_data_index) = control_totevptrns_tmp(control_data_index) + tmp(control_data_index);
%        control_totevptrns_count(control_data_index) = control_totevptrns_count(control_data_index) + 1;
%        %control_totevptrns_tmp(control_data_index) = control_totevptrns_tmp(control_data_index).*lm_tmp(control_index);
%
%        tmp = control_sfmc_tavg24_mat(:,:,i-j);
%        control_sfmc_tmp(control_data_index) = control_sfmc_tmp(control_data_index) + tmp(control_data_index);
%        control_sfmc_count(control_data_index) = control_sfmc_count(control_data_index) + 1;
%        %control_sfmc_tmp(control_data_index) = control_sfmc_tmp(control_data_index).*lm_tmp(control_index);
%
%        tmp = control_rzmc_tavg24_mat(:,:,i-j);
%        control_rzmc_tmp(control_data_index) = control_rzmc_tmp(control_data_index) + tmp(control_data_index);
%        control_rzmc_count(control_data_index) = control_rzmc_count(control_data_index) + 1;
%        %control_rzmc_tmp(control_data_index) = control_rzmc_tmp(control_data_index).*lm_tmp(control_index);
%
%        tmp = control_prmc_tavg24_mat(:,:,i-j);
%        control_prmc_tmp(control_data_index) = control_prmc_tmp(control_data_index) + tmp(control_data_index);
%        control_prmc_count(control_data_index) = control_prmc_count(control_data_index) + 1;
%        %control_rzmc_tmp(control_data_index) = control_rzmc_tmp(control_data_index).*lm_tmp(control_index);
%
%        exp_data_index = (exp_land_mask_time_daily(:,:,i)==j);
%     %   exp_data_index = exp_data_index == 2;
%
%        tmp = exp_totlh_ttot24_mat(:,:,i-j);
%        exp_lh_tmp(exp_data_index) = exp_lh_tmp(exp_data_index) + tmp(exp_data_index);
%        exp_lh_count(exp_data_index) = exp_lh_count(exp_data_index) + 1;
%        %exp_lh_tmp(exp_data_index) = exp_lh_tmp(exp_data_index).*lm_tmp_exp(exp_index);
%
%        tmp = exp_totsh_ttot24_mat(:,:,i-j);
%        exp_sh_tmp(exp_data_index) = exp_sh_tmp(exp_data_index) + tmp(exp_data_index);
%        exp_sh_count(exp_data_index) = exp_sh_count(exp_data_index) + 1;
%        %exp_sh_tmp(exp_data_index) = exp_sh_tmp(exp_data_index).*lm_tmp_exp(exp_index);
%
%        tmp = exp_totevptrns_ttot24_mat(:,:,i-j);
%        exp_totevptrns_tmp(exp_data_index) = exp_totevptrns_tmp(exp_data_index) + tmp(exp_data_index);
%        exp_totevptrns_count(exp_data_index) = exp_totevptrns_count(exp_data_index) + 1;
%        %exp_totevptrns_tmp(exp_data_index) = exp_totevptrns_tmp(exp_data_index).*lm_tmp_exp(exp_index);
%
%        tmp = exp_sfmc_tavg24_mat(:,:,i-j);
%        exp_sfmc_tmp(exp_data_index) = exp_sfmc_tmp(exp_data_index) + tmp(exp_data_index);
%        exp_sfmc_count(exp_data_index) = exp_sfmc_count(exp_data_index) + 1;
%        % exp_sfmc_tmp(exp_data_index) = exp_sfmc_tmp(exp_data_index).*lm_tmp_exp(exp_index);
%
%        tmp = exp_rzmc_tavg24_mat(:,:,i-j);
%        exp_rzmc_tmp(exp_data_index) = exp_rzmc_tmp(exp_data_index) + tmp(exp_data_index);
%        exp_rzmc_count(exp_data_index) = exp_rzmc_count(exp_data_index) + 1;
%        %exp_rzmc_tmp(exp_data_index) = exp_rzmc_tmp(exp_data_index).*lm_tmp_exp(exp_index);
%
%        tmp = exp_prmc_tavg24_mat(:,:,i-j);
%        exp_prmc_tmp(exp_data_index) = exp_prmc_tmp(exp_data_index) + tmp(exp_data_index);
%        exp_prmc_count(exp_data_index) = exp_prmc_count(exp_data_index) + 1;
%        %exp_rzmc_tmp(exp_data_index) = exp_rzmc_tmp(exp_data_index).*lm_tmp_exp(exp_index);
%
%    end % j 
%
%    control_lh_tmp(control_lh_tmp==0) = NaN;
%    control_sh_tmp(control_sh_tmp==0) = NaN;
%    control_totevptrns_tmp(control_totevptrns_tmp==0) = NaN;
%    control_sfmc_tmp(control_sfmc_tmp==0) = NaN;
%    control_rzmc_tmp(control_rzmc_tmp==0) = NaN;
%    control_prmc_tmp(control_prmc_tmp==0) = NaN;
%
%    exp_lh_tmp(exp_lh_tmp==0) = NaN;
%    exp_sh_tmp(exp_sh_tmp==0) = NaN;
%    exp_totevptrns_tmp(exp_totevptrns_tmp==0) = NaN;
%    exp_sfmc_tmp(exp_sfmc_tmp==0) = NaN;
%    exp_rzmc_tmp(exp_rzmc_tmp==0) = NaN;
%    exp_prmc_tmp(exp_prmc_tmp==0) = NaN;
%
%    control_lh_count(control_lh_count==0) = NaN;
%    control_sh_count(control_sh_count==0) = NaN;
%    control_totevptrns_count(control_totevptrns_count==0) = NaN;
%    control_sfmc_count(control_sfmc_count==0) = NaN;
%    control_rzmc_count(control_rzmc_count==0) = NaN;
%    control_prmc_count(control_prmc_count==0) = NaN;
%
%    exp_lh_count(exp_lh_count==0) = NaN;
%    exp_sh_count(exp_sh_count==0) = NaN;
%    exp_totevptrns_count(exp_totevptrns_count==0) = NaN;
%    exp_sfmc_count(exp_sfmc_count==0) = NaN;
%    exp_rzmc_count(exp_rzmc_count==0) = NaN;
%    exp_prmc_count(exp_prmc_count==0) = NaN;
%
%    control_lh_tmp = control_lh_tmp./control_lh_count;
%    control_sh_tmp = control_sh_tmp./control_sh_count;
%    control_totevptrns_tmp = control_totevptrns_tmp./control_totevptrns_count;
%    control_sfmc_tmp = control_sfmc_tmp./control_sfmc_count;
%    control_rzmc_tmp = control_rzmc_tmp./control_rzmc_count;
%    control_prmc_tmp = control_prmc_tmp./control_prmc_count;
%
%    exp_lh_tmp = exp_lh_tmp./exp_lh_count;
%    exp_sh_tmp = exp_sh_tmp./exp_sh_count;
%    exp_totevptrns_tmp = exp_totevptrns_tmp./exp_totevptrns_count;
%    exp_sfmc_tmp = exp_sfmc_tmp./exp_sfmc_count;
%    exp_rzmc_tmp = exp_rzmc_tmp./exp_rzmc_count;
%    exp_prmc_tmp = exp_prmc_tmp./exp_prmc_count;
%
%    control_totlh_ttot24_mat_corr(:,:,i) = control_lh_tmp;
%    exp_totlh_ttot24_mat_corr(:,:,i) = exp_lh_tmp;
%    control_totsh_ttot24_mat_corr(:,:,i) = control_sh_tmp;
%    exp_totsh_ttot24_mat_corr(:,:,i) = exp_sh_tmp;
%    control_totevptrns_ttot24_mat_corr(:,:,i) = control_totevptrns_tmp;
%    exp_totevptrns_ttot24_mat_corr(:,:,i) = exp_totevptrns_tmp;
%    control_sfmc_tavg24_mat_corr(:,:,i) = control_sfmc_tmp;
%    exp_sfmc_tavg24_mat_corr(:,:,i) = exp_sfmc_tmp;
%    control_rzmc_tavg24_mat_corr(:,:,i) = control_rzmc_tmp;
%    exp_rzmc_tavg24_mat_corr(:,:,i) = exp_rzmc_tmp;
%    control_prmc_tavg24_mat_corr(:,:,i) = control_prmc_tmp;
%    exp_prmc_tavg24_mat_corr(:,:,i) = exp_prmc_tmp;
%
%    end % if 
%
%
%end 
%
%
%% compute differences between experiment and control
%
%diff_sfmc_tavg24_mat = exp_sfmc_tavg24_mat_corr - control_sfmc_tavg24_mat_corr;
%diff_rzmc_tavg24_mat = exp_rzmc_tavg24_mat_corr - control_rzmc_tavg24_mat_corr;
%diff_prmc_tavg24_mat = exp_prmc_tavg24_mat_corr - control_prmc_tavg24_mat_corr;
%diff_totevptrns_ttot24_mat = exp_totevptrns_ttot24_mat_corr - control_totevptrns_ttot24_mat_corr;
%diff_KEint_tavg24_mat = exp_KEint_tavg24_mat - control_KEint_tavg24_mat;
%diff_QVint_tavg24_mat = exp_QVint_tavg24_mat - control_QVint_tavg24_mat;

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_prmc_tavg24_mat = exp_prmc_tavg24_mat - control_prmc_tavg24_mat;
diff_totevptrns_ttot24_mat = exp_totevptrns_ttot24_mat - control_totevptrns_ttot24_mat;
diff_KEint_tavg24_mat = exp_KEint_tavg24_mat - control_KEint_tavg24_mat;
diff_QVint_tavg24_mat = exp_QVint_tavg24_mat - control_QVint_tavg24_mat;

diff_totevptrns_ttot24_mat(isnan(diff_sfmc_tavg24_mat)) = NaN;

% 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

cmap_mat_KEint = cmocean('balance',24);
cmap_mat_KEint = flipud(cmap_mat_KEint);
cmap_mat_KEint(length(cmap_mat_KEint)/2:length(cmap_mat_KEint)/2+1,:) = ones(2,3);
s1 = quantile(diff_KEint_tavg24_mat(:),0.01);
s2 = quantile(diff_KEint_tavg24_mat(:),0.99);
caxis_min_KEint = -1*max(abs(s1),abs(s2));
caxis_max_KEint = max(abs(s1),abs(s2));
%caxis_min_KEint = -15000;
%caxis_max_KEint = 15000;


cmap_mat_QVint = cmocean('balance',24);
cmap_mat_QVint = flipud(cmap_mat_QVint);
cmap_mat_QVint(length(cmap_mat_QVint)/2:length(cmap_mat_QVint)/2+1,:) = ones(2,3);
s1 = quantile(diff_QVint_tavg24_mat(:),0.01);
s2 = quantile(diff_QVint_tavg24_mat(:),0.99);
caxis_min_QVint = -1*max(abs(s1),abs(s2));
caxis_max_QVint = max(abs(s1),abs(s2));


cmap_mat_KEint_tot = cmocean('matter',12);
%cmap_mat_KEint_tot = flipud(cmap_mat_KEint_tot);
caxis_min_KEint_tot = min(exp_KEint_tavg24_mat(:));
caxis_max_KEint_tot = max(exp_KEint_tavg24_mat(:));
caxis_min_KEint_tot = 10000;
caxis_max_KEint_tot = 60000;


cmap_mat_QVint_tot = cmocean('ice',12);
cmap_mat_QVint_tot = flipud(cmap_mat_QVint_tot);
caxis_min_QVint_tot = min(exp_QVint_tavg24_mat(:));
caxis_max_QVint_tot = max(exp_QVint_tavg24_mat(:));
%caxis_min_KEint_tot = 10000;
%caxis_max_KEint_tot = 60000;


tmp_mat = [control_sfmc_tavg24_mat exp_sfmc_tavg24_mat];

s1 = quantile(tmp_mat(:),0.99);
s2 = quantile(tmp_mat(:),0.01);
caxis_min_sfmc_raw = s2;
caxis_min_sfmc_raw = round(caxis_min_sfmc_raw,2);
caxis_max_sfmc_raw = s1;
caxis_max_sfmc_raw = round(caxis_max_sfmc_raw,2);

caxis_min_sfmc_raw = 0.05;
caxis_max_sfmc_raw = 0.45;

c_length = caxis_max_sfmc_raw/0.05;
c_length = c_length*2;

cmap_mat_sfmc_raw = cmocean('rain',c_length);
%cmap_mat_sfmc_raw = 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_raw = cmap_mat_sfmc_raw;

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.01;
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_prmc_tavg24_mat;

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

c_length = caxis_max_prmc/0.01;
c_length = c_length*2;

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

tmp_mat = diff_totevptrns_ttot24_mat;

s1 = quantile(tmp_mat(:),0.99);
s2 = quantile(tmp_mat(:),0.01);
s1 = max(tmp_mat(:));
s2 = min(tmp_mat(:));
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.5;
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 = [control_totevptrns_ttot24_mat exp_totevptrns_ttot24_mat];

s1 = quantile(tmp_mat(:),0.99);
s2 = quantile(tmp_mat(:),0.01);
s1 = max(tmp_mat(:));
s2 = min(tmp_mat(:));
caxis_min_totevptrns_raw = s2;
caxis_min_totevptrns_raw = round(caxis_min_totevptrns_raw,1);
caxis_max_totevptrns_raw = s1;
caxis_max_totevptrns_raw = round(caxis_max_totevptrns_raw,1);

c_length = caxis_max_totevptrns_raw/1;
c_length = c_length*2;

cmap_mat_totevptrns_raw = cmocean('rain',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_raw = cmap_mat_totevptrns_raw;

mean_KEint_control = NaN*ones(length(datenum_vec),1);
mean_QVint_control = NaN*ones(length(datenum_vec),1);
mean_totlh_control = NaN*ones(length(datenum_vec),1);
mean_totsh_control = NaN*ones(length(datenum_vec),1);
mean_totevptrns_control = NaN*ones(length(datenum_vec),1);
mean_sfmc_control = NaN*ones(length(datenum_vec),1);
mean_rzmc_control = NaN*ones(length(datenum_vec),1);

mean_KEint_exp = NaN*ones(length(datenum_vec),1);
mean_QVint_exp = NaN*ones(length(datenum_vec),1);
mean_totlh_exp = NaN*ones(length(datenum_vec),1);
mean_totsh_exp = NaN*ones(length(datenum_vec),1);
mean_totevptrns_exp = NaN*ones(length(datenum_vec),1);
mean_sfmc_exp = NaN*ones(length(datenum_vec),1);
mean_rzmc_exp = NaN*ones(length(datenum_vec),1);

cum_KEint_control = zeros(length(datenum_vec),1);
cum_QVint_control = zeros(length(datenum_vec),1);
cum_totlh_control = zeros(length(datenum_vec),1);
cum_totsh_control = zeros(length(datenum_vec),1);
cum_totevptrns_control = zeros(length(datenum_vec),1);
cum_sfmc_control = zeros(length(datenum_vec),1);
cum_rzmc_control = zeros(length(datenum_vec),1);

cum_KEint_exp = zeros(length(datenum_vec),1);
cum_QVint_exp = zeros(length(datenum_vec),1);
cum_totlh_exp = zeros(length(datenum_vec),1);
cum_totsh_exp = zeros(length(datenum_vec),1);
cum_totevptrns_exp = zeros(length(datenum_vec),1);
cum_sfmc_exp = zeros(length(datenum_vec),1);
cum_rzmc_exp = zeros(length(datenum_vec),1);

for i = 1:length(datenum_vec)
   
    disp('plotting')

%    cmap_mat_KEint = cmocean('balance',24);
%    cmap_mat_KEint = flipud(cmap_mat_KEint);
%    cmap_mat_KEint(length(cmap_mat_KEint)/2:length(cmap_mat_KEint)/2+1,:) = ones(2,3);
%    tmp = diff_KEint_tavg24_mat(:,:,i);
%    s1 = min(tmp(:));
%    s2 = max(tmp(:));
%    caxis_min_KEint = -1*max(abs(s1),abs(s2));
%    caxis_max_KEint = max(abs(s1),abs(s2));
%
%    cmap_mat_KEint_tot = cmocean('matter',10);
%    cmap_mat_KEint_tot = flipud(cmap_mat_KEint_tot);
%    tmp = exp_KEint_tavg24_mat(:,:,i);
%    caxis_min_KEint_tot = min(tmp(:));
%    caxis_max_KEint_tot = max(tmp(:));


    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','off',...
         '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),diff_KEint_tavg24_mat(:,:,i)')
    caxis([caxis_min_KEint caxis_max_KEint]);
    set(gcf,'Colormap',cmap_mat_KEint);
    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 vert. int. KE [kg m^{-1} s^{-2}] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    h = colorbar;
    set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,[save_path storm_name '_ana_KEint_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(lat_index),lon_vec(lon_index),control_KEint_tavg24_mat(:,:,i)')
    caxis([caxis_min_KEint_tot caxis_max_KEint_tot]);
    set(gcf,'Colormap',cmap_mat_KEint_tot);
    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(['vert. int. KE [kg m^{-1} s^{-2}] CNTRL ' 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 '_ana_KEint_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(lat_index),lon_vec(lon_index),exp_KEint_tavg24_mat(:,:,i)')
    caxis([caxis_min_KEint_tot caxis_max_KEint_tot]);
    set(gcf,'Colormap',cmap_mat_KEint_tot);
    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(['vert. int. KE [kg m^{-1} s^{-2}] EXP ' 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 '_ana_KEint_EXP_' 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','off',...
         '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),diff_QVint_tavg24_mat(:,:,i)')
    caxis([caxis_min_QVint caxis_max_QVint]);
    set(gcf,'Colormap',cmap_mat_QVint);
    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 vert. int. QV [kg/kg] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    h = colorbar;
    set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,[save_path storm_name '_ana_QVint_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(lat_index),lon_vec(lon_index),control_QVint_tavg24_mat(:,:,i)')
    caxis([caxis_min_QVint_tot caxis_max_QVint_tot]);
    set(gcf,'Colormap',cmap_mat_QVint_tot);
    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(['vert. int. QV [kg/kg] CNTRL ' 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 '_ana_QVint_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(lat_index),lon_vec(lon_index),exp_QVint_tavg24_mat(:,:,i)')
    caxis([caxis_min_QVint_tot caxis_max_QVint_tot]);
    set(gcf,'Colormap',cmap_mat_QVint_tot);
    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(['vert. int. QV [kg/kg] EXP ' 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 '_ana_QVint_EXP_' 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,...
         'ParallelLabel','on','PLabelLocation',5,...
          'FontSize',12)
%    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),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 (SMAP\_ANA - CTRL\_ANA) [m^{3} m^{-3}] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    h = colorbar;
    set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,[save_path storm_name '_diag_SFMC_diff_landMask_interp_' bt_time_str '_'  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,...
         'ParallelLabel','on','PLabelLocation',5,...
          'FontSize',12)
%    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),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\_ANA [m^{3} m^{-3}] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    h = colorbar;
    set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,[save_path storm_name '_diag_SFMC_CNTRL_landMask_interp_' bt_time_str '_'  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,...
         'ParallelLabel','on','PLabelLocation',5,...
          'FontSize',12)
%    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),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 SMAP\_ANA [m^{3} m^{-3}] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    h = colorbar;
    set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,[save_path storm_name '_diag_SFMC_EXP_landMask_interp_' bt_time_str '_'  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)
    setm(ax,'MLineLocation',5,'PLineLocation',5,...
         'ParallelLabel','on','PLabelLocation',5,...
          'FontSize',12)
    geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
    pcolorm(lat_vec(lat_index),lon_vec(lon_index),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 (SMAP\_ANA - CTRL\_ANA) [m^{3} m^{-3}] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    h = colorbar;
    set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,[save_path storm_name '_diag_RZMC_diff_landMask_interp_' bt_time_str '_' 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)
    setm(ax,'MLineLocation',5,'PLineLocation',5,...
         'ParallelLabel','on','PLabelLocation',5,...
          'FontSize',12)
    geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
    pcolorm(lat_vec(lat_index),lon_vec(lon_index),diff_prmc_tavg24_mat(:,:,i)')
    caxis([caxis_min_prmc caxis_max_prmc]);
    set(gcf,'Colormap',cmap_prmc);
    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 PRMC (SMAP\_ANA - CTRL\_ANA) [m^{3} m^{-3}] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    h = colorbar;
    set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,[save_path storm_name '_diag_PRMC_diff_landMask_interp_' bt_time_str '_' 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)
    setm(ax,'MLineLocation',5,'PLineLocation',5,...
         'ParallelLabel','on','PLabelLocation',5,...
          'FontSize',12)
    geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
    pcolorm(lat_vec(lat_index),lon_vec(lon_index),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 (SMAP\_ANA - CTRL\_ANA) [mm] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    h = colorbar;
    set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,[save_path storm_name '_diag_TOTEVPTRNS_diff_landMask_interp_' bt_time_str '_' 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)
    setm(ax,'MLineLocation',5,'PLineLocation',5,...
         'ParallelLabel','on','PLabelLocation',5,...
          'FontSize',12)
    geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
    pcolorm(lat_vec(lat_index),lon_vec(lon_index),control_totevptrns_ttot24_mat(:,:,i)')
    caxis([caxis_min_totevptrns_raw caxis_max_totevptrns_raw]);
    set(gcf,'Colormap',cmap_totevptrns_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(['total daily evapotranspiration CTRL\_ANA [mm] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    h = colorbar;
    set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,[save_path storm_name '_diag_TOTEVPTRNS_CNTRL_landMask_interp_' bt_time_str '_' 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)
    setm(ax,'MLineLocation',5,'PLineLocation',5,...
         'ParallelLabel','on','PLabelLocation',5,...
          'FontSize',12)
    geoshow('landareas.shp', 'FaceColor', [0.7 0.7 0.7])
    pcolorm(lat_vec(lat_index),lon_vec(lon_index),exp_totevptrns_ttot24_mat(:,:,i)')
    caxis([caxis_min_totevptrns_raw caxis_max_totevptrns_raw]);
    set(gcf,'Colormap',cmap_totevptrns_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(['total daily evapotranspiration SMAP\_ANA [mm] ' datestr(datenum_vec(i),'yyyy-mm-dd')])
    h = colorbar;
    set(h,'Location','southoutside','Orientation','horizontal','FontSize',14);
    set(gcf,'paperunits','centimeters')
    set(gcf,'paperposition',[1 1 17 15])
    saveas(gcf,[save_path storm_name '_diag_TOTEVPTRNS_EXP_landMask_interp_' bt_time_str '_' datestr(datenum_vec(i),'yyyy-mm-dd')],'png')

    if (i>3)
    lm_tmp = control_land_mask_heat_map_daily(:,:,i);
    sum_control_heat_mask(i) = sum(lm_tmp(:),'omitnan'); 

    tmp = control_KEint_tavg24_mat(:,:,i);
    mean_KEint_control(i) = nanmean(tmp(:));
    tmp = control_QVint_tavg24_mat(:,:,i);
    mean_QVint_control(i) = nanmean(tmp(:));

    control_lh_tmp = zeros(size(control_totlh_ttot24_mat(:,:,i)));
    control_sh_tmp = zeros(size(control_totsh_ttot24_mat(:,:,i)));
    control_totevptrns_tmp = zeros(size(control_totevptrns_ttot24_mat(:,:,i)));
    control_sfmc_tmp = zeros(size(control_sfmc_tavg24_mat(:,:,i)));
    control_rzmc_tmp = zeros(size(control_rzmc_tavg24_mat(:,:,i)));

    for j = 1:3

        control_data_index = control_land_mask_daily(:,:,i) + (control_land_mask_time_daily(:,:,i)==j);
        control_data_index = control_data_index == 2;

%        tmp = control_totlh_ttot24_mat(:,:,i-j);
%        control_lh_tmp(control_data_index) = control_lh_tmp(control_data_index) + tmp(control_data_index);
%        control_lh_tmp(control_data_index) = control_lh_tmp(control_data_index).*lm_tmp(control_data_index);
%
%        tmp = control_totsh_ttot24_mat(:,:,i-j);
%        control_sh_tmp(control_data_index) = control_sh_tmp(control_data_index) + tmp(control_data_index);
%        control_sh_tmp(control_data_index) = control_sh_tmp(control_data_index).*lm_tmp(control_data_index);
%
%        tmp = control_totevptrns_ttot24_mat(:,:,i-j);
%        control_totevptrns_tmp(control_data_index) = control_totevptrns_tmp(control_data_index) + tmp(control_data_index);
%        control_totevptrns_tmp(control_data_index) = control_totevptrns_tmp(control_data_index).*lm_tmp(control_data_index);
%
%        tmp = control_sfmc_tavg24_mat(:,:,i-j);
%        control_sfmc_tmp(control_data_index) = control_sfmc_tmp(control_data_index) + tmp(control_data_index);
%        control_sfmc_tmp(control_data_index) = control_sfmc_tmp(control_data_index).*lm_tmp(control_data_index);
%
%        tmp = control_rzmc_tavg24_mat(:,:,i-j);
%        control_rzmc_tmp(control_data_index) = control_rzmc_tmp(control_data_index) + tmp(control_data_index);
%        control_rzmc_tmp(control_data_index) = control_rzmc_tmp(control_data_index).*lm_tmp(control_data_index);

        tmp = control_totlh_ttot24_mat(:,:,i-j);
        control_lh_tmp = control_lh_tmp + tmp.*lm_tmp;
        control_lh_tmp = control_lh_tmp;

        tmp = control_totsh_ttot24_mat(:,:,i-j);
        control_sh_tmp = control_sh_tmp + tmp.*lm_tmp;
        control_sh_tmp = control_sh_tmp;

        tmp = control_totevptrns_ttot24_mat(:,:,i-j);
        control_totevptrns_tmp = control_totevptrns_tmp + tmp.*lm_tmp;
        control_totevptrns_tmp = control_totevptrns_tmp;

        tmp = control_sfmc_tavg24_mat(:,:,i-j);
        control_sfmc_tmp = control_sfmc_tmp + tmp.*lm_tmp;
        control_sfmc_tmp = control_sfmc_tmp;

        tmp = control_rzmc_tavg24_mat(:,:,i-j);
        control_rzmc_tmp = control_rzmc_tmp + tmp.*lm_tmp;
        control_rzmc_tmp = control_rzmc_tmp;

    end 

    control_lh_tmp(control_lh_tmp==0) = NaN;
    control_sh_tmp(control_sh_tmp==0) = NaN;
    control_totevptrns_tmp(control_totevptrns_tmp==0) = NaN;
    control_sfmc_tmp(control_sfmc_tmp==0) = NaN;
    control_rzmc_tmp(control_rzmc_tmp==0) = NaN;

    mean_totlh_control(i) = sum(control_lh_tmp(:),'omitnan')./sum(3*lm_tmp(:),'omitnan');
    mean_totsh_control(i) = sum(control_sh_tmp(:),'omitnan')./sum(3*lm_tmp(:),'omitnan');
    mean_totevptrns_control(i) = sum(control_totevptrns_tmp(:),'omitnan')./sum(3*lm_tmp(:),'omitnan');
    mean_sfmc_control(i) = sum(control_sfmc_tmp(:),'omitnan')./sum(3*lm_tmp(:),'omitnan');
    mean_rzmc_control(i) = sum(control_rzmc_tmp(:),'omitnan')./sum(3*lm_tmp(:),'omitnan');

    lm_tmp = exp_land_mask_heat_map_daily(:,:,i);
    sum_exp_heat_mask(i) = sum(lm_tmp(:),'omitnan');

    tmp = exp_KEint_tavg24_mat(:,:,i);
    mean_KEint_exp(i) = nanmean(tmp(:));
    tmp = exp_QVint_tavg24_mat(:,:,i);
    mean_QVint_exp(i) = nanmean(tmp(:));

    exp_lh_tmp = zeros(size(exp_totlh_ttot24_mat(:,:,i)));
    exp_sh_tmp = zeros(size(exp_totsh_ttot24_mat(:,:,i)));
    exp_totevptrns_tmp = zeros(size(exp_totevptrns_ttot24_mat(:,:,i)));
    exp_sfmc_tmp = zeros(size(exp_sfmc_tavg24_mat(:,:,i)));
    exp_rzmc_tmp = zeros(size(exp_rzmc_tavg24_mat(:,:,i)));

    for j = 1:3

        exp_data_index = exp_land_mask_daily(:,:,i) + (exp_land_mask_time_daily(:,:,i)==j);
        exp_data_index = exp_data_index == 2;

%        tmp = exp_totlh_ttot24_mat(:,:,i-j);
%        exp_lh_tmp(exp_data_index) = exp_lh_tmp(exp_data_index) + tmp(exp_data_index);
%        exp_lh_tmp(exp_data_index) = exp_lh_tmp(exp_data_index).*lm_tmp(exp_data_index);
%
%        tmp = exp_totsh_ttot24_mat(:,:,i-j);
%        exp_sh_tmp(exp_data_index) = exp_sh_tmp(exp_data_index) + tmp(exp_data_index);
%        exp_sh_tmp(exp_data_index) = exp_sh_tmp(exp_data_index).*lm_tmp(exp_data_index);
%
%        tmp = exp_totevptrns_ttot24_mat(:,:,i-j);
%        exp_totevptrns_tmp(exp_data_index) = exp_totevptrns_tmp(exp_data_index) + tmp(exp_data_index);
%        exp_totevptrns_tmp(exp_data_index) = exp_totevptrns_tmp(exp_data_index).*lm_tmp(exp_data_index);
%
%        tmp = exp_sfmc_tavg24_mat(:,:,i-j);
%        exp_sfmc_tmp(exp_data_index) = exp_sfmc_tmp(exp_data_index) + tmp(exp_data_index);
%        exp_sfmc_tmp(exp_data_index) = exp_sfmc_tmp(exp_data_index).*lm_tmp(exp_data_index);
%
%        tmp = exp_rzmc_tavg24_mat(:,:,i-j);
%        exp_rzmc_tmp(exp_data_index) = exp_rzmc_tmp(exp_data_index) + tmp(exp_data_index);
%        exp_rzmc_tmp(exp_data_index) = exp_rzmc_tmp(exp_data_index).*lm_tmp(exp_data_index);

        tmp = exp_totlh_ttot24_mat(:,:,i-j);
        exp_lh_tmp = exp_lh_tmp + tmp.*lm_tmp;
        exp_lh_tmp = exp_lh_tmp;

        tmp = exp_totsh_ttot24_mat(:,:,i-j);
        exp_sh_tmp = exp_sh_tmp + tmp.*lm_tmp;
        exp_sh_tmp = exp_sh_tmp;

        tmp = exp_totevptrns_ttot24_mat(:,:,i-j);
        exp_totevptrns_tmp = exp_totevptrns_tmp + tmp.*lm_tmp;
        exp_totevptrns_tmp = exp_totevptrns_tmp;

        tmp = exp_sfmc_tavg24_mat(:,:,i-j);
        exp_sfmc_tmp = exp_sfmc_tmp + tmp.*lm_tmp;
        exp_sfmc_tmp = exp_sfmc_tmp;

        tmp = exp_rzmc_tavg24_mat(:,:,i-j);
        exp_rzmc_tmp = exp_rzmc_tmp + tmp.*lm_tmp;
        exp_rzmc_tmp = exp_rzmc_tmp;

    end

    exp_lh_tmp(exp_lh_tmp==0) = NaN;
    exp_sh_tmp(exp_sh_tmp==0) = NaN;
    exp_totevptrns_tmp(exp_totevptrns_tmp==0) = NaN;
    exp_sfmc_tmp(exp_sfmc_tmp==0) = NaN;
    exp_rzmc_tmp(exp_rzmc_tmp==0) = NaN;

    mean_totlh_exp(i) = sum(exp_lh_tmp(:),'omitnan')./sum(3*lm_tmp(:),'omitnan');
    mean_totsh_exp(i) = sum(exp_sh_tmp(:),'omitnan')./sum(3*lm_tmp(:),'omitnan');
    mean_totevptrns_exp(i) = sum(exp_totevptrns_tmp(:),'omitnan')./sum(3*lm_tmp(:),'omitnan');
    mean_sfmc_exp(i) = sum(exp_sfmc_tmp(:),'omitnan')./sum(3*lm_tmp(:),'omitnan');
    mean_rzmc_exp(i) = sum(exp_rzmc_tmp(:),'omitnan')./sum(3*lm_tmp(:),'omitnan');

    if (i==1)
     
        cum_KEint_control(i) = mean_KEint_control(i);
        cum_QVint_control(i) = mean_QVint_control(i);
        cum_totlh_control(i) = mean_totlh_control(i);
        cum_totsh_control(i) = mean_totsh_control(i);
        cum_totevptrns_control(i) = mean_totevptrns_control(i);
        cum_sfmc_control(i) = mean_sfmc_control(i);
        cum_rzmc_control(i) = mean_rzmc_control(i);

        cum_KEint_exp(i) = mean_KEint_exp(i);
        cum_QVint_exp(i) = mean_QVint_exp(i);
        cum_totlh_exp(i) = mean_totlh_exp(i);
        cum_totsh_exp(i) = mean_totsh_exp(i);
        cum_totevptrns_exp(i) = mean_totevptrns_exp(i);
        cum_sfmc_exp(i) = mean_sfmc_exp(i);
        cum_rzmc_exp(i) = mean_rzmc_exp(i);

    else

        cum_KEint_control(i) = cum_KEint_control(i-1) + mean_KEint_control(i);
        cum_QVint_control(i) = cum_QVint_control(i-1) + mean_QVint_control(i);
        cum_totlh_control(i) = cum_totlh_control(i-1) + mean_totlh_control(i);
        cum_totsh_control(i) = cum_totsh_control(i-1) + mean_totsh_control(i);
        cum_totevptrns_control(i) = cum_totevptrns_control(i-1) + mean_totevptrns_control(i);
        cum_sfmc_control(i) = cum_sfmc_control(i-1) + mean_sfmc_control(i);
        cum_rzmc_control(i) = cum_rzmc_control(i-1) + mean_rzmc_control(i);
    
        cum_KEint_exp(i) = cum_KEint_exp(i-1) + mean_KEint_exp(i);
        cum_QVint_exp(i) = cum_QVint_exp(i-1) + mean_QVint_exp(i);
        cum_totlh_exp(i) = cum_totlh_exp(i-1) + mean_totlh_exp(i);
        cum_totsh_exp(i) = cum_totsh_exp(i-1) + mean_totsh_exp(i);
        cum_totevptrns_exp(i) = cum_totevptrns_exp(i-1) + mean_totevptrns_exp(i);
        cum_sfmc_exp(i) = cum_sfmc_exp(i-1) + mean_sfmc_exp(i);
        cum_rzmc_exp(i) = cum_rzmc_exp(i-1) + mean_rzmc_exp(i);

    end 
    end 
end 

mean_KEint_plot_mat = [mean_KEint_exp mean_KEint_control];

figure
bar(datenum_vec,mean_KEint_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('mean  vert. int. KE [kg m^{-1} s^{-2}]')
legend('SMAP DA','Control','Location','northwest')
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 '_mean_KEint_bar'],'png')

mean_QVint_plot_mat = [mean_QVint_exp mean_QVint_control];

figure
bar(datenum_vec,mean_QVint_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('mean  vert. int. QV [kg/kg]')
legend('SMAP DA','Control','Location','northwest')
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 '_mean_QVint_bar'],'png')


mean_totlh_plot_mat = [mean_totlh_exp mean_totlh_control];

figure
bar(datenum_vec,mean_totlh_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('mean latent heat [mm]')
legend('SMAP DA','Control','Location','northwest')
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 '_mean_LH_bar_interp'],'png')

mean_totsh_plot_mat = [mean_totsh_exp mean_totsh_control];

figure
bar(datenum_vec,mean_totsh_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('mean sensible heat [mm]')
legend('SMAP DA','Control','Location','northwest')
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 '_mean_SH_bar_interp'],'png')

cum_KEint_plot_mat = [cum_KEint_exp cum_KEint_control];

figure
bar(datenum_vec,cum_KEint_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('cumulative vert. int. KE [kg m^{-1} s^{-2}]')
legend('SMAP DA','Control','Location','northwest')
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 '_cum_KEint_bar'],'png')

cum_QVint_plot_mat = [cum_QVint_exp cum_QVint_control];

figure
bar(datenum_vec,cum_QVint_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('cumulative vert. int. QV [kg/kg]')
legend('SMAP DA','Control','Location','northwest')
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 '_cum_QVint_bar'],'png')


cum_totlh_plot_mat = [cum_totlh_exp cum_totlh_control];

figure
bar(datenum_vec,cum_totlh_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('cumulative latent heat [mm]')
legend('SMAP DA','Control','Location','northwest')
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 '_cum_LH_bar_interp'],'png')

cum_totsh_plot_mat = [cum_totsh_exp cum_totsh_control];

figure
bar(datenum_vec,cum_totsh_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('cumulative sensible heat [mm]')
legend('SMAP DA','Control','Location','northwest')
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 '_cum_SH_bar_interp'],'png')

cum_totevptrns_plot_mat = [cum_totevptrns_exp cum_totevptrns_control];

figure
bar(datenum_vec,cum_totevptrns_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('cumulative evapotranspiration [mm]')
legend('SMAP DA','Control','Location','northwest')
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 '_cum_EVPTRNS_bar_interp'],'png')

sum_heat_mask = [sum_exp_heat_mask' sum_control_heat_mask'];

figure
bar(datenum_vec,sum_heat_mask)
datetick('x','dd-mmm')
xlabel('Day')
title('number of gridcells [-]')
legend('SMAP DA','Control','Location','northeast')
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 '_heat_map_bar_interp'],'png')

sum_heat_mask_area = sum_heat_mask*209.7;

figure
plot(datenum_vec(4:end),sum_heat_mask_area(4:end,2),'--','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on 
plot(datenum_vec(4:end),sum_heat_mask_area(4:end,1),'--k','LineWidth',2)
grid on
datetick('x','dd')
ax= gca;
ax.FontSize = 14;
xlabel('March', 'FontSize',14)
ylabel('Contributing land area [km^2]','FontSize',14)
title('(a) Contributing land area [km^2]','FontSize',14)
legend('CNTRL\_ANA','SMAP\_ANA','Location','northeast')
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 '_heat_map_area_line_interp'],'png')

cum_diff_KEint = cum_KEint_exp - cum_KEint_control;

figure
bar(datenum_vec,cum_diff_KEint)
datetick('x','dd-mmm')
xlabel('Day')
title('\Delta cumulative vert. int. KE [kg m^{-1} s^{-2}]')
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 '_cum_diff_KEint_bar'],'png')

cum_sfmc_plot_mat = [cum_sfmc_exp cum_sfmc_control];

figure
bar(datenum_vec,cum_sfmc_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('cumulative surface soil moisture [m^3 m^{-3}]')
legend('SMAP DA','Control','Location','northwest')
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 '_cum_SFMC_bar_interp'],'png')

cum_rzmc_plot_mat = [cum_rzmc_exp cum_rzmc_control];

figure
bar(datenum_vec,cum_rzmc_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('cumulative root-zone soil moisture [m^3 m^{-3}]')
legend('SMAP DA','Control','Location','northwest')
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 '_cum_RZMC_bar_interp'],'png')

mean_sfmc_plot_mat = [mean_sfmc_exp mean_sfmc_control];

figure
bar(datenum_vec,mean_sfmc_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('mean surface soil moisture [m^3 m^{-3}]')
legend('SMAP DA','Control','Location','northwest')
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 '_mean_SFMC_bar_interp'],'png')

mean_rzmc_plot_mat = [mean_rzmc_exp mean_rzmc_control];

figure
bar(datenum_vec,mean_rzmc_plot_mat)
datetick('x','dd-mmm')
xlabel('Day')
title('mean root-zone soil moisture [m^3 m^{-3}]')
legend('SMAP DA','Control','Location','northwest')
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 '_mean_RZMC_bar_interp'],'png')


figure
colororder({'k','k'})
yyaxis left
plot(datenum_vec(4:end),mean_QVint_plot_mat(4:end,2),'--','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on
plot(datenum_vec(4:end),mean_QVint_plot_mat(4:end,1),'--k','LineWidth',2)
ylabel('mean QV [kg/kg]', 'FontSize',14)
yyaxis right
plot(datenum_vec(4:end),cum_QVint_plot_mat(4:end,2),'-','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on
plot(datenum_vec(4:end),cum_QVint_plot_mat(4:end,1),'-k','LineWidth',2)
ylabel('cumulative QV [kg/kg]', 'FontSize',14)
grid on
datetick('x','dd')
ax= gca;
ax.FontSize = 14;
xlabel('March', 'FontSize',14)
title('(e) vert. int. specific humidity [kg/kg]', 'FontSize',14)
%legend('mean SMAP\_FCST','mean CNTRL\_FCST','cumulative SMAP\_FCST','cumulative CNTRL\_FCST','Location','northwest','FontSize',14)
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 '_mean_cum_QVint_line'],'png')

%figure
%colororder({'k','k'})
%h1 = yyaxis left;
%plot(datenum_vec(4:end),mean_QVint_plot_mat(4:end,2),'--','Color',[0.5 0.5 0.5],'LineWidth',3)
%hold on
%h2 = yyaxis right;
%plot(datenum_vec(4:end),cum_QVint_plot_mat(4:end,2),'-','Color',[0.5 0.5 0.5],'LineWidth',3)
%hold on 
%h1 = yyaxis left;
%plot(datenum_vec(4:end),mean_QVint_plot_mat(4:end,1),'--k','LineWidth',2)
%h2 = yyaxis right;
%plot(datenum_vec(4:end),cum_QVint_plot_mat(4:end,1),'-k','LineWidth',2)
%set(h1,'YLabel','mean QV [kg/kg]', 'FontSize',14)
%set(h2,'YLabel','cumulative QV [kg/kg]', 'FontSize',14)
%grid on
%datetick('x','dd')
%ax= gca;
%ax.FontSize = 14;
%xlabel('March', 'FontSize',14)
%title('(e) vert. int. specific humidity [kg/kg]', 'FontSize',14)
%%legend('mean SMAP\_FCST','mean CNTRL\_FCST','cumulative SMAP\_FCST','cumulative CNTRL\_FCST','Location','northwest','FontSize',14)
%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 '_mean_cum_QVint_line_test'],'png')

figure
colororder({'k','k'})
yyaxis left
plot(datenum_vec(4:end),mean_KEint_plot_mat(4:end,2),'--','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on
plot(datenum_vec(4:end),mean_KEint_plot_mat(4:end,1),'--k','LineWidth',2)
ylabel('mean KE [kg m^{-1} s^{-2}]', 'FontSize',14)
yyaxis right
plot(datenum_vec(4:end),cum_KEint_plot_mat(4:end,2),'-','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on
plot(datenum_vec(4:end),cum_KEint_plot_mat(4:end,1),'-k','LineWidth',2)
ylabel('cumulative  KE [kg m^{-1} s^{-2}]', 'FontSize',14)
grid on
datetick('x','dd')
ax= gca;
ax.FontSize = 14;
xlabel('March', 'FontSize',14)
title('(f) vert. int. kinetic energy [kg m^{-1} s^{-2}]', 'FontSize',14)
%legend('mean SMAP\_FCST','mean CNTRL\_FCST','cumulative SMAP\_FCST','cumulative CNTRL\_FCST','Location','northwest','FontSize',14)
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 '_mean_cum_KEint_line'],'png')

figure 
colororder({'k','k'})
yyaxis left 
plot(datenum_vec(4:end),mean_totlh_plot_mat(4:end,2),'--','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on 
plot(datenum_vec(4:end),mean_totlh_plot_mat(4:end,1),'--k','LineWidth',2)
ylabel('mean daily latent heat [mm]', 'FontSize',14)
yyaxis right
plot(datenum_vec(4:end),cum_totlh_plot_mat(4:end,2),'-','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on
plot(datenum_vec(4:end),cum_totlh_plot_mat(4:end,1),'-k','LineWidth',2)
ylabel('cumulative daily latent heat [mm]', 'FontSize',14)
grid on
datetick('x','dd')
ax= gca;
ax.FontSize = 14;
xlabel('March', 'FontSize',14)
title('(c) latent heat [mm]', 'FontSize',14)
%legend('mean SMAP\_FCST','mean CNTRL\_FCST','cumulative SMAP\_FCST','cumulative CNTRL\_FCST','Location','northwest','FontSize',14)
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 '_mean_cum_LH_line'],'png')

figure 
colororder({'k','k'})
yyaxis left 
plot(datenum_vec(4:end),mean_totsh_plot_mat(4:end,2),'--','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on 
plot(datenum_vec(4:end),mean_totsh_plot_mat(4:end,1),'--k','LineWidth',2)
ylabel('mean daily sensible heat [mm]', 'FontSize',14)
yyaxis right
plot(datenum_vec(4:end),cum_totsh_plot_mat(4:end,2),'-','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on
plot(datenum_vec(4:end),cum_totsh_plot_mat(4:end,1),'-k','LineWidth',2)
ylabel('cumulative daily sensible heat [mm]', 'FontSize',14)
grid on
datetick('x','dd')
ax= gca;
ax.FontSize = 14;
xlabel('March', 'FontSize',14)
title('(d) sensible heat [mm]', 'FontSize',14)
%legend('mean SMAP\_FCST','mean CNTRL\_FCST','cumulative SMAP\_FCST','cumulative CNTRL\_FCST','Location','northwest','FontSize',14)
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 '_mean_cum_SH_line'],'png')

figure
colororder({'k','k'})
yyaxis left
plot(datenum_vec(4:end),mean_sfmc_plot_mat(4:end,2),'--','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on
plot(datenum_vec(4:end),mean_sfmc_plot_mat(4:end,1),'--k','LineWidth',2)
ylabel('mean daily SFMC [m^3 m^{-3}]', 'FontSize',14)
yyaxis right
plot(datenum_vec(4:end),cum_sfmc_plot_mat(4:end,2),'-','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on
plot(datenum_vec(4:end),cum_sfmc_plot_mat(4:end,1),'-k','LineWidth',2)
ylabel('cumulative daily SFMC [m^3 m^{-3}]', 'FontSize',14)
grid on
datetick('x','dd')
ax= gca;
ax.FontSize = 14;
xlabel('March', 'FontSize',14)
title('(b) surface soil moisture [m^3 m^{-3}]', 'FontSize',14)
legend('mean CNTRL\_ANA','mean SMAP\_ANA','cumulative CNTRL\_ANA','cumulative SMAP\_ANA','Location','northwest','FontSize',14)
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 '_mean_cum_SFMC_line'],'png')

figure
colororder({'k','k'})
yyaxis left
plot(datenum_vec(4:end),mean_rzmc_plot_mat(4:end,2),'--','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on
plot(datenum_vec(4:end),mean_rzmc_plot_mat(4:end,1),'--k','LineWidth',2)
ylabel('mean daily RZMC [m^3 m^{-3}]', 'FontSize',14)
yyaxis right
plot(datenum_vec(4:end),cum_rzmc_plot_mat(4:end,2),'-','Color',[0.7 0.7 0.7],'LineWidth',4)
hold on
plot(datenum_vec(4:end),cum_rzmc_plot_mat(4:end,1),'-k','LineWidth',2)
ylabel('cumulative daily RZMC [m^3 m^{-3}]', 'FontSize',14)
grid on
datetick('x','dd')
ax= gca;
ax.FontSize = 14;
xlabel('March', 'FontSize',14)
title('(b) root zone soil moisture [m^3 m^{-3}]', 'FontSize',14)
%legend('mean SMAP\_FCST','mean CNTRL\_FCST','cumulative SMAP\_FCST','cumulative CNTRL\_FCST','Location','northwest','FontSize',14)
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 '_mean_cum_RZMC_line'],'png')
