diff --git a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
index c418dda2579c25ffb6c70f8d0f15fe8be0b5f91e..4d95d941c08143caa3e0a54b35ee4617f6fe87a7 100644
--- a/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
+++ b/extreme_data/meteo_france_data/scm_models_data/abstract_study.py
@@ -650,6 +650,10 @@ class AbstractStudy(object):
 
     """  Spatial properties """
 
+    @cached_property
+    def massif_name_to_massif_id(self):
+        return {name: i for i, name in enumerate(self.study_massif_names)}
+
     @classproperty
     def original_safran_massif_id_to_massif_name(cls) -> Dict[int, str]:
         return {massif_id: massif_name for massif_id, massif_name in enumerate(cls.all_massif_names)}
diff --git a/projects/contrasting_trends_in_snow_loads/altitudes_fit/altitudes_studies.py b/projects/contrasting_trends_in_snow_loads/altitudes_fit/altitudes_studies.py
index 84a6c8f640b19e134850a46dcd466b36319f8334..e0772e2eccbf3d86493308ecd81eacc2dab4c695 100644
--- a/projects/contrasting_trends_in_snow_loads/altitudes_fit/altitudes_studies.py
+++ b/projects/contrasting_trends_in_snow_loads/altitudes_fit/altitudes_studies.py
@@ -113,7 +113,7 @@ class AltitudesStudies(object):
         ax.xaxis.set_ticks(x[1::10])
         ax.tick_params(axis='both', which='major', labelsize=13)
         ax.legend()
-        plot_name = 'Annual maxima of {}\nin {}'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[self.study_class], massif_name)
+        plot_name = 'Annual maxima of {} in {}'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[self.study_class], massif_name)
         ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
         ax.set_xlabel('years', fontsize=15)
         self.show_or_save_to_file(plot_name=plot_name, show=show)
@@ -152,6 +152,7 @@ class AltitudesStudies(object):
                     moment = after - before
                     if change is None:
                         moment /= before
+                        moment *= 100
                 else:
                     moment = function(annual_maxima)
                 mean_moment.append(moment)
diff --git a/projects/contrasting_trends_in_snow_loads/altitudes_fit/main_altitudes_studies.py b/projects/contrasting_trends_in_snow_loads/altitudes_fit/main_altitudes_studies.py
index 759b35c8a02b6fa236afff734adcda746389f19a..08566ff0c7978a1d9da5d6d38452dee3d0ba0940 100644
--- a/projects/contrasting_trends_in_snow_loads/altitudes_fit/main_altitudes_studies.py
+++ b/projects/contrasting_trends_in_snow_loads/altitudes_fit/main_altitudes_studies.py
@@ -9,17 +9,17 @@ def main_plots_moments():
     study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
     study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
                      SafranPrecipitation7Days][:]
-    study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day]
+    study_classes = [SafranSnowfall1Day, SafranPrecipitation1Day][:]
 
     for study_class in study_classes:
 
         studies = AltitudesStudies(study_class, altitudes)
         # massifs_names = ['Vercors', 'Chartreuse', 'Belledonne']
         # studies.plot_mean_maxima_against_altitude(massif_names=massifs_names, show=True)
-        # studies.plot_maxima_time_series()
-        for std in [True, False]:
-            for change in [True, False, None]:
-                studies.plot_mean_maxima_against_altitude(std=std, change=change)
+        studies.plot_maxima_time_series()
+        # for std in [True, False]:
+        #     for change in [True, False, None]:
+        #         studies.plot_mean_maxima_against_altitude(std=std, change=change)
 
 
 if __name__ == '__main__':
diff --git a/projects/contrasting_trends_in_snow_loads/gorman_figures/figure  3/__init__.py b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure  3/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/projects/contrasting_trends_in_snow_loads/gorman_figures/figure3_daily snowfall fraction.py b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure  3/figure3_daily snowfall fraction.py
similarity index 79%
rename from projects/contrasting_trends_in_snow_loads/gorman_figures/figure3_daily snowfall fraction.py
rename to projects/contrasting_trends_in_snow_loads/gorman_figures/figure  3/figure3_daily snowfall fraction.py
index fd72bd0a2495b7fcbd7b61f87e41e2d12c40128c..672b84086c380a951f1825f70b3e87f5825960d0 100644
--- a/projects/contrasting_trends_in_snow_loads/gorman_figures/figure3_daily snowfall fraction.py	
+++ b/projects/contrasting_trends_in_snow_loads/gorman_figures/figure  3/figure3_daily snowfall fraction.py	
@@ -8,10 +8,11 @@ from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranR
 Altitudes_groups = [[300, 600, 900], [1200, 1500, 1800], [2100, 2400, 2700], [3000, 3300, 3600]]
 
 
-def compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim):
-    all_time_series = np.concatenate([get_time_series(altitude, year_min, year_max) for altitude in altitudes], axis=1)
+def compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim, massif_names=None):
+    all_time_series = np.concatenate([get_time_series(altitude, year_min, year_max, massif_names) for altitude in altitudes], axis=1)
     temperature_array, snowfall_array, rainfall_array = all_time_series
-    nb_bins_for_one_celsius = 4
+    # nb_bins_for_one_celsius = 4
+    nb_bins_for_one_celsius = 2
     bins = np.linspace(-lim, lim, num=int((nb_bins_for_one_celsius * 2 * lim) + 1))
     inds = np.digitize(temperature_array, bins)
     snowfall_fractions = []
@@ -23,9 +24,10 @@ def compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim):
     return np.array(snowfall_fractions), x
 
 
-def daily_snowfall_fraction(ax, altitudes, year_min=1959, year_max=2019):
+def daily_snowfall_fraction(ax, altitudes, year_min=1959, year_max=2019, massif_names=None):
     lim = 6
-    snowfall_fractions, x = compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim=lim)
+    snowfall_fractions, x = compute_smoother_snowfall_fraction(altitudes, year_min, year_max, lim=lim,
+                                                               massif_names=massif_names)
     sigma = 1.0
     snowfall_fractions = 100 * gaussian_filter1d(snowfall_fractions, sigma=sigma)
     # Plot results
@@ -62,12 +64,20 @@ def end_plot(ax, sigma):
     ax.legend()
 
 
-def get_time_series(altitude, year_min, year_max):
+def get_time_series(altitude, year_min, year_max, massif_names=None):
     temperature_study = SafranTemperature(altitude=altitude, year_min=year_min, year_max=year_max)
     study_rainfall = SafranRainfall1Day(altitude=altitude, year_min=year_min, year_max=year_max)
     study_snowfall = SafranSnowfall1Day(altitude=altitude, year_min=year_min, year_max=year_max)
-    all_time_series = [temperature_study.all_daily_series, study_snowfall.all_daily_series,
+    studies = [temperature_study, study_rainfall, study_snowfall]
+    if massif_names is None:
+        all
+    all_time_series = [temperature_study.all_daily_series,
+                       study_snowfall.all_daily_series,
                        study_rainfall.all_daily_series]
+    if massif_names is not None:
+        # Select only the index corresponding to the massif of interest
+        massif_ids = [temperature_study.massif_name_to_massif_id[name] for name in massif_names]
+        all_time_series = [a[:, massif_ids] for a in all_time_series]
     all_time_series = [np.concatenate(t) for t in all_time_series]
     all_time_series = np.array(all_time_series)
     return all_time_series
@@ -76,8 +86,11 @@ def get_time_series(altitude, year_min, year_max):
 def daily_snowfall_fraction_wrt_altitude(fast=True):
     ax = plt.gca()
     groups = [[1800]] if fast else Altitudes_groups
+    groups = [[a] for a in [900, 1200, 1500]]
+    massif_names = None
+    massif_names = ['Vercors']
     for altitudes in groups:
-        daily_snowfall_fraction(ax, altitudes)
+        daily_snowfall_fraction(ax, altitudes, massif_names=massif_names, year_max=1999)
     plt.show()