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';
ecmwf_dir = '/discover/nobackup/projects/gmao/g6dev/ltakacs/ecmwf/';
ecmwf_name = 'ecmwf';

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

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

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

%---------------------------------------------
% 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;
f_index = 0;
f_count = zeros(1,21);

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

            control_dir_name = [control_dir control_name '/prog/Y' y_str '/M' m_str '/D' d_str '/H21/'];
            control_file_list = dir([control_dir_name '*inst1_2d_asm_Nx*']);
            control_3D_list = dir([control_dir_name '*inst3_3d_asm_Np*']);

            exp_dir_name = [exp_dir exp_name '/prog/Y' y_str '/M' m_str '/D' d_str '/H21/'];
            exp_file_list = dir([exp_dir_name '*inst1_2d_asm_Nx*']);
            exp_3D_list = dir([exp_dir_name '*inst3_3d_asm_Np*']);

            if ((length(control_file_list)==21) & (length(exp_file_list)==21)) % only include days where all forecasts are available for both experiments       
            % loop over time steps
            for f = 1:length(control_file_list)

                f_index = f_index + 1;
                f_count(f) = f_count(f) + 1;

                control_file = [control_dir_name control_file_list(f).name];
                exp_file = [exp_dir_name exp_file_list(f).name];

                date_str = control_file(end-17:end-10);
                time_str = control_file(end-8:end-7);

                ver_dir_name = [ecmwf_dir 'Y' date_str(1:4) '/M' date_str(5:6) '/'];
                ver_file = [ver_dir_name 'ecmwf.inst3_3d_wxm_Np.' date_str '_' time_str 'z+' date_str '_' time_str 'z.nc4'];

                % on first day with data, initialize matrices
                if (f_index == 1)

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

                   lon_vec(lon_vec<0) = lon_vec(lon_vec<0) + 360;
                   [lat_mat,lon_mat] = meshgrid(lat_vec,lon_vec);

                   ver_lat_vec = ncread(ver_file,'lat');
                   ver_lon_vec = ncread(ver_file,'lon');
                  % ver_lon_vec(ver_lon_vec>180) = ver_lon_vec(ver_lon_vec>180) - 360; %adjust longitude from [0 360] to [-180 180]
                   [ver_lat_mat,ver_lon_mat] = meshgrid(ver_lat_vec,ver_lon_vec);
                

                   control_t2m_dat = cell(length(control_file_list),1);
                   exp_t2m_dat = cell(length(exp_file_list),1);
                   ver_t2m_dat = cell(length(control_file_list),1);

                   control_rh2m_dat = cell(length(control_file_list),1);
                   exp_rh2m_dat = cell(length(exp_file_list),1);
                   ver_rh2m_dat = cell(length(control_file_list),1);

                 end % if 

                % read 2m data

                control_t2m_mat = ncread(control_file,'T2M');
                control_q2m_mat = ncread(control_file,'QV2M');
                control_ps_mat = ncread(control_file,'PS');
                exp_t2m_mat = ncread(exp_file,'T2M');
                exp_q2m_mat = ncread(exp_file,'QV2M');
                exp_ps_mat = ncread(exp_file,'PS');
                ver_t2m_mat = ncread(ver_file,'N2_metre_temperature');
                ver_t2m_mat_c = ver_t2m_mat - 273.15;
                ver_dewpoint_mat = ncread(ver_file,'N2_metre_dewpoint_temperature');
                ver_dewpoint_mat_c = ver_dewpoint_mat - 273.15;

                % convert dewpoint temperature to relative humidity 

                tmp_1 = (17.625*ver_dewpoint_mat_c)./(243.04+ver_dewpoint_mat_c);
                tmp_2 = (17.625*ver_t2m_mat_c)./(243.04+ver_t2m_mat_c);
                tmp_3 = exp(tmp_1);
                tmp_4 = exp(tmp_2);

                ver_rh2m_mat = 100*(tmp_3./tmp_4);

                % convert specific humidity to relative humidity

                tmp_1 = (exp(17.67*(control_t2m_mat-273.15)./(control_t2m_mat-29.65))).^-1;
                tmp_2 = 0.263*control_ps_mat.*control_q2m_mat;

                control_rh2m_mat = tmp_2.*tmp_1;

                tmp_1 = (exp(17.67*(exp_t2m_mat-273.15)./(exp_t2m_mat-29.65))).^-1;
                tmp_2 = 0.263*exp_ps_mat.*exp_q2m_mat;

                exp_rh2m_mat = tmp_2.*tmp_1;

                % interpolate ECMWF data

                ver_t2m_mat_interp = interp2(ver_lat_mat,ver_lon_mat,ver_t2m_mat,lat_mat,lon_mat,'spline');
                ver_rh2m_mat_interp = interp2(ver_lat_mat,ver_lon_mat,ver_rh2m_mat,lat_mat,lon_mat,'spline');        
                % store data in cells

                control_t2m_dat{f}(:,:,f_count(f)) = control_t2m_mat(lon_index,lat_index);
                exp_t2m_dat{f}(:,:,f_count(f)) = exp_t2m_mat(lon_index,lat_index);
                ver_t2m_dat{f}(:,:,f_count(f)) = ver_t2m_mat_interp(lon_index,lat_index);

                control_rh2m_dat{f}(:,:,f_count(f)) = control_rh2m_mat(lon_index,lat_index);
                exp_rh2m_dat{f}(:,:,f_count(f)) = exp_rh2m_mat(lon_index,lat_index);
                ver_rh2m_dat{f}(:,:,f_count(f)) = ver_rh2m_mat_interp(lon_index,lat_index);

            end % f

            end % if 

        end % d
      end %m
   end %y

%-----------------------------------------
% COMPUTE METRICS
%-----------------------------------------

for f = 1:length(f_count)

    % compute time series mean

    control_t2m_loc_mat = control_t2m_dat{f};
    exp_t2m_loc_mat = exp_t2m_dat{f};
    ver_t2m_loc_mat = ver_t2m_dat{f};

    control_rh2m_loc_mat = control_rh2m_dat{f};
    exp_rh2m_loc_mat = exp_rh2m_dat{f};
    ver_rh2m_loc_mat = ver_rh2m_dat{f};

    control_mean_t2m = nanmean(control_t2m_loc_mat,3);
    exp_mean_t2m = nanmean(exp_t2m_loc_mat,3);
    ver_mean_t2m = nanmean(ver_t2m_loc_mat,3);
  
    control_mean_rh2m = nanmean(control_rh2m_loc_mat,3);
    exp_mean_rh2m = nanmean(exp_rh2m_loc_mat,3);
    ver_mean_rh2m = nanmean(ver_rh2m_loc_mat,3);

    % compute anomaly matrices

    control_anom_t2m = control_t2m_loc_mat - control_mean_t2m;
    exp_anom_t2m = exp_t2m_loc_mat - exp_mean_t2m;
    ver_anom_t2m = ver_t2m_loc_mat - ver_mean_t2m;

    control_anom_rh2m = control_rh2m_loc_mat - control_mean_rh2m;
    exp_anom_rh2m = exp_rh2m_loc_mat - exp_mean_rh2m;
    ver_anom_rh2m = ver_rh2m_loc_mat - ver_mean_rh2m;

    % compute ubRMSE

    tmp =  (control_anom_t2m - ver_anom_t2m).^2 
    tmp_2 = sum(tmp,3,'omitnan')./size(control_anom_t2m,3);
    control_t2m_ubrmse(:,:,f) = sqrt(tmp_2);
    tmp = control_t2m_ubrmse(:,:,f)
    mean_control_t2m_ubrmse(f) = nanmean(tmp(:));

    tmp =  (exp_anom_t2m - ver_anom_t2m).^2
    tmp_2 = sum(tmp,3,'omitnan')./size(exp_anom_t2m,3);
    exp_t2m_ubrmse(:,:,f) = sqrt(tmp_2);
    tmp = exp_t2m_ubrmse(:,:,f)
    mean_exp_t2m_ubrmse(f) = nanmean(tmp(:));

    t2m_closeness(:,:,f) = exp_t2m_ubrmse(:,:,f) - control_t2m_ubrmse(:,:,f);

    tmp =  (control_anom_rh2m - ver_anom_rh2m).^2
    tmp_2 = sum(tmp,3,'omitnan')./size(control_anom_rh2m,3);
    control_rh2m_ubrmse(:,:,f) = sqrt(tmp_2);
    tmp = control_rh2m_ubrmse(:,:,f)
    mean_control_rh2m_ubrmse(f) = nanmean(tmp(:));

    tmp =  (exp_anom_rh2m - ver_anom_rh2m).^2
    tmp_2 = sum(tmp,3,'omitnan')./size(exp_anom_rh2m,3);
    exp_rh2m_ubrmse(:,:,f) = sqrt(tmp_2);
    tmp = exp_rh2m_ubrmse(:,:,f)
    mean_exp_rh2m_ubrmse(f) = nanmean(tmp(:));

    rh2m_closeness(:,:,f) = exp_rh2m_ubrmse(:,:,f) - control_rh2m_ubrmse(:,:,f);

    % compute anomaly correlations

    corr_anom_t2m_control_ver = NaN*ones(size(control_t2m_loc_mat,1),size(control_t2m_loc_mat,2));
    corr_anom_t2m_exp_ver = NaN*ones(size(exp_t2m_loc_mat,1),size(exp_t2m_loc_mat,2));

    corr_anom_rh2m_control_ver = NaN*ones(size(control_t2m_loc_mat,1),size(control_t2m_loc_mat,2));
    corr_anom_rh2m_exp_ver = NaN*ones(size(exp_t2m_loc_mat,1),size(exp_t2m_loc_mat,2));

    for i = 1:size(control_t2m_loc_mat,1)
        for j = 1:size(control_t2m_loc_mat,2)

            control_t2m_anom_vec = control_anom_t2m(i,j,:);
            control_t2m_anom_vec = control_t2m_anom_vec(:);
            exp_t2m_anom_vec = exp_anom_t2m(i,j,:);
            exp_t2m_anom_vec = exp_t2m_anom_vec(:);
            ver_t2m_anom_vec = ver_anom_t2m(i,j,:);
            ver_t2m_anom_vec = ver_t2m_anom_vec(:);

            control_rh2m_anom_vec = control_anom_rh2m(i,j,:);
            control_rh2m_anom_vec = control_rh2m_anom_vec(:);
            exp_rh2m_anom_vec = exp_anom_rh2m(i,j,:);
            exp_rh2m_anom_vec = exp_rh2m_anom_vec(:);
            ver_rh2m_anom_vec = ver_anom_rh2m(i,j,:);
            ver_rh2m_anom_vec = ver_rh2m_anom_vec(:);

            data_index = (isnan(control_t2m_anom_vec)+isnan(ver_t2m_anom_vec))==0;
            tmp = corrcoef(control_t2m_anom_vec(data_index),ver_t2m_anom_vec(data_index));
            corr_anom_t2m_control_ver(i,j) = tmp(1,2);
            data_index = (isnan(exp_t2m_anom_vec)+isnan(ver_t2m_anom_vec))==0;
            tmp = corrcoef(exp_t2m_anom_vec(data_index),ver_t2m_anom_vec(data_index));
            corr_anom_t2m_exp_ver(i,j) = tmp(1,2);

            data_index = (isnan(control_rh2m_anom_vec)+isnan(ver_rh2m_anom_vec))==0;
            tmp = corrcoef(control_rh2m_anom_vec(data_index),ver_rh2m_anom_vec(data_index));
            corr_anom_rh2m_control_ver(i,j) = tmp(1,2);
            data_index = (isnan(exp_rh2m_anom_vec)+isnan(ver_rh2m_anom_vec))==0;
            tmp = corrcoef(exp_rh2m_anom_vec(data_index),ver_rh2m_anom_vec(data_index));
            corr_exp_anom_rh2m_control_ver(i,j) = tmp(1,2);
        end 
    end

    mean_corr_anom_t2m_control_ver(f) = nanmean(corr_anom_t2m_control_ver(:));
    mean_corr_anom_t2m_exp_ver(f) = nanmean(corr_anom_t2m_exp_ver(:));
    mean_corr_anom_rh2m_control_ver(f) = nanmean(corr_anom_rh2m_control_ver(:));
    mean_corr_anom_rh2m_exp_ver(f) = nanmean(corr_anom_rh2m_exp_ver(:));

    corr_spat_anom_t2m_control_ver = NaN*ones(size(control_t2m_loc_mat,3),1);
    corr_spat_anom_t2m_exp_ver = NaN*ones(size(exp_t2m_loc_mat,3),1);
    corr_spat_anom_rh2m_control_ver = NaN*ones(size(control_rh2m_loc_mat,3),1);
    corr_spat_anom_rh2m_exp_ver = NaN*ones(size(exp_rh2m_loc_mat,3),1);

    for i = 1:size(control_t2m_loc_mat,3)

        control_t2m_tmp_vec = control_anom_t2m(:,:,i);
        control_t2m_tmp_vec = control_t2m_tmp_vec(:);
        exp_t2m_tmp_vec = exp_anom_t2m(:,:,i);
        exp_t2m_tmp_vec = exp_t2m_tmp_vec(:);
        ver_t2m_tmp_vec = ver_anom_t2m(:,:,i);
        ver_t2m_tmp_vec = ver_t2m_tmp_vec(:);

        control_rh2m_tmp_vec = control_anom_rh2m(:,:,i);
        control_rh2m_tmp_vec = control_rh2m_tmp_vec(:);
        exp_rh2m_tmp_vec = exp_anom_rh2m(:,:,i);
        exp_rh2m_tmp_vec = exp_rh2m_tmp_vec(:);
        ver_rh2m_tmp_vec = ver_anom_rh2m(:,:,i);
        ver_rh2m_tmp_vec = ver_rh2m_tmp_vec(:);

        data_index = (isnan(control_t2m_tmp_vec) + isnan(ver_t2m_tmp_vec)) == 0;
        tmp = corrcoef(control_t2m_tmp_vec(data_index),ver_t2m_tmp_vec(data_index));
        corr_spat_anom_t2m_control_ver(i) = tmp(1,2);
        data_index = (isnan(exp_t2m_tmp_vec) + isnan(ver_t2m_tmp_vec)) == 0;
        tmp = corrcoef(exp_t2m_tmp_vec(data_index),ver_t2m_tmp_vec(data_index));
        corr_spat_anom_t2m_exp_ver(i) = tmp(1,2);
        data_index = (isnan(control_rh2m_tmp_vec) + isnan(ver_rh2m_tmp_vec)) == 0;
        tmp = corrcoef(control_rh2m_tmp_vec(data_index),ver_rh2m_tmp_vec(data_index));
        corr_spat_anom_rh2m_control_ver(i) = tmp(1,2);
        data_index = (isnan(exp_rh2m_tmp_vec) + isnan(ver_rh2m_tmp_vec)) == 0;
        tmp = corrcoef(exp_rh2m_tmp_vec(data_index),ver_rh2m_tmp_vec(data_index)); 
        corr_spat_anom_rh2m_exp_ver(i) = tmp(1,2);

    end

    mean_corr_spat_anom_t2m_control_ver(f) = nanmean(corr_spat_anom_t2m_control_ver);
    mean_corr_spat_anom_t2m_exp_ver(f) = nanmean(corr_spat_anom_t2m_exp_ver);
    mean_corr_spat_anom_rh2m_control_ver(f) = nanmean(corr_spat_anom_rh2m_control_ver);
    mean_corr_spat_anom_rh2m_exp_ver(f) = nanmean(corr_spat_anom_rh2m_exp_ver);


end

% closeness plots

% convert longitudes back to [-180 180]

lon_vec(lon_vec>180) = lon_vec(lon_vec>180) - 360;

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

tmp_mat_t2m = t2m_closeness(:,:,:);

s1 = quantile(tmp_mat_t2m(:),0.99);
s2 = quantile(tmp_mat_t2m(:),0.01);
caxis_min_t2m = -1*max(abs(s1),abs(s2));
caxis_max_t2m = max(abs(s1),abs(s2));

cmap_mat_t2m = cmocean('balance',12);
%cmap_mat_t2m = flipud(cmap_mat_t2m);
cmap_mat_t2m(length(cmap_mat_t2m)/2:length(cmap_mat_t2m)/2+1,:) = 0.85*ones(2,3);
cmap_t2m = cmap_mat_t2m;

tmp_mat_rh2m = rh2m_closeness(:,:,:);

s1 = quantile(tmp_mat_rh2m(:),0.99);
s2 = quantile(tmp_mat_rh2m(:),0.01);
caxis_min_rh2m = -1*max(abs(s1),abs(s2));
caxis_max_rh2m = max(abs(s1),abs(s2));

cmap_mat_rh2m = cmocean('balance',12);
%cmap_mat_t2m = flipud(cmap_mat_t2m);
cmap_mat_rh2m(length(cmap_mat_rh2m)/2:length(cmap_mat_rh2m)/2+1,:) = 0.85*ones(2,3);
cmap_rh2m = cmap_mat_rh2m;

for i = 1:length(f_count)

    lead_time = i*6;

    % plot closeness maps wrt to ECMWF

    tmp_vec_t2m = tmp_mat_t2m(:,:,i);
    mean_closeness_t2m = nanmean(tmp_vec_t2m(:));
    mean_closeness_t2m = round(mean_closeness_t2m,3);

    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),t2m_closeness(:,:,i)')
    caxis([caxis_min_t2m caxis_max_t2m]);
    set(gcf,'Colormap',cmap_t2m);
    plotm(coastlat, coastlon,'k')
    hold off
    set(gca,'XLim',lonlim);
    set(gca,'YLim',latlim);
    tightmap
    title(['T2M closeness to ECMWF at ' num2str(lead_time) 'hr lead time (avg=' num2str(mean_closeness_t2m) ')'])
    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 '_prog_T2M_closeness_ubRMSE_' num2str(lead_time) 'hr'],'png')

    tmp_vec_rh2m = tmp_mat_rh2m(:,:,i);
    mean_closeness_rh2m = nanmean(tmp_vec_rh2m(:));
    mean_closeness_rh2m = round(mean_closeness_rh2m,3);

    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),rh2m_closeness(:,:,i)')
    caxis([caxis_min_rh2m caxis_max_rh2m]);
    set(gcf,'Colormap',cmap_rh2m);
    plotm(coastlat, coastlon,'k')
    hold off
    set(gca,'XLim',lonlim);
    set(gca,'YLim',latlim);
    tightmap
    title(['RH2M closeness to ECMWF at ' num2str(lead_time) 'hr lead time (avg=' num2str(mean_closeness_rh2m) ')'])
    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 '_prog_RH2M_closeness_ubRMSE_' num2str(lead_time) 'hr'],'png')
end

lead_time_vec = 6:6:126;

figure
plot(lead_time_vec,mean_corr_anom_t2m_control_ver,'LineWidth',2,'Color','m')
hold on
plot(lead_time_vec,mean_corr_anom_t2m_control_ver,'LineWidth',2,'Color','b')
hold off
grid on
xlabel('lead time [hrs]')
ylabel('anom corr T2M')
title('T2M anomaly correlation against ECMWF')
legend('control','SMAP DA exp')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_T2M_anomCorr'],'png')

figure
plot(lead_time_vec,mean_corr_anom_rh2m_control_ver,'LineWidth',2,'Color','m')
hold on
plot(lead_time_vec,mean_corr_anom_rh2m_control_ver,'LineWidth',2,'Color','b')
hold off
grid on
xlabel('lead time [hrs]')
ylabel('anom corr T2M')
title('RH2M anomaly correlation against ECMWF')
legend('control','SMAP DA exp')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_RH2M_anomCorr'],'png')

figure
plot(lead_time_vec,mean_corr_spat_anom_t2m_control_ver,'LineWidth',2,'Color','m')
hold on
plot(lead_time_vec,mean_corr_spat_anom_t2m_control_ver,'LineWidth',2,'Color','b')
hold off
grid on
xlabel('lead time [hrs]')
ylabel('anom corr T2M')
title('T2M anomaly correlation against ECMWF')
legend('control','SMAP DA exp')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_T2M_spatAnomCorr'],'png')

figure
plot(lead_time_vec,mean_corr_spat_anom_rh2m_control_ver,'LineWidth',2,'Color','m')
hold on
plot(lead_time_vec,mean_corr_spat_anom_rh2m_control_ver,'LineWidth',2,'Color','b')
hold off
grid on
xlabel('lead time [hrs]')
ylabel('anom corr T2M')
title('RH2M anomaly correlation against ECMWF')
legend('control','SMAP DA exp')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_RH2M_spatAnomCorr'],'png')

mean_spat_corr = [mean_corr_spat_anom_t2m_control_ver;mean_corr_spat_anom_t2m_exp_ver];

figure
b1=bar(lead_time_vec,mean_spat_corr);
grid on
set(b1(1),'FaceColor',[0.5 0.5 0.5])
set(b1(2),'FaceColor',[0.1 0.1 0.1])
set(gca,'YLim',[0.3 0.7])
xlabel('lead time [hrs]')
ylabel('anom corr T2M')
title('T2M anomaly correlation against ECMWF')
legend('control','SMAP DA exp')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_T2M_spatAnomCorr_bar'],'png')

mean_spat_corr_rh = [mean_corr_spat_anom_rh2m_control_ver;mean_corr_spat_anom_rh2m_exp_ver];

figure
b1=bar(lead_time_vec,mean_spat_corr_rh);
grid on
set(b1(1),'FaceColor',[0.5 0.5 0.5])
set(b1(2),'FaceColor',[0.1 0.1 0.1])
set(gca,'YLim',[0.39 0.71])
xlabel('lead time [hrs]')
ylabel('anom corr RH2M')
title('RH2M anomaly correlation against ECMWF')
legend('control','SMAP DA exp')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_RH2M_spatAnomCorr_bar'],'png')

mean_t2m_ubrmse = [mean_control_t2m_ubrmse; mean_exp_t2m_ubrmse];

figure
b1=bar(lead_time_vec,mean_t2m_ubrmse);
grid on
set(b1(1),'FaceColor',[0.5 0.5 0.5])
set(b1(2),'FaceColor',[0.1 0.1 0.1])
set(gca,'YLim',[0.8 1.6])
xlabel('lead time [hrs]')
ylabel('ubRMSE T2M [K]')
title('T2M ubRMSE against ECMWF')
legend('control','SMAP DA exp')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_T2M_ubRMSE_bar'],'png')

mean_rh2m_ubrmse = [mean_control_rh2m_ubrmse; mean_exp_rh2m_ubrmse];

figure
b1=bar(lead_time_vec,mean_rh2m_ubrmse);
grid on
set(b1(1),'FaceColor',[0.5 0.5 0.5])
set(b1(2),'FaceColor',[0.1 0.1 0.1])
set(gca,'YLim',[5 9])
xlabel('lead time [hrs]')
ylabel('ubRMSE RH2M [%]')
title('RH2M ubRMSE against ECMWF')
legend('control','SMAP DA exp')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_RH2M_ubRMSE_bar'],'png')
