diff --git a/SYNOP/STATS/python/plot_p_values_map.py b/SYNOP/STATS/python/plot_p_values_map.py
index 7b6e3527bb32cf416e20330327ff36769ca5f2bb..4c17da7b5b1e6d0e9ad8e494230425161e861f2d 100644
--- a/SYNOP/STATS/python/plot_p_values_map.py
+++ b/SYNOP/STATS/python/plot_p_values_map.py
@@ -2,6 +2,8 @@
 
 '''
 Written by Madeleine Ekblom, 2023-04-06
+updates 2023-11-27
+- Option of moving legend into plot when plotting on a global map(lower left corner)
 
 Based on Jouni Räisänen's script for plotting p values on the map of Finland
 
@@ -38,43 +40,95 @@ def plot_p_values(p_values):
     Input: numpy array  
     Output: saves a figure in the running directory
     '''
-
-    lon_min=19
-    lon_max=32
-    lat_min=59
-    lat_max=71
-
+    # add more variables and units when needed
+    variable_shortnames={'39': '2t'}
+    variable_units={'39': '$^\circ$C'}
+    
+    # variables to add, shortnames and units need to be checked
+    # total amount of clouds, 91
+    # sea level pressure, 108 
+    # 1-hour precipitation, 80
+    # relative humidity, 58
+    # 10-minute precipitation intensity, 999
+    # snow depth, 71 
+    # 2m temp, 39 --> 2t, degrees Celcius
+    # dew point temperature, 40 
+    # visibility, 6
+    # wind direction, 111 
+    # max wind gust in last 10 mins,261
+    # surface wind speed, 221
+    
+    variable='39' # change to be an input value (as string)
+    
+    # if data only from SYNOP stations in Finland:
+    only_Finland_data=True
+    
+    # Define area of map:
+    if only_Finland_data:
+        # Finland:
+        lon_min=19
+        lon_max=32
+        lat_min=59
+        lat_max=71
+    else:
+        # Global:
+        lon_min=0
+        lon_max=359
+        lat_min=-90
+        lat_max=90
+
+    # define colours and limits for the colours:
     colors = ['darkviolet', 'blue', 'skyblue', 'lime', 'yellow', 'orange', 'red', 'grey']
     limits = [0,0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5,1]
 
-    fig = plt.figure(figsize=(12,8))
+    # create figure and add a subplot with a map:
+    fig = plt.figure(figsize=(8,12))
     ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree())
     ax.coastlines()
     ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())
     ax.gridlines(draw_labels=True)
     
-    # basic plot, not changing colors of p:       
-    #ax.scatter(p_values['lon'], p_values['lat'], c=p_values['p'])
-
+    # Plotting details:
+    fs1=20 # fontsize for title
+    fs2=14 # fontsize for legend
+    fs3=10 # fontisze for text box explanation
+    
     #loop over limits:
     legend_elements = []
     for n in range(8):
+        # find indices 
         p_ind = ((p_values['p'] > limits[n]) & (p_values['p'] <= limits[n+1]))
+        # scatter plot
         ax.scatter(p_values['lon'][p_ind], p_values['lat'][p_ind], c=colors[n])
+        # add to legend_elements, color and values
         legend_elements.append(Line2D([0], [0], marker='o', color='w',markerfacecolor=colors[n], label=str(limits[n])+'-'+str(limits[n+1])))
 
-    # legend table with customized values
-
-    #print(legend_elements)
-    
-    ax.legend(handles=legend_elements, loc='center left', bbox_to_anchor=(1.1, 0.5))
-    ax.set_title('P-values for T2m quantiles in ' + timestring, fontsize=20)
-
+    # explanation for figure:
+    text_expl ='p-value: the probability that the MSD statistics for the quantile rank\nhistogram would be exceeded by pure chance\nwhen comparing a perfect model simulation with observations.'
+
+    #  add legend table with customized values
+    if only_Finland_data:
+        # If Finland map, put legend outside map
+        ax.legend(handles=legend_elements, loc='center left', fontsize=fs3) #, bbox_to_anchor=(1.1, 0.5), fontsize=fs2)
+        # add text box with explanation to the bottom of the figure: 0.1, -0.25
+        ax.text(0.1, -0.25, text_expl, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3, transform=ax.transAxes)
+    else:
+        # If global map, put Legend in the lower left corner
+        ax.legend(handles=legend_elements, loc='lower left', fontsize=fs3)
+        # add text box with explanation to the bottom of the figure, 0.1, -0.45
+        ax.text(0.1, -0.45, text_expl, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3, transform=ax.transAxes)	
+    	
+    # set title:
+    fig.suptitle('SYNOP stations: p-values for ' + variable_shortnames[variable] +'\nquantiles in ' + timestring, fontsize=fs1)
+
+    # save figure
     plt.savefig('p_values.png', dpi=300)
-#    plt.show()
     
 
 def main(p_values_file):
+    '''
+    Main file: (1) reads in p values from file, (2) plot p values on a map
+    '''
     
     p_values_data = read_data(p_values_file)
 
diff --git a/SYNOP/STATS/python/plot_quantiles_rankhist.py b/SYNOP/STATS/python/plot_quantiles_rankhist.py
index b5fd9958ffa4fe85ca38d8a56ea5c8962e8ed7ef..aec03020e2c76c48bce952812d05a1fd0ac46ac8 100644
--- a/SYNOP/STATS/python/plot_quantiles_rankhist.py
+++ b/SYNOP/STATS/python/plot_quantiles_rankhist.py
@@ -2,6 +2,8 @@
 
 '''
 Written by Madeleine Ekblom, 2023-05-10
+updates 2023-11-27:
+- visual changes of the plot: labels, legends , title, and colours updated
 
 Based on Jouni Räisänen's script for plotting quantiles/time series data (00, 06, 12, 18 UTC) and rank histogram (24UTC=combination of 00,06,12,and 18 UTC data)
 
@@ -35,35 +37,109 @@ import sys
 import configparser
 import datetime
 from matplotlib.dates import drange, MonthLocator, DateFormatter
+import matplotlib.lines as mlines
+import matplotlib.patches as mpatches
+import matplotlib.ticker as ticker
 
 config = configparser.ConfigParser()
 
-#plot_data(quantiles_data, rank_hist_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end, True)
-def plot_data(quantiles_data, rank_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end, savefig=False):
+def plot_data(quantiles_data, rank_data, time_series_data, variable, station_id, lon, lat, year_beg, year_end, month_beg, month_end, figname, savefig=False):
     '''
-    quantiles data give the lines and the time series data give the dots
+    Plot data: quantiles data give the lines and the time series data give the dots
     '''
+    # add more variables and units when needed
+    variable_shortnames={'39': '2t'}
+    variable_units={'39': '$^\circ$C'}
+    
+    # variables to add, shortnames and units need to be checked
+    # total amount of clouds, 91
+    # sea level pressure, 108 
+    # 1-hour precipitation, 80
+    # relative humidity, 58
+    # 10-minute precipitation intensity, 999
+    # snow depth, 71 
+    # 2m temp, 39 --> 2t, degrees Celcius
+    # dew point temperature, 40 
+    # visibility, 6
+    # wind direction, 111 
+    # max wind gust in last 10 mins,261
+    # surface wind speed, 221
+    
+    # experiment id: now use 'test data' -  to be replaced with real exp_id (as input argument)
+    exp_id = 'test data' 
+        
     # set up date array consisting of all days in beg year:
     days = np.arange(1,366)
-#    day1 = datetime.date(year_beg, 1, 1)
-#    day2 = datetime.date(year_beg+1, 1, 1)
+    #day1 = datetime.date(year_beg, 1, 1)
+    #day2 = datetime.date(year_beg+1, 1, 1)
     day1 = datetime.date(2001, 1, 1)
     day2 = datetime.date(2002, 1, 1)
     delta = datetime.timedelta(days=1)
     dates = drange(day1, day2, delta)
-
+    
     # arrays for quantile names, axes titles, hours to plot
     q_names = ['q01', 'q10', 'q25', 'q50', 'q75', 'q90', 'q99']
     sub_titles = ['00 UTC', '06 UTC', '12 UTC', '18 UTC']
     hours = [0, 6, 12, 18]
-
+    
+    # plotting details:
+    alpha = [0.5, 0.3, 0.1] # transparency for shaded area plotting
+    color='blue' # color for shading
+    linecolor='red' # color for solid line
+    scatter_size=5 # size for scattered dots
+    fs1=20 # fontsize for title
+    fs2=14 # fontsize for labels, subplot titles
+    fs3 = 12 # fontsize for
+    fs4=10 # fontsize for extra text box, legend
+    ndigits = 3 # round off to ndigits (e.g., lon, lat)
+    
+    # use variable to get variable and corresponding unit, used as ylabel in plots:
+    ylabel = variable_shortnames[variable] + ' (' + variable_units[variable] + ')'
+        
+    # legend elemnents for quantile data, one legend for model and one legend for observations
+    legend_elements1 = [mlines.Line2D([],[], marker='.', linestyle='', color='k', label='Model data (' + exp_id + ')')]
+    legend_elements2= [mpatches.Patch(color=color, alpha=alpha[0], label='25-50%'),
+                                               mpatches.Patch(color=color, alpha=alpha[1], label='10-25%'),
+                                               mpatches.Patch(color=color, alpha=alpha[2], label='1-10%'),
+                                               mlines.Line2D([],[], color=linecolor, label='50%'),
+                                               mpatches.Patch(color=color, alpha=alpha[0], label='50-75%'),
+                                               mpatches.Patch(color=color, alpha=alpha[1], label='75-90%'),
+                                               mpatches.Patch(color=color, alpha=alpha[2], label='90-99%')]
+    
+    # round lon, lat to three digits:
+    lon = np.round(float(lon), ndigits)
+    lat = np.round(float(lat), ndigits)
+    
+    # define lon/lat text strings for title of the form lat: XX N/S, lat: XX E/W
+    if (lat > 0):
+        # if lat positive -> N
+        lat_text = ' (lat: ' + str(lat) + '$^\circ$N, '
+    else:
+        # if negative -> S
+        lat_text = ' (lat ' + str(-lat) + '$^\circ$S, '
+    
+    if (0 <= lon) and (lon <= 180):
+        # lon btw 0 an 180 -> E
+        lon_text = ' lon: ' + str(lon) + '$^\circ$E)'
+    else:
+        if lon < 0:
+            # lon < 0 --> - lon W
+            lon_text = ' lon: ' + str(-lon) + '$^\circ$W)'
+        else:
+            # 180 < lon < 360 --> - (lon - 360) W
+            lon = -(lon - 360)
+            lon_text = ' lon: ' + str(lon) + '$^\circ$W)'
+    
+    # define figure title:
+    fig_title = 'SYNOP station: ' + station_id + lat_text + lon_text + '\nModel data from ' + str(year_beg) + month_beg + '-' + str(year_end) + month_end 
+    
     # calculate number of years:
     nyears = year_end - year_beg + 1
 
     # set up figure
     fig = plt.figure(figsize=(10,10), layout='constrained')
-    spec = fig.add_gridspec(3,2)
-    fig.suptitle(fig_title, fontsize=20)
+    spec = fig.add_gridspec(4,2)
+    fig.suptitle(fig_title, fontsize=fs1)
 
     # counter
     c = 0
@@ -71,40 +147,95 @@ def plot_data(quantiles_data, rank_data, time_series_data, fig_title, fig_name,
     # plot quantiles/time series for times 00, 06, 12, 18:
     for i in range(2):
         for j in range(2):
-            # set up axis: title, xaxis
+            # set up axis: title, xaxis, x- and y-labels
             ax = fig.add_subplot(spec[i,j])
-            ax.set_title(sub_titles[c])
+            ax.set_title(sub_titles[c], fontsize=fs2)
             ax.set_xlim(dates[0], dates[-1])
             ax.xaxis.set_major_locator(MonthLocator())
             ax.xaxis.set_major_formatter(DateFormatter('%b'))
-    
+            ax.set_xlabel('Time', fontsize=fs2)
+            ax.set_ylabel(ylabel, fontsize=fs2)
+            
             # quantile data:
             # find quantile data hour = 0, 6, 12, 18
             qh_ind = (quantiles_data['hour'][:] == hours[c])
-            for q in q_names:
-                q_data = quantiles_data[q][qh_ind]
-                ax.plot(dates, q_data[:365], 'k-') # --> change q only contains one year of data!!
+            
+            # get quantile data, qXX_data = quantile XX data at hour qh_ind
+            q01_data = quantiles_data[q_names[0]][qh_ind]
+            q10_data = quantiles_data[q_names[1]][qh_ind]
+            q25_data = quantiles_data[q_names[2]][qh_ind]
+            q50_data = quantiles_data[q_names[3]][qh_ind]
+            q75_data = quantiles_data[q_names[4]][qh_ind]
+            q90_data = quantiles_data[q_names[5]][qh_ind]
+            q99_data = quantiles_data[q_names[6]][qh_ind]
+            
+            # solid line at q50:
+            ax.plot(dates, q50_data, c=linecolor, ls='-')
+            
+            # shading will be the most opaque close to 50 and more transparent close to 1/99
+            
+            # fill area btw q25-q50, q50-q75: alpha = alpha[0] 
+            ax.fill_between(dates, q25_data[:], q50_data[:], color=color, alpha=alpha[0])
+            ax.fill_between(dates, q50_data[:], q75_data[:], color=color, alpha=alpha[0])
+            
+            # fill area btw q10-q25, q75-q90, alpha=alpha[1]
+            ax.fill_between(dates, q10_data[:], q25_data[:], color=color, alpha=alpha[1])
+            ax.fill_between(dates, q75_data[:], q90_data[:], color=color, alpha=alpha[1])
+
+            # fill area btw q01-q10, q90-q99, alpha=alpha[2]            
+            ax.fill_between(dates, q01_data[:], q10_data[:], color=color, alpha=alpha[2])
+            ax.fill_between(dates, q90_data[:], q99_data[:], color=color, alpha=alpha[2])            
+            
             # plot time series data:
             # find time series where hour = 0, 6, 12, 18
             th_ind = (time_series_data['hour'][:] == hours[c])
             t_data = time_series_data['value'][th_ind] # all years
+            
             for n in range(nyears):
-                ax.scatter(dates, t_data[n*365:(n+1)*365], marker='.')
+                # plot data with black, small dots
+                ax.scatter(dates, t_data[n*365:(n+1)*365], s=scatter_size, c='k', marker='.') 
             c = c + 1
 
     # plot rank histogram data:
-    ax2 = fig.add_subplot(spec[2,:])
-    ax2.set_title('Rank histogram\n' + 'MSD:' + str(rank_data[4,4]) + ' p-value:' + str(rank_data[4,5]))
-    ax2.bar(np.arange(100), rank_data[4,6:]*100)
+    ax2 = fig.add_subplot(spec[2,0])
+    ax2.set_title('Rank histogram', fontsize=fs2)
+    ax2.bar(np.arange(100)+0.5, rank_data[4,6:]*100) # add + 0.5 to center bars
     ax2.set_xlim(0,100)
-    ax2.axhline(y=1)
+    ax2.set_xlabel('Percentile', fontsize=fs2)
+    ax2.set_ylabel('Relative frequency ', fontsize=fs2)
+    ax2.axhline(y=1, c='k', ls='--') # horizontal line at y=1.0 
+    ax2.xaxis.set_major_locator(ticker.FixedLocator([10,25,50,75,90]))
+    
+    # add a sixth axis for legends and text data:
+    ax3 = fig.add_subplot(spec[2,1])
+    ax3.axis('off')
+        
+    # legend with two columns at the upper centre of the figure:                                   
+    leg1 = ax3.legend(loc='upper center', handles=legend_elements1, fontsize=fs4)
+    leg2 = ax3.legend(loc='lower center', handles=legend_elements2, ncols=2, fontsize=fs4, title='Quantiles for observations', title_fontsize=fs4)
+    ax3.add_artist(leg1)
+    ax3.add_artist(leg2)
+
+    # add text boxes with (1) rank histogram data and (2) explanation
+    rank_hist_text = r'$\bf{MSD}$:' + str(np.round(rank_data[4,4], ndigits)) + '\n' + r'$\bf{ p-value}$:' + str(np.round(rank_data[4,5], ndigits))
+   
+    # add text boxes to axis in right lower corner
+    ax3.text(0.35, -0.3, rank_hist_text, bbox={'pad': 10, 'facecolor':'white'}, fontsize=fs3)
     
+    # figure explanation: explaining what quantiles and model data are
+    expl_msd_pvalue = 'Mean Square Deviation (' + r'$\bf{MSD}$' + ' relative to flat histogram\n' + r'$\bf{p-value}$' + ' of MSD relative to sampling distribution\n'
+    expl_quantile='Quantiles are estimated based on observations at the station ' + station_id +' (period: 2010-2021)\n' 
+    expl_model='Model data from experiment (' + exp_id + ') interpolated to location of this station '
+    fig_text = expl_msd_pvalue + expl_quantile + expl_model
+    
+    # add subplot at the bottom of the figure, and add figure explanation as text to this subplot
+    ax4 = fig.add_subplot(spec[3,:])
+    ax4.axis('off')
+    ax4.text(0.1, 0.4, fig_text, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs4)
+    
+    # save figure:
     if savefig:
-        plt.savefig(fig_name, dpi=300)
-
-    #plt.savefig(fig_name, dpi=300)
-
-#    plt.show( )
+        plt.savefig(figname, dpi=300)
 
 
 def read_quantiles_data(quantile_sel_file):
@@ -161,8 +292,8 @@ def main(config_file):
     # data to be plotted
     variable = config['data']['variable']
     station_id = config['data']['station_id']
-    longitude = config['data']['longitude']
-    latitude = config['data']['latitude']
+    lon = config['data']['longitude']
+    lat = config['data']['latitude']
     year_beg = int(config['data']['year_beg'])
     year_end = int(config['data']['year_end'])
     month_beg = config['data']['month_beg']
@@ -174,7 +305,7 @@ def main(config_file):
     rank_hist_file = config['input']['rank_hists']
     
     # output files
-    fig_name = config['output']['fig_name']
+    figname = config['output']['fig_name']
     
     ## Read input data ##
     time_series_data = read_time_series_data(sim_file)
@@ -182,8 +313,7 @@ def main(config_file):
     rank_hist_data = read_rank_data(rank_hist_file)
 
     ## Plot data ##
-    fig_title = 'Station: ' + station_id + ', ' + str(year_beg) + month_beg + '-' + str(year_end) + month_end + '\nLat=' + latitude + ' Lon=' + longitude
-    plot_data(quantiles_data, rank_hist_data, time_series_data, fig_title, fig_name, year_beg, year_end, month_beg, month_end, True)
+    plot_data(quantiles_data, rank_hist_data, time_series_data, variable, station_id, lon, lat, year_beg, year_end, month_beg, month_end, figname, True)
     
     
 if __name__=='__main__':
diff --git a/SYNOP/STATS/python/plot_rank_hist_sum_all_stations.py b/SYNOP/STATS/python/plot_rank_hist_sum_all_stations.py
index 602caac02c2942ceacec6554b58eac26fb107971..10ebd3235cdfecdd53cf1eebf1339d3a8f5d9b5e 100644
--- a/SYNOP/STATS/python/plot_rank_hist_sum_all_stations.py
+++ b/SYNOP/STATS/python/plot_rank_hist_sum_all_stations.py
@@ -2,6 +2,8 @@
 
 '''
 Written by Madeleine Ekblom, 2023-09-12
+updates 2023-11-27:
+- moved text from title to a new text box
 
 Based on Jouni Räisänen's script for plotting rank histogram summart statistics 
 
@@ -15,27 +17,57 @@ import matplotlib.pyplot as plt
 from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)
 import sys
 
-#def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, figname):
 def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q):
     '''
-    Input: p_freq, q_freq, number of stations, p01, p1, p5, max freq p, max freq q, figname 
+    Input: p_freq, q_freq, number of stations, p01, p1, p5, max freq p, max freq q
     Description: Plots p and q frequencies and saves the figure with figname
     '''
-    title='Rank Histogram Summary Statistics\n Number of stations: ' + str(number_of_stations) + '\nFrequency, p<0.001: ' + str(p01) + '\nFrequency, p<0.01:' + str(p1) + '\nFrequency, p<0.05:' + str(p5)
+    # add more variables and units when needed
+    variable_shortnames={'39': '2t'}
+    variable_units={'39': '$^\circ$C'}
     
-    # can also be input arguments if the number of bins are not constant
+    # variables to add, shortnames and units need to be checked
+    # total amount of clouds, 91
+    # sea level pressure, 108 
+    # 1-hour precipitation, 80
+    # relative humidity, 58
+    # 10-minute precipitation intensity, 999
+    # snow depth, 71 
+    # 2m temp, 39 --> 2t, degrees Celcius
+    # dew point temperature, 40 
+    # visibility, 6
+    # wind direction, 111 
+    # max wind gust in last 10 mins,261
+    # surface wind speed, 221
+    
+    variable='39' # change to be an input value, as string 
+    
+    # title
+    title='SYNOP: Rank histogram summary statistics for '  + variable_shortnames[variable]
+    
+    # text box information
+    text_box = 'Number of stations: ' + str(number_of_stations) + '\nFrequency, p<0.001: ' + str(p01) + '\nFrequency, p<0.01: ' + str(p1) + '\nFrequency, p<0.05: ' + str(p5)
+
+    # these can also be input arguments if the number of bins are not constant
     number_of_p_bins=20
     number_of_q_bins=100
 
+    # set up figure
     fig = plt.figure(figsize=(8,12))
-    gs = fig.add_gridspec(2, hspace=0.4)
-    ax1, ax2 = gs.subplots(sharex=False, sharey=False)
+    gs = fig.add_gridspec(3, hspace=0.5)
+    ax1, ax2, ax3 = gs.subplots(sharex=False, sharey=False)
+    
+    # plotting details
+    fs1=20 # fontsize for title
+    fs2=14 # fontsize for labels, subplot titles, textbox
+    fs3=10 # fontsize for text explanation
 
     # Plot p values
-    ax1.set_title('Normalized p-value frequencies')
-    ax1.set_xlabel('p-value')
+    ax1.set_title('Normalized p-value frequencies', fontsize=fs2)
+    ax1.set_xlabel('p-value', fontsize=fs2)
+    ax1.set_ylabel('Relative frequency', fontsize=fs2)
     ax1.set_ylim([0,max_freq_p])
-    ax1.bar(np.arange(number_of_p_bins)+0.5, number_of_p_bins*p_freq)
+    ax1.bar(np.arange(number_of_p_bins)+0.5, number_of_p_bins*p_freq) # add +0.5 to center the bars
     ax1.xaxis.set_major_locator(MultipleLocator(2))
     ax1.xaxis.set_minor_locator(MultipleLocator(1))
     p_xlabels = ax1.get_xticks()
@@ -43,15 +75,34 @@ def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p,
     ax1.set_xlim([0,number_of_p_bins])
 
     # Plot q values
-    ax2.set_title('Normalized quantile frequencies')
-    ax2.set_xlabel('Quantile (%)')
-    ax2.bar(np.arange(number_of_q_bins)+0.5, number_of_q_bins*q_freq)
+    ax2.set_title('Normalized quantile frequencies', fontsize=fs2)
+    ax2.set_xlabel('Quantile (%)', fontsize=fs2)
+    ax2.set_ylabel('Relative frequency', fontsize=fs2)
+    ax2.bar(np.arange(number_of_q_bins)+0.5, number_of_q_bins*q_freq) # add +0.5 to center the bars
     ax2.xaxis.set_major_locator(MultipleLocator(10))
     ax2.xaxis.set_minor_locator(MultipleLocator(5))
     ax2.set_ylim([0,max_freq_q])
     ax2.set_xlim([0,number_of_q_bins])
 
-    fig.suptitle(title)
+    # remove axis from ax3
+    ax3.axis('off')
+    
+    # add text box with additional infromation:
+    ax3.text(0.3, 0.6, text_box, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs2)
+    
+    # explanation for figure
+    text_pvalue= 'Normalized p-values: Frequency distribution of p-values at\nall individual stations, normalized to give an expected value of 1.\n'
+    text_quantile='Normalized quantiles: Frequency distribution of model data\nrelative to quantiles of observations, averaged over all stations\nand divided by the expected value of 0.01.'
+    
+    text_expl = text_pvalue + text_quantile
+    
+    # add text box with explanation to the bottom of the figure
+    ax3.text(0.1, -0.05, text_expl, bbox={'pad':10, 'facecolor':'white'}, fontsize=fs3)
+
+    # add title to figure
+    fig.suptitle(title, fontsize=fs1)
+    
+    # save figure
     figname='rank_hist_sumstats.png'
     plt.savefig(figname, dpi=300)
 
@@ -80,9 +131,7 @@ if __name__=='__main__':
     p5 = sys.argv[5]
     max_freq_p = float(sys.argv[6])
     max_freq_q = float(sys.argv[7])
-    #figname = sys.argv[8]
 
     p_freq, q_freq = read_data(rh_summary_file)
-    #plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, figname)
     plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q)
 
diff --git a/SYNOP/STATS/python/plot_rank_hist_sum_all_stations_v2.py b/SYNOP/STATS/python/plot_rank_hist_sum_all_stations_v2.py
deleted file mode 100644
index e31630ae90ef0bc4fa5b4ea28a76d09f8bbca899..0000000000000000000000000000000000000000
--- a/SYNOP/STATS/python/plot_rank_hist_sum_all_stations_v2.py
+++ /dev/null
@@ -1,85 +0,0 @@
-#!/usr/bin/env python
-
-'''
-Written by Madeleine Ekblom, 2023-09-12
-
-Based on Jouni Räisänen's script for plotting rank histogram summart statistics 
-
-Example:
-$ python3 plot_rank_hist_sum_all_stations.py rh_summary_file nstat p01 p1 p5 max_freq_p max_freq_q
-
-'''
-
-import numpy as np
-import matplotlib.pyplot as plt
-from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)
-import sys
-
-def plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, figname):
-    '''
-    Input: p_freq, q_freq, number of stations, p01, p1, p5, max freq p, max freq q, figname 
-    Description: Plots p and q frequencies and saves the figure with figname
-    '''
-    title='Rank histogram summary statistics\n Number of stations: ' + str(number_of_stations) + '\nFrequency, p<0.001: ' + str(p01) + '\nFrequency, p<0.01: ' + str(p1) + '\nFrequency, p<0.05: ' + str(p5)
-
-    # can also be input arguments if the number of bins are not constant
-    number_of_p_bins=20
-    number_of_q_bins=100
-
-    fig = plt.figure(figsize=(8,12))
-    gs = fig.add_gridspec(2, hspace=0.4)
-    ax1, ax2 = gs.subplots(sharex=False, sharey=False)
-
-    # Plot p values
-    ax1.set_title('Normalized p-value frequencies')
-    ax1.set_xlabel('p-value')
-    ax1.set_ylim([0,max_freq_p])
-    ax1.bar(np.arange(number_of_p_bins)+0.5, number_of_p_bins*p_freq)
-    ax1.xaxis.set_major_locator(MultipleLocator(2))
-    ax1.xaxis.set_minor_locator(MultipleLocator(1))
-    p_xlabels = ax1.get_xticks()
-    ax1.set_xticks(p_xlabels, p_xlabels/number_of_p_bins)
-    ax1.set_xlim([0,number_of_p_bins])
-
-    # Plot q values
-    ax2.set_title('Normalized quantile frequencies')
-    ax2.set_xlabel('Quantile (%)')
-    ax2.bar(np.arange(number_of_q_bins)+0.5, number_of_q_bins*q_freq)
-    ax2.xaxis.set_major_locator(MultipleLocator(10))
-    ax2.xaxis.set_minor_locator(MultipleLocator(5))
-    ax2.set_ylim([0,max_freq_q])
-    ax2.set_xlim([0,number_of_q_bins])
-
-    fig.suptitle(title)
-    plt.savefig(figname)
-
-def read_data(rh_summary_file):
-    '''
-    Input: path/to/rh_summary_file
-    Output: p_freq, q_freq
-    Description: reads in p frequencies and q frequencies
-    '''
-
-    # first line contains a header and first column contains hour
-    p_freq_all = np.loadtxt(rh_summary_file + '_p-freq.txt', skiprows=1)
-    q_freq_all = np.loadtxt(rh_summary_file + '_q-freq.txt', skiprows=1)
-
-    # p_freq and q_freq contain data 00, 06, 12, 18, 24 (the last is used in these plots)
-    p_freq = p_freq_all[-1,1:]
-    q_freq = q_freq_all[-1,1:]
-
-    return p_freq, q_freq
-
-if __name__=='__main__':
-    rh_summary_file = sys.argv[1]
-    number_of_stations = sys.argv[2]
-    p01 = sys.argv[3]
-    p1 = sys.argv[4]
-    p5 = sys.argv[5]
-    max_freq_p = float(sys.argv[6])
-    max_freq_q = float(sys.argv[7])
-    figname = sys.argv[8]
-
-    p_freq, q_freq = read_data(rh_summary_file)
-    plot_rank_hist(p_freq, q_freq, number_of_stations, p01, p1, p5, max_freq_p, max_freq_q, figname)
-