clear all;
close all;

format long
scrnsz = get(0,'ScreenSize');
fontsize = 16;
line_width_plot = 1.5;
set(0,'DefaultAxesFontSize',fontsize);
set(0,'DefaultFigureColor','White');
set(0,'DefaultFigurePosition',[100,100,0.6*(scrnsz(3)-100),0.9*(scrnsz(4)-100)]');
%set(0,'DefaultFigurePosition',[100,100,0.65*(scrnsz(3)-100),0.9*(scrnsz(4)-100)]');
%set(0,'DefaultFigurePosition',[100,100,0.65*(scrnsz(3)-100),0.9*(scrnsz(4)-100)]');
set(0,'DefaultLineLineWidth',line_width_plot);

% Germany data (not part of original data set, downloaded from 
% https://gadm.org/download_country.html )
addpath('/Volumes/DATA3/FESSTHH2020data/gadm40_DEU_shp');

month_to_examine = '08';
day_to_examine = '10';
apollo_folder_name = strcat('/Volumes/DATA3/FESSTHH2020data/FESSTHH2020_',month_to_examine,'_APOLLO_L2/',day_to_examine,'/');
wxt_folder_name = strcat('/Volumes/DATA3/FESSTHH2020data/FESSTHH2020_',month_to_examine,'_WXT_L2/',day_to_examine,'/');

% For Figure 9 reproduction:
stations_to_plot_apollo = {'018LMa','054OGa','082PGa','104PGa'};
stations_to_plot_wxt = {};
hour_to_start_plotting = 12;
minute_to_start_plotting = 30;
hour_to_end_plotting = 15;
minute_to_end_plotting = 0;

% For gridded data:
hour_to_start_plotting_grid = 12;
minute_to_start_plotting_grid = 30;
hour_to_end_plotting_grid = 15;
minute_to_end_plotting_grid = 0;
plot_interval_minutes_grid = 5;
hour_reference_grid = 12;
minute_reference_grid = 30;

% TASKS TO RUN:
% 0 = do not do this
% 1 = do this
reproduce_figure_9 = 1;
plot_all_temps = 0;
plot_all_pressures = 0;
plot_wxt_hail = 0;
plot_wxt_precip = 0;
plot_wxt_rv = 0;
plot_wxt_winds = 0;

darkblue = [1 0 137]./255;
lightblue = [135 206 235]./255;
redcolor = [220 20 60]./255;
orange = [255 165 0]./255;
set(groot,'DefaultAxesColorOrder',[redcolor;orange;lightblue;darkblue])



apollo_vars = {'pa','ta'};
num_apollo_vars = numel(apollo_vars);

wxt_vars = {'hail','precip','hur','pa','ta','wdir','wspeed_max','wspeed'};
num_wxt_vars = numel(wxt_vars);

sample_nc_file_apollo = strcat(apollo_folder_name,'fessthh_uhh_apollo_l2_pa_v00_2020',month_to_examine,day_to_examine,'000000.nc');
sample_nc_file_wxt = strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_pa_v00_2020',month_to_examine,day_to_examine,'000000.nc');
%ncdisp(strcat(apollo_folder_name,'fessthh_uhh_apollo_l2_pa_v00_2020',month_to_examine,day_to_examine,'000000.nc'));
%ncdisp(strcat(apollo_folder_name,'fessthh_uhh_apollo_l2_ta_v00_2020',month_to_examine,day_to_examine,'000000.nc'));

%ncdisp(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_pa_v00_2020',month_to_examine,day_to_examine,'000000.nc'));
%ncdisp(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_ta_v00_2020',month_to_examine,day_to_examine,'000000.nc'));

%start_datetime = datetime(1970,1,1,0,0,0,0,'TimeZone','UTCLeapSeconds');
start_datetime = datetime(1970,1,1,0,0,0,0,'TimeZone','UTC');
% Times
timedata_apollo = ncread(sample_nc_file_apollo,'time');
timedata_wxt = ncread(sample_nc_file_wxt,'time');
% Start and end of time averaging intervals: (2 x num times)
time_bnds_apollo = ncread(sample_nc_file_apollo,'time_bnds');
time_bnds_wxt = ncread(sample_nc_file_wxt,'time_bnds');
% Station identifier code:
station_id_apollo = ncread(sample_nc_file_apollo,'station_id');
station_id_wxt = ncread(sample_nc_file_wxt,'station_id');
% Station name:
station_name_apollo = ncread(sample_nc_file_apollo,'station_name');
station_name_wxt = ncread(sample_nc_file_wxt,'station_name');
% Latitude:
lat_apollo = ncread(sample_nc_file_apollo,'lat');
lat_wxt = ncread(sample_nc_file_wxt,'lat');
% Longitude:
lon_apollo = ncread(sample_nc_file_apollo,'lon');
lon_wxt = ncread(sample_nc_file_wxt,'lon');
% Number of stations:
num_stations_apollo = numel(lat_apollo);
num_stations_wxt = numel(lat_wxt);
% Altitude of ground level relative to mean sea level:
zsl_apollo = ncread(sample_nc_file_apollo,'zsl');
zsl_wxt = ncread(sample_nc_file_wxt,'zsl');
% Height of sensor above ground:
zag_apollo = ncread(sample_nc_file_apollo,'zag');
zag_wxt = ncread(sample_nc_file_wxt,'zag');
% Local climate zone of station: (40 x num_stations)
lcz_apollo = ncread(sample_nc_file_apollo,'lcz');
lcz_wxt = ncread(sample_nc_file_wxt,'lcz');


full_datetimes_apollo = start_datetime + seconds(timedata_apollo);
full_datetimes_wxt = start_datetime + seconds(timedata_wxt);

% Read in the meteorological data:

for apollo_var_num = 1:num_apollo_vars
    apollo_var_name = apollo_vars{1,apollo_var_num};
    switch apollo_var_name
        case 'pa' % Pressure in Pa
            pa_apollo = ncread(strcat(apollo_folder_name,'fessthh_uhh_apollo_l2_',apollo_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),apollo_var_name);
        case 'ta' % Temperature in K
            ta_apollo = ncread(strcat(apollo_folder_name,'fessthh_uhh_apollo_l2_',apollo_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),apollo_var_name);
        otherwise
            error(strcat('Need to define what to do with the variable named:',32,apollo_var_name));
    end
end

for wxt_var_num = 1:num_wxt_vars
    wxt_var_name = wxt_vars{1,wxt_var_num};
    switch wxt_var_name
        case 'hail' % Amount of hail: Number of hail hits within the measurement interval, m^-2
            hail_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
        case 'precip' % Accumulated amount of precipitation within measurement interval, kg/m^2
            precip_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
        case 'hur' % Relative humidity (fraction)
            hur_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
        case 'pa'
            pa_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
        case 'ta'
            ta_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
        case 'wdir' % Wind direction (average of 4-Hz data within measurement interval), degrees
            wdir_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
        case 'wspeed_max' % Wind gust: maximum 3-s average of 4-Hz data within measurement interval, m/s
            wspeed_max_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
        case 'wspeed' % Wind speed: average of 4-Hz data within measurement interval, m/s
            wspeed_wxt = ncread(strcat(wxt_folder_name,'fessthh_uhh_wxt_l2_',wxt_var_name,'_v00_2020',month_to_examine,day_to_examine,'000000.nc'),wxt_var_name);
        otherwise
            error(strcat('Need to define what to do with the variable named:',32,wxt_var_name));
    end
end

% Identify the indices corresponding to the apollo stations to plot
num_stations_to_plot = numel(stations_to_plot_apollo);
station_id_apollo_short = station_id_apollo(1:6,:);
station_id_apollo_short_rotated = station_id_apollo_short';
stations_to_plot_indices_apollo = zeros(1,num_stations_to_plot);
for station_to_plot = 1:num_stations_to_plot
    for station_to_compare = 1:num_stations_apollo
        if strcmp(stations_to_plot_apollo{1,station_to_plot},station_id_apollo_short_rotated(station_to_compare,:))
            stations_to_plot_indices_apollo(1,station_to_plot) = station_to_compare;
            break;
        end
    end
    if stations_to_plot_indices_apollo(1,station_to_plot) == 0
        error(strcat('Could not find apollo station',32,stations_to_plot_apollo{1,station_to_plot},', #',num2str(station_to_plot)));
    end
end

% Identify the indices corresponding to the wxt stations to plot
num_stations_to_plot = numel(stations_to_plot_wxt);
station_id_wxt_short = station_id_wxt(1:6,:);
station_id_wxt_short_rotated = station_id_wxt_short';
stations_to_plot_indices_wxt = zeros(1,num_stations_to_plot);
for station_to_plot = 1:num_stations_to_plot
    for station_to_compare = 1:num_stations_wxt
        if strcmp(stations_to_plot_wxt{1,station_to_plot},station_id_wxt_short_rotated(station_to_compare,:))
            stations_to_plot_indices_wxt(1,station_to_plot) = station_to_compare;
            break;
        end
    end
    if stations_to_plot_indices_wxt(1,station_to_plot) == 0
        error(strcat('Could not find wxt station',32,stations_to_plot_wxt{1,station_to_plot},', #',num2str(station_to_plot)));
    end
end

% Reproduce Figure 9:
start_datetime_plot = datetime(2020,str2double(month_to_examine),str2double(day_to_examine),...
    hour_to_start_plotting,minute_to_start_plotting,0,0,'TimeZone','UTC');
end_datetime_plot = datetime(2020,str2double(month_to_examine),str2double(day_to_examine),...
    hour_to_end_plotting,minute_to_end_plotting,0,0,'TimeZone','UTC');

start_plot_time_index_apollo = find(full_datetimes_apollo == start_datetime_plot);
end_plot_time_index_apollo = find(full_datetimes_apollo == end_datetime_plot);

start_plot_time_index_wxt = find(full_datetimes_wxt == start_datetime_plot);
end_plot_time_index_wxt = find(full_datetimes_wxt == end_datetime_plot);

datetime_plot_apollo = full_datetimes_apollo(start_plot_time_index_apollo:end_plot_time_index_apollo);
datetime_plot_wxt = full_datetimes_wxt(start_plot_time_index_wxt:end_plot_time_index_wxt);

ta_apollo_to_plot = -273.15 + ta_apollo(stations_to_plot_indices_apollo,start_plot_time_index_apollo:end_plot_time_index_apollo);
pa_apollo_to_plot = 0.01.*pa_apollo(stations_to_plot_indices_apollo,start_plot_time_index_apollo:end_plot_time_index_apollo);

ta_apollo_to_plot = ta_apollo_to_plot - repmat(ta_apollo_to_plot(:,1),1,numel(start_plot_time_index_apollo:end_plot_time_index_apollo));
pa_apollo_to_plot = pa_apollo_to_plot - repmat(pa_apollo_to_plot(:,1),1,numel(start_plot_time_index_apollo:end_plot_time_index_apollo));

ta_wxt_to_plot = -273.15 + ta_wxt(stations_to_plot_indices_wxt,start_plot_time_index_wxt:end_plot_time_index_wxt);
pa_wxt_to_plot = 0.01.*pa_wxt(stations_to_plot_indices_wxt,start_plot_time_index_wxt:end_plot_time_index_wxt);

ta_wxt_to_plot = ta_wxt_to_plot - repmat(ta_wxt_to_plot(:,1),1,numel(start_plot_time_index_wxt:end_plot_time_index_wxt));
pa_wxt_to_plot = pa_wxt_to_plot - repmat(pa_wxt_to_plot(:,1),1,numel(start_plot_time_index_wxt:end_plot_time_index_wxt));

if reproduce_figure_9
    figure(1)
    subplot(2,1,1)
    plot(datetime_plot_apollo,ta_apollo_to_plot)
    %xlabel('Time (UTC)');
    ylabel('Temperature perturbation (K)');
    legend(stations_to_plot_apollo,'Location','SouthWest');
    %legend boxoff;
    ylim([-14 2])
    grid on;
    subplot(2,1,2)
    plot(datetime_plot_apollo,pa_apollo_to_plot)
    xlabel('Time (UTC)');
    ylabel('Pressure perturbation (hPa)');
    %legend(stations_to_plot_apollo,'Location','SouthWest');
    %legend boxoff;
    ylim([-1.5 1.5])
    grid on;
end
% Interpolate to a grid:
if plot_wxt_winds
    loni = 9.64:0.02:10.40;
    lati = 53.35:0.01:53.68;
else
    loni = 9.64:0.001:10.40;
    lati = 53.35:0.001:53.68;
end
[lati_mesh,loni_mesh]=meshgrid(lati,loni);
lat_all = double([lat_apollo;lat_wxt]);
lon_all = double([lon_apollo;lon_wxt]);

start_datetime_plot_grid = datetime(2020,str2double(month_to_examine),str2double(day_to_examine),...
    hour_to_start_plotting_grid,minute_to_start_plotting_grid,0,0,'TimeZone','UTC');
end_datetime_plot_grid = datetime(2020,str2double(month_to_examine),str2double(day_to_examine),...
    hour_to_end_plotting_grid,minute_to_end_plotting_grid,0,0,'TimeZone','UTC');

datetime_plot_apollo_grid = start_datetime_plot_grid:minutes(plot_interval_minutes_grid):end_datetime_plot_grid;
datetime_plot_wxt_grid = start_datetime_plot_grid:minutes(plot_interval_minutes_grid):end_datetime_plot_grid;

start_plot_time_index_apollo_grid = find(full_datetimes_apollo == start_datetime_plot_grid);
end_plot_time_index_apollo_grid = find(full_datetimes_apollo == end_datetime_plot_grid);

start_plot_time_index_wxt_grid = find(full_datetimes_wxt == start_datetime_plot_grid);
end_plot_time_index_wxt_grid = find(full_datetimes_wxt == end_datetime_plot_grid);

time_spacing_index_apollo_grid = plot_interval_minutes_grid*60;
time_spacing_index_wxt_grid = plot_interval_minutes_grid*6;

indices_plot_apollo_grid = start_plot_time_index_apollo_grid: ...
    time_spacing_index_apollo_grid:end_plot_time_index_apollo_grid;

indices_plot_wxt_grid = start_plot_time_index_wxt_grid: ...
    time_spacing_index_wxt_grid:end_plot_time_index_wxt_grid;

num_indices_plot_grid = numel(indices_plot_apollo_grid);
if num_indices_plot_grid ~= numel(indices_plot_wxt_grid)
    error('The number of apollo time indices does not match the number of wxt time indices for gridded plots');
end

datetime_reference_grid = datetime(2020,str2double(month_to_examine),str2double(day_to_examine),...
    hour_reference_grid,minute_reference_grid,0,0,'TimeZone','UTC');

reference_plot_time_index_apollo_grid = find(full_datetimes_apollo == datetime_reference_grid);
reference_plot_time_index_wxt_grid = find(full_datetimes_wxt == datetime_reference_grid);

germany_map_data = shaperead('/Volumes/DATA3/FESSTHH2020data/gadm40_DEU_shp/gadm40_DEU_1.shp', 'UseGeoCoords', true);

% Plot temperature perturbation:
if plot_all_temps
    ref_temp_data_apollo = ta_apollo(:,reference_plot_time_index_apollo_grid);
    ref_temp_data_wxt = ta_wxt(:,reference_plot_time_index_wxt_grid);
    ref_temp_data_all = [ref_temp_data_apollo;ref_temp_data_wxt];
    ref_indices_to_include = ~isnan(ref_temp_data_all);
    ref_temp_interpolant = scatteredInterpolant(lat_all(ref_indices_to_include),...
        lon_all(ref_indices_to_include),ref_temp_data_all(ref_indices_to_include),'linear','none');
    ref_temp_gridded = ref_temp_interpolant(lati_mesh,loni_mesh);
    for t_ind = 1:num_indices_plot_grid
        t_ind_of_entire_day_apollo = indices_plot_apollo_grid(t_ind);
        t_ind_of_entire_day_wxt = indices_plot_wxt_grid(t_ind);
        temp_data_apollo = ta_apollo(:,t_ind_of_entire_day_apollo);
        temp_data_wxt = ta_wxt(:,t_ind_of_entire_day_wxt);
        temp_data_all = [temp_data_apollo;temp_data_wxt];
        indices_to_include = ~isnan(temp_data_all);
        temp_interpolant = scatteredInterpolant(lat_all(indices_to_include),...
            lon_all(indices_to_include),temp_data_all(indices_to_include),'linear','none');
        temp_gridded = temp_interpolant(lati_mesh,loni_mesh);
        temp_perturbation_gridded = temp_gridded - ref_temp_gridded;

        fig2 = figure(2);
        ax = worldmap('Germany');
        C1=pcolorm(lati_mesh,loni_mesh,temp_perturbation_gridded);
        caxis([-14 4]);
        setm(ax,'MapLatLimit',[min(lati),max(lati)]);
        setm(ax,'MapLonLimit',[min(loni),max(loni)]);
        gridm('MeridianLabel','on', ...
            'ParallelLabel','on',...
            'glinestyle',':','gcolor','k','glinewidth',1,...
            'mlinelocation',9.7:0.1:10.4,'plinelocation',53.4:0.1:53.7)
        setm(ax,'FontSize',16,...
            'MLabelLocation',9.7:0.1:10.4,'MLabelRound',-1,...
            'PLabelLocation',53.4:0.1:53.7,'PLabelRound',-1)
        geoshow(germany_map_data,'FaceColor','none')
        str1 = char(string(datetime_reference_grid));
        str2 = char(string(datetime_plot_apollo_grid(t_ind)));
        title({strcat('Temperature perturbation (K) relative to',32,str1,32,'UTC');strcat('at',32,str2,32,'UTC')});
        colorbar EastOutside
        image_here = getframe(fig2);
        imwrite(image_here.cdata,...
            strcat('/Volumes/DATA3/FESSTHH2020data/images/temp_pert_',...
            string(datetime_plot_apollo_grid(t_ind)),'.png'),'png');
        gifname = '/Volumes/DATA3/FESSTHH2020data/images/temp_pert.gif';
        gif_delaytime = 0.75;
        [A,map] = rgb2ind(image_here.cdata,256);
        if t_ind == 1
            imwrite(A,map,gifname,'gif','LoopCount',Inf,'DelayTime',gif_delaytime);
        else
            imwrite(A,map,gifname,'gif','WriteMode','append','DelayTime',gif_delaytime);
        end
    end
end
% Plot pressure perturbation:
if plot_all_pressures
    ref_press_data_apollo = pa_apollo(:,reference_plot_time_index_apollo_grid);
    ref_press_data_wxt = pa_wxt(:,reference_plot_time_index_wxt_grid);
    ref_press_data_all = [ref_press_data_apollo;ref_press_data_wxt];
    ref_indices_to_include = ~isnan(ref_press_data_all);
    ref_press_interpolant = scatteredInterpolant(lat_all(ref_indices_to_include),...
        lon_all(ref_indices_to_include),ref_press_data_all(ref_indices_to_include),'linear','none');
    ref_press_gridded = ref_press_interpolant(lati_mesh,loni_mesh);
    for t_ind = 1:num_indices_plot_grid
        t_ind_of_entire_day_apollo = indices_plot_apollo_grid(t_ind);
        t_ind_of_entire_day_wxt = indices_plot_wxt_grid(t_ind);
        press_data_apollo = pa_apollo(:,t_ind_of_entire_day_apollo);
        press_data_wxt = pa_wxt(:,t_ind_of_entire_day_wxt);
        press_data_all = [press_data_apollo;press_data_wxt];
        indices_to_include = ~isnan(press_data_all);
        press_interpolant = scatteredInterpolant(lat_all(indices_to_include),...
            lon_all(indices_to_include),press_data_all(indices_to_include),'linear','none');
        press_gridded = press_interpolant(lati_mesh,loni_mesh);
        press_perturbation_gridded = press_gridded - ref_press_gridded;

        fig3 = figure(3);
        ax = worldmap('Germany');
        C1=pcolorm(lati_mesh,loni_mesh,0.01.*press_perturbation_gridded);
        caxis([-2 2]);
        setm(ax,'MapLatLimit',[min(lati),max(lati)]);
        setm(ax,'MapLonLimit',[min(loni),max(loni)]);
        gridm('MeridianLabel','on', ...
            'ParallelLabel','on',...
            'glinestyle',':','gcolor','k','glinewidth',1,...
            'mlinelocation',9.7:0.1:10.4,'plinelocation',53.4:0.1:53.7)
        setm(ax,'FontSize',16,...
            'MLabelLocation',9.7:0.1:10.4,'MLabelRound',-1,...
            'PLabelLocation',53.4:0.1:53.7,'PLabelRound',-1)
        geoshow(germany_map_data,'FaceColor','none')
        str1 = char(string(datetime_reference_grid));
        str2 = char(string(datetime_plot_apollo_grid(t_ind)));
        title({strcat('Pressure perturbation (hPa) relative to',32,str1,32,'UTC');strcat('at',32,str2,32,'UTC')});
        colorbar EastOutside
        image_here = getframe(fig3);
        imwrite(image_here.cdata,...
            strcat('/Volumes/DATA3/FESSTHH2020data/images/press_pert_',...
            string(datetime_plot_apollo_grid(t_ind)),'.png'),'png');
        
        gifname = '/Volumes/DATA3/FESSTHH2020data/images/press_pert.gif';
        gif_delaytime = 0.75;
        [A,map] = rgb2ind(image_here.cdata,256);
        if t_ind == 1
            imwrite(A,map,gifname,'gif','LoopCount',Inf,'DelayTime',gif_delaytime);
        else
            imwrite(A,map,gifname,'gif','WriteMode','append','DelayTime',gif_delaytime);
        end
    end
end
% Plot WXT variables:
% max hail rate for color axis, in impacts/(m^2 second)
hail_wxt_max = 0.1.*max(hail_wxt(:),[],'omitnan');
% max precip rate for color axis, in mm/hour
% precip_wxt_max = 360.*max(precip_wxt(:),[],'omitnan');
precip_wxt_max = 50; 
% ignore this - used to assess best max and min values to use for plotting mixing ratio
rv_min = 1000;
rv_max = 0;
for t_ind = 1:num_indices_plot_grid
    
    t_ind_of_entire_day_wxt = indices_plot_wxt_grid(t_ind);
    
    if plot_wxt_hail
        hail_data_wxt = 0.1.*hail_wxt(:,t_ind_of_entire_day_wxt);
        indices_to_include_hail = ~isnan(hail_data_wxt);
        hail_interpolant = scatteredInterpolant(lat_all(indices_to_include_hail),...
            lon_all(indices_to_include_hail),hail_data_wxt(indices_to_include_hail),'linear','none');
        hail_gridded = hail_interpolant(lati_mesh,loni_mesh);
        
        fig4 = figure(4);
        ax = worldmap('Germany');
        C1=pcolorm(lati_mesh,loni_mesh,hail_gridded);
        caxis([0 hail_wxt_max]);
        setm(ax,'MapLatLimit',[min(lati),max(lati)]);
        setm(ax,'MapLonLimit',[min(loni),max(loni)]);
        gridm('MeridianLabel','on', ...
            'ParallelLabel','on',...
            'glinestyle',':','gcolor','k','glinewidth',1,...
            'mlinelocation',9.7:0.1:10.4,'plinelocation',53.4:0.1:53.7)
        setm(ax,'FontSize',16,...
            'MLabelLocation',9.7:0.1:10.4,'MLabelRound',-1,...
            'PLabelLocation',53.4:0.1:53.7,'PLabelRound',-1)
        geoshow(germany_map_data,'FaceColor','none')
        %str1 = char(string(datetime_reference_grid));
        str2 = char(string(datetime_plot_wxt_grid(t_ind)));
        title(strcat('Hail strikes (m^{-2} s^{-1}) at',32,str2,32,'UTC'));
        colorbar EastOutside
        image_here = getframe(fig4);
        imwrite(image_here.cdata,...
            strcat('/Volumes/DATA3/FESSTHH2020data/images/hail_',...
            string(datetime_plot_wxt_grid(t_ind)),'.png'),'png');
        gifname = '/Volumes/DATA3/FESSTHH2020data/images/hail.gif';
        gif_delaytime = 0.75;
        [A,map] = rgb2ind(image_here.cdata,256);
        if t_ind == 1
            imwrite(A,map,gifname,'gif','LoopCount',Inf,'DelayTime',gif_delaytime);
        else
            imwrite(A,map,gifname,'gif','WriteMode','append','DelayTime',gif_delaytime);
        end
    end
    if plot_wxt_precip
        precip_data_wxt = 360.*precip_wxt(:,t_ind_of_entire_day_wxt);
        indices_to_include_precip = ~isnan(precip_data_wxt);
        precip_interpolant = scatteredInterpolant(lat_all(indices_to_include_precip),...
            lon_all(indices_to_include_precip),precip_data_wxt(indices_to_include_precip),'linear','none');
        precip_gridded = precip_interpolant(lati_mesh,loni_mesh);
        
        fig5 = figure(5);
        ax = worldmap('Germany');
        C1=pcolorm(lati_mesh,loni_mesh,precip_gridded);
        caxis([0 precip_wxt_max]);
        setm(ax,'MapLatLimit',[min(lati),max(lati)]);
        setm(ax,'MapLonLimit',[min(loni),max(loni)]);
        gridm('MeridianLabel','on', ...
            'ParallelLabel','on',...
            'glinestyle',':','gcolor','k','glinewidth',1,...
            'mlinelocation',9.7:0.1:10.4,'plinelocation',53.4:0.1:53.7)
        setm(ax,'FontSize',16,...
            'MLabelLocation',9.7:0.1:10.4,'MLabelRound',-1,...
            'PLabelLocation',53.4:0.1:53.7,'PLabelRound',-1)
        geoshow(germany_map_data,'FaceColor','none')
        %str1 = char(string(datetime_reference_grid));
        str2 = char(string(datetime_plot_wxt_grid(t_ind)));
        title(strcat('Precip rate (mm h^{-1}) at',32,str2,32,'UTC'));
        colorbar EastOutside
        image_here = getframe(fig5);
        imwrite(image_here.cdata,...
            strcat('/Volumes/DATA3/FESSTHH2020data/images/precip_',...
            string(datetime_plot_wxt_grid(t_ind)),'.png'),'png');
        gifname = '/Volumes/DATA3/FESSTHH2020data/images/precip.gif';
        gif_delaytime = 0.75;
        [A,map] = rgb2ind(image_here.cdata,256);
        if t_ind == 1
            imwrite(A,map,gifname,'gif','LoopCount',Inf,'DelayTime',gif_delaytime);
        else
            imwrite(A,map,gifname,'gif','WriteMode','append','DelayTime',gif_delaytime);
        end
    end
    if plot_wxt_rv
        
        rh_data_wxt = hur_wxt(:,t_ind_of_entire_day_wxt);
        pa_data_wxt = pa_wxt(:,t_ind_of_entire_day_wxt);
        ta_data_wxt = ta_wxt(:,t_ind_of_entire_day_wxt);
        
        % Water vapor mixing ratio:
        rv_data_wxt = calc_mix_ratio(ta_data_wxt,rh_data_wxt,pa_data_wxt);
        
        indices_to_include_rv = ~isnan(rv_data_wxt);
        rv_interpolant = scatteredInterpolant(lat_all(indices_to_include_rv),...
            lon_all(indices_to_include_rv),rv_data_wxt(indices_to_include_rv),'linear','none');
        rv_gridded = rv_interpolant(lati_mesh,loni_mesh);
        
        % Uncomment these if needing to assess what bounds to use for color
        % axis:
        %
        % rv_min = min(min(rv_gridded(:)),rv_min);
        % rv_max = max(max(rv_gridded(:)),rv_max);
        
        fig6 = figure(6);
        ax = worldmap('Germany');
        C1=pcolorm(lati_mesh,loni_mesh,rv_gridded);
        caxis([25 45]);
        setm(ax,'MapLatLimit',[min(lati),max(lati)]);
        setm(ax,'MapLonLimit',[min(loni),max(loni)]);
        gridm('MeridianLabel','on', ...
            'ParallelLabel','on',...
            'glinestyle',':','gcolor','k','glinewidth',1,...
            'mlinelocation',9.7:0.1:10.4,'plinelocation',53.4:0.1:53.7)
        setm(ax,'FontSize',16,...
            'MLabelLocation',9.7:0.1:10.4,'MLabelRound',-1,...
            'PLabelLocation',53.4:0.1:53.7,'PLabelRound',-1)
        geoshow(germany_map_data,'FaceColor','none')
        %str1 = char(string(datetime_reference_grid));
        str2 = char(string(datetime_plot_wxt_grid(t_ind)));
        title(strcat('Water vapor mixing ratio (g kg^{-1}) at',32,str2,32,'UTC'));
        colorbar EastOutside
        image_here = getframe(fig6);
        imwrite(image_here.cdata,...
            strcat('/Volumes/DATA3/FESSTHH2020data/images/rv_',...
            string(datetime_plot_wxt_grid(t_ind)),'.png'),'png');
        gifname = '/Volumes/DATA3/FESSTHH2020data/images/rv.gif';
        gif_delaytime = 0.75;
        [A,map] = rgb2ind(image_here.cdata,256);
        if t_ind == 1
            imwrite(A,map,gifname,'gif','LoopCount',Inf,'DelayTime',gif_delaytime);
        else
            imwrite(A,map,gifname,'gif','WriteMode','append','DelayTime',gif_delaytime);
        end
    end
    
    if plot_wxt_winds
        wdir_data_wxt = wdir_wxt(:,t_ind_of_entire_day_wxt);
        wspeed_max_data_wxt = wspeed_max_wxt(:,t_ind_of_entire_day_wxt);
        wspeed_data_wxt = wspeed_wxt(:,t_ind_of_entire_day_wxt);
        % u-wind and v-wind:
        [u_data_wxt,v_data_wxt] = convert_wind(wspeed_data_wxt,wdir_data_wxt);
        % gust ratio:
        gust_ratio = wspeed_max_data_wxt./wspeed_data_wxt;
        gust_ratio_data_wxt = min(gust_ratio,10);

        indices_to_include_gust_ratio = ~isnan(gust_ratio_data_wxt);
        gust_ratio_interpolant = scatteredInterpolant(lat_all(indices_to_include_gust_ratio),...
            lon_all(indices_to_include_gust_ratio),gust_ratio_data_wxt(indices_to_include_gust_ratio),'linear','none');
        gust_ratio_gridded = gust_ratio_interpolant(lati_mesh,loni_mesh);

        indices_to_include_v = ~isnan(v_data_wxt);
        v_interpolant = scatteredInterpolant(lat_all(indices_to_include_v),...
            lon_all(indices_to_include_v),v_data_wxt(indices_to_include_v),'linear','none');
        v_gridded = v_interpolant(lati_mesh,loni_mesh);

        indices_to_include_u = ~isnan(u_data_wxt);
        u_interpolant = scatteredInterpolant(lat_all(indices_to_include_u),...
            lon_all(indices_to_include_u),u_data_wxt(indices_to_include_u),'linear','none');
        u_gridded = u_interpolant(lati_mesh,loni_mesh);

        fig7 = figure(7);
        ax = worldmap('Germany');
        C1=pcolorm(lati_mesh,loni_mesh,sqrt(u_gridded.^2 + v_gridded.^2));
        hold on
        C2=quiverm(lati_mesh,loni_mesh,v_gridded,u_gridded,'w');
        hold off
        caxis([0 8]);
        setm(ax,'MapLatLimit',[min(lati),max(lati)]);
        setm(ax,'MapLonLimit',[min(loni),max(loni)]);
        gridm('MeridianLabel','on', ...
            'ParallelLabel','on',...
            'glinestyle',':','gcolor','k','glinewidth',1,...
            'mlinelocation',9.7:0.1:10.4,'plinelocation',53.4:0.1:53.7)
        setm(ax,'FontSize',16,...
            'MLabelLocation',9.7:0.1:10.4,'MLabelRound',-1,...
            'PLabelLocation',53.4:0.1:53.7,'PLabelRound',-1)
        geoshow(germany_map_data,'FaceColor','none')
        %str1 = char(string(datetime_reference_grid));
        str2 = char(string(datetime_plot_wxt_grid(t_ind)));
        title(strcat('Wind speed (m s^{-1}) at',32,str2,32,'UTC'));
        colorbar EastOutside
        image_here = getframe(fig7);
        imwrite(image_here.cdata,...
            strcat('/Volumes/DATA3/FESSTHH2020data/images/wind_',...
            string(datetime_plot_wxt_grid(t_ind)),'.png'),'png');
        gifname = '/Volumes/DATA3/FESSTHH2020data/images/wind.gif';
        gif_delaytime = 0.75;
        [A,map] = rgb2ind(image_here.cdata,256);
        if t_ind == 1
            imwrite(A,map,gifname,'gif','LoopCount',Inf,'DelayTime',gif_delaytime);
        else
            imwrite(A,map,gifname,'gif','WriteMode','append','DelayTime',gif_delaytime);
        end

        fig8 = figure(8);
        ax = worldmap('Germany');
        C1=pcolorm(lati_mesh,loni_mesh,gust_ratio_gridded);
        caxis([1 10]);
        hold on
        C2=quiverm(lati_mesh,loni_mesh,v_gridded,u_gridded,'w');
        hold off
        setm(ax,'MapLatLimit',[min(lati),max(lati)]);
        setm(ax,'MapLonLimit',[min(loni),max(loni)]);
        gridm('MeridianLabel','on', ...
            'ParallelLabel','on',...
            'glinestyle',':','gcolor','k','glinewidth',1,...
            'mlinelocation',9.7:0.1:10.4,'plinelocation',53.4:0.1:53.7)
        setm(ax,'FontSize',16,...
            'MLabelLocation',9.7:0.1:10.4,'MLabelRound',-1,...
            'PLabelLocation',53.4:0.1:53.7,'PLabelRound',-1)
        geoshow(germany_map_data,'FaceColor','none')
        %str1 = char(string(datetime_reference_grid));
        str2 = char(string(datetime_plot_wxt_grid(t_ind)));
        title(strcat('Wind gust speed/wind speed at',32,str2,32,'UTC'));
        colorbar EastOutside
        image_here = getframe(fig8);
        imwrite(image_here.cdata,...
            strcat('/Volumes/DATA3/FESSTHH2020data/images/wind_gust_',...
            string(datetime_plot_wxt_grid(t_ind)),'.png'),'png');
        gifname = '/Volumes/DATA3/FESSTHH2020data/images/wind_gust.gif';
        gif_delaytime = 0.75;
        [A,map] = rgb2ind(image_here.cdata,256);
        if t_ind == 1
            imwrite(A,map,gifname,'gif','LoopCount',Inf,'DelayTime',gif_delaytime);
        else
            imwrite(A,map,gifname,'gif','WriteMode','append','DelayTime',gif_delaytime);
        end
    
    end
end

function rv = calc_mix_ratio(T,RH,p)
    % T is temperature in Kelvin
    % RH relative humidity as a fraction
    % p is pressure in Pa
    % rv is water vapor mixing ratio in g/kg
    esat_value = satvap_water_MK05(T);
    e_value = RH.*esat_value;
    epsilon_value = 461.5./287.047;
    rv = 1000.*epsilon_value.*e_value./(p - e_value);
end

function es = satvap_water_MK05(T)
    % Calculate saturation vapor pressure over liquid water, using an
    % equation from Murphy and Koop (2005, QJRMS)
    es = exp(54.842763 - 6763.22./T - 4.210.*log(T) + 0.000367.*T ...
        + tanh(0.0415.*(T - 218.8)).*(53.878 - 1331.22./T ...
        - 9.44523.*log(T) + 0.014025.*T));
end

function [u,v] = convert_wind(wspeed,wdir)
    math_direction = 270 - wdir;
    u = wspeed.*cosd(math_direction);
    v = wspeed.*sind(math_direction);
end