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/emcgrath/archive/'; 
control_dir = ['/discover/nobackup/projects/smaptc/jkolassa/archive/' storm_name '/'];
%control_name = 'Gulab_CNTRL_4DVAR'; 
control_name =  'LADAS_control_4D_new';
%exp_dir = '/discover/nobackup/projects/smaptc/emcgrath/archive/';
exp_dir = ['/discover/nobackup/projects/smaptc/jkolassa/archive/' storm_name '/'];
%exp_name = 'Gulab_SMAP_4DVAR';
exp_name = 'Idai_SMAP_4D';
ecmwf_dir = '/discover/nobackup/projects/gmao/g6dev/ltakacs/ecmwf/';
ecmwf_name = 'ecmwf';

%lat_range = [14 18]; 
lat_range = [-24 -10];
%lon_range = [52 97]; ]
lon_range = [24 52];

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

%start_year = 2021;
%start_month = 9;
%start_day = 24;

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

%end_year = 2021;
%end_month = 10;
%end_day = 4;

%---------------------------------------------
% 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) & ...
                (length(control_3D_list)==21) & (length(exp_3D_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];

                control_3D_file = [control_dir_name control_3D_list(f).name];
                exp_3D_file = [exp_dir_name exp_3D_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'));
                   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_mse_mat = zeros(length(lon_vec),length(lat_vec),21);
                   exp_t2m_mse_mat = zeros(length(lon_vec),length(lat_vec),21);
                   control_rh2m_mse_mat = zeros(length(lon_vec),length(lat_vec),21);
                   exp_rh2m_mse_mat = zeros(length(lon_vec),length(lat_vec),21);

                   control_tlml_mse_mat = zeros(length(lon_vec),length(lat_vec),21);
                   exp_tlml_mse_mat = zeros(length(lon_vec),length(lat_vec),21);
                   control_rhlml_mse_mat = zeros(length(lon_vec),length(lat_vec),21);
                   exp_rhlml_mse_mat = zeros(length(lon_vec),length(lat_vec),21);
                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;
               
                % read lowest model layer data

                control_level = ncread(control_3D_file,'lev');
                control_lev_index = (control_level==850); 
                ver_level = ncread(ver_file,'lev');
                ver_lev_index = (ver_level==850);

                tmp_mat = ncread(control_3D_file,'T');
                control_tlml_mat = tmp_mat(:,:,control_lev_index);
                tmp_mat = ncread(exp_3D_file,'T');
                exp_tlml_mat = tmp_mat(:,:,control_lev_index);
                tmp_mat = ncread(ver_file,'Temperature');
                ver_tlml_mat = tmp_mat(:,:,ver_lev_index);

                tmp_mat = ncread(control_3D_file,'RH');
                control_rhlml_mat = tmp_mat(:,:,control_lev_index);
                tmp_mat = ncread(exp_3D_file,'RH');
                exp_rhlml_mat = tmp_mat(:,:,control_lev_index);
                tmp_mat = ncread(ver_file,'Relative_humidity');
                ver_rhlml_mat = tmp_mat(:,:,ver_lev_index);

                % 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');         
                ver_tlml_mat_interp = interp2(ver_lat_mat,ver_lon_mat,ver_tlml_mat,lat_mat,lon_mat,'spline');
                ver_rhlml_mat_interp = interp2(ver_lat_mat,ver_lon_mat,ver_rhlml_mat,lat_mat,lon_mat,'spline');

                % compute sum MSE

                tmp_mat = (control_t2m_mat - ver_t2m_mat_interp).^2;
                control_t2m_mse_mat(:,:,f) = control_t2m_mse_mat(:,:,f) + tmp_mat; 

                tmp_mat = (exp_t2m_mat - ver_t2m_mat_interp).^2;
                exp_t2m_mse_mat(:,:,f) = exp_t2m_mse_mat(:,:,f) + tmp_mat;

                tmp_mat = (control_rh2m_mat - ver_rh2m_mat_interp).^2;
                control_rh2m_mse_mat(:,:,f) = control_rh2m_mse_mat(:,:,f) + tmp_mat;

                tmp_mat = (exp_rh2m_mat - ver_rh2m_mat_interp).^2;
                exp_rh2m_mse_mat(:,:,f) = exp_rh2m_mse_mat(:,:,f) + tmp_mat;
 
                tmp_mat = (control_tlml_mat - ver_tlml_mat_interp).^2;
                control_tlml_mse_mat(:,:,f) = control_tlml_mse_mat(:,:,f) + tmp_mat;

                tmp_mat = (exp_tlml_mat - ver_tlml_mat_interp).^2;
                exp_tlml_mse_mat(:,:,f) = exp_tlml_mse_mat(:,:,f) + tmp_mat;

                tmp_mat = (control_rhlml_mat - ver_rhlml_mat_interp).^2;
                control_rhlml_mse_mat(:,:,f) = control_rhlml_mse_mat(:,:,f) + tmp_mat;

                tmp_mat = (exp_rhlml_mat - ver_rhlml_mat_interp).^2;
                exp_rhlml_mse_mat(:,:,f) = exp_rhlml_mse_mat(:,:,f) + tmp_mat;

            end % f

            end % if 

        end % d
      end %m
   end %y

% convert longitudes back to [-180 180]

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

[lat_mat,lon_mat] = meshgrid(lat_vec,lon_vec);
[ver_lat_mat,ver_lon_mat] = meshgrid(ver_lat_vec,ver_lon_vec);

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;

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

for i = 1:length(f_count)

    control_t2m_rmse_mat(:,:,i) = sqrt(control_t2m_mse_mat(:,:,i)./f_count(i));
    exp_t2m_rmse_mat(:,:,i) = sqrt(exp_t2m_mse_mat(:,:,i)./f_count(i));

    control_rh2m_rmse_mat(:,:,i) = sqrt(control_rh2m_mse_mat(:,:,i)./f_count(i));
    exp_rh2m_rmse_mat(:,:,i) = sqrt(exp_rh2m_mse_mat(:,:,i)./f_count(i));

    t2m_closeness_mat(:,:,i) = exp_t2m_rmse_mat(:,:,i) - control_t2m_rmse_mat(:,:,i);
    rh2m_closeness_mat(:,:,i) = exp_rh2m_rmse_mat(:,:,i) - control_rh2m_rmse_mat(:,:,i);

    control_tlml_rmse_mat(:,:,i) = sqrt(control_tlml_mse_mat(:,:,i)./f_count(i));
    exp_tlml_rmse_mat(:,:,i) = sqrt(exp_tlml_mse_mat(:,:,i)./f_count(i));

    control_rhlml_rmse_mat(:,:,i) = sqrt(control_rhlml_mse_mat(:,:,i)./f_count(i));
    exp_rhlml_rmse_mat(:,:,i) = sqrt(exp_rhlml_mse_mat(:,:,i)./f_count(i));

    tlml_closeness_mat(:,:,i) = exp_tlml_rmse_mat(:,:,i) - control_tlml_rmse_mat(:,:,i);
    rhlml_closeness_mat(:,:,i) = exp_rhlml_rmse_mat(:,:,i) - control_rhlml_rmse_mat(:,:,i);
end

%-----------------------------------------
% PLOT METRICS
%-----------------------------------------

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

tmp_mat_t2m = t2m_closeness_mat(lon_index,lat_index,:);

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_mat(lon_index,lat_index,:);

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;

tmp_mat_tlml = tlml_closeness_mat(lon_index,lat_index,:);

s1 = quantile(tmp_mat_tlml(:),0.99);
s2 = quantile(tmp_mat_tlml(:),0.01);
caxis_min_tlml = -1*max(abs(s1),abs(s2));
caxis_max_tlml = max(abs(s1),abs(s2));

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

tmp_mat_rhlml = rhlml_closeness_mat(lon_index,lat_index,:);

s1 = quantile(tmp_mat_rhlml(:),0.99);
s2 = quantile(tmp_mat_rhlml(:),0.01);
caxis_min_rhlml = -1*max(abs(s1),abs(s2));
caxis_max_rhlml = max(abs(s1),abs(s2));

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


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,lon_vec,t2m_closeness_mat(:,:,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_' 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,lon_vec,rh2m_closeness_mat(:,:,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_' num2str(lead_time) 'hr'],'png')

    tmp_vec_tlml = tmp_mat_tlml(:,:,i);
    mean_closeness_tlml = nanmean(tmp_vec_tlml(:));
    mean_closeness_tlml = round(mean_closeness_tlml,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,lon_vec,tlml_closeness_mat(:,:,i)')
    caxis([caxis_min_tlml caxis_max_tlml]);
    set(gcf,'Colormap',cmap_tlml);
    plotm(coastlat, coastlon,'k')
    hold off
    set(gca,'XLim',lonlim);
    set(gca,'YLim',latlim);
    tightmap
    title(['T@850mb closeness to ECMWF at ' num2str(lead_time) 'hr lead time (avg=' num2str(mean_closeness_tlml) ')'])
    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_TLML_closeness_' num2str(lead_time) 'hr'],'png')

    tmp_vec_rhlml = tmp_mat_rhlml(:,:,i);
    mean_closeness_rhlml = nanmean(tmp_vec_rhlml(:));
    mean_closeness_rhlml = round(mean_closeness_rhlml,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,lon_vec,rhlml_closeness_mat(:,:,i)')
    caxis([caxis_min_rhlml caxis_max_rhlml]);
    set(gcf,'Colormap',cmap_rhlml);
    plotm(coastlat, coastlon,'k')
    hold off
    set(gca,'XLim',lonlim);
    set(gca,'YLim',latlim);
    tightmap
    title(['RH@850mb closeness to ECMWF at ' num2str(lead_time) 'hr lead time (avg=' num2str(mean_closeness_rhlml) ')'])
    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_RHLML_closeness_' num2str(lead_time) 'hr'],'png')


    % compute area-average RMSE

    tmp_vec = control_t2m_rmse_mat(lon_index,lat_index,i);
    control_t2m_rmse_aavg(i) = nanmean(tmp_vec(:));
    tmp_vec = control_rh2m_rmse_mat(lon_index,lat_index,i);
    control_rh2m_rmse_aavg(i) = nanmean(tmp_vec(:));

    tmp_vec = exp_t2m_rmse_mat(lon_index,lat_index,i);
    exp_t2m_rmse_aavg(i) = nanmean(tmp_vec(:));
    tmp_vec = exp_rh2m_rmse_mat(lon_index,lat_index,i);
    exp_rh2m_rmse_aavg(i) = nanmean(tmp_vec(:));

    tmp_vec = control_tlml_rmse_mat(lon_index,lat_index,i);
    control_tlml_rmse_aavg(i) = nanmean(tmp_vec(:));
    tmp_vec = control_rhlml_rmse_mat(lon_index,lat_index,i);
    control_rhlml_rmse_aavg(i) = nanmean(tmp_vec(:));

    tmp_vec = exp_tlml_rmse_mat(lon_index,lat_index,i);
    exp_tlml_rmse_aavg(i) = nanmean(tmp_vec(:));
    tmp_vec = exp_rhlml_rmse_mat(lon_index,lat_index,i);
    exp_rhlml_rmse_aavg(i) = nanmean(tmp_vec(:));
end

lead_time_vec = 6:6:126;
full_day_index = [4:4:length(lead_time_vec)];

mean_rmse_t2m = [control_t2m_rmse_aavg; exp_t2m_rmse_aavg];

figure
b1=bar(lead_time_vec,mean_rmse_t2m);
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',[1.2 2])
xlabel('lead time [hrs]')
ylabel('RMSE T2M [K]')
title('T2M RMSE against ECMWF')
legend('CNTRL\_FCST','SMAP\_FCST','Location','northwest')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_T2M_RMSE_bar'],'png')

day_vec = lead_time_vec(full_day_index)/24;

figure
b1=bar(day_vec,mean_rmse_t2m(:,full_day_index));
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',[1.2 1.5])
set(gca,'FontSize',14)
xlabel('Forecast lead time [days]','FontSize',14')
ylabel('RMSE T2M [K]','FontSize',14')
title('T2M RMSE against ECMWF','FontSize',14')
legend('CNTRL\_FCST','SMAP\_FCST','Location','northwest', 'FontSize',14)
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_T2M_RMSE_bar_fullDay'],'png')


mean_rmse_rh2m = [control_rh2m_rmse_aavg; exp_rh2m_rmse_aavg];

figure
b1=bar(lead_time_vec,mean_rmse_rh2m);
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',[6 10])
xlabel('lead time [hrs]')
ylabel('RMSE RH2M [%]')
title('RH2M RMSE against ECMWF')
legend('CNTRL\_FCST','SMAP\_FCST','Location','northwest')
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_RH2M_RMSE_bar'],'png')

day_vec = lead_time_vec(full_day_index)/24;

figure
b1=bar(day_vec,mean_rmse_rh2m(:,full_day_index));
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',[6 9])
set(gca,'FontSize',14')
xlabel('Forecast lead time [days]','FontSize',14')
ylabel('RMSE RH2M [%]','FontSize',14')
title('RH2M RMSE against ECMWF','FontSize',14')
legend('CNTRL\_FCST','SMAP\_FCST','Location','northwest', 'FontSize',14)
set(gcf,'paperunits','centimeters')
set(gcf,'paperposition',[1 1 17 15])
saveas(gcf,[save_path storm_name '_prog_RH2M_RMSE_bar_fullDay'],'png')


figure 
plot(lead_time_vec,control_t2m_rmse_aavg,'LineWidth',2,'Color','m')
hold on
plot(lead_time_vec,exp_t2m_rmse_aavg,'LineWidth',2,'Color','b')
hold off
grid on
xlabel('lead time [hrs]')
ylabel('T2M RMSE [K]')
title('area-average T2M RMSE [K] 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_RMSE_aavg'],'png')

figure
plot(lead_time_vec,control_rh2m_rmse_aavg,'LineWidth',2,'Color','m')
hold on
plot(lead_time_vec,exp_rh2m_rmse_aavg,'LineWidth',2,'Color','b')
hold off
grid on
xlabel('lead time [hrs]')
ylabel('RH2M RMSE [K]')
title('area-average RH2M RMSE [-] 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_RMSE_aavg'],'png')

figure
plot(lead_time_vec,control_tlml_rmse_aavg,'LineWidth',2,'Color','m')
hold on
plot(lead_time_vec,exp_tlml_rmse_aavg,'LineWidth',2,'Color','b')
hold off
grid on
xlabel('lead time [hrs]')
ylabel('T RMSE [K]')
title('area-average T@850mb RMSE [K] 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_TLML_RMSE_aavg'],'png')

figure
plot(lead_time_vec,control_rhlml_rmse_aavg,'LineWidth',2,'Color','m')
hold on
plot(lead_time_vec,exp_rhlml_rmse_aavg,'LineWidth',2,'Color','b')
hold off
grid on
xlabel('lead time [hrs]')
ylabel('RH RMSE [K]')
title('area-average RH@850mb RMSE [-] 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_RHLML_RMSE_aavg'],'png')
