diff --git a/extreme_data/meteo_france_data/adamont_data/abstract_simulation_study.py b/extreme_data/meteo_france_data/adamont_data/abstract_simulation_study.py
index add2f79000d96097c23283c22a721318b3563af5..b852d342ce54b08feb35999574db854b3876d177 100644
--- a/extreme_data/meteo_france_data/adamont_data/abstract_simulation_study.py
+++ b/extreme_data/meteo_france_data/adamont_data/abstract_simulation_study.py
@@ -29,6 +29,10 @@ class SimulationStudy(AbstractStudy):
                  scenario="HISTO", ensemble_idx=0):
         super().__init__(variable_class, altitude, year_min, year_max, multiprocessing, orientation, slope, season,
                          french_region, split_years)
+        if scenario == 'HISTO':
+            self.year_min, self.year_max = 1950, 2004
+        else:
+            self.year_min, self.year_max = 2006, 2100
         assert scenario in self.scenarios
         assert 0 <= ensemble_idx <= 13
         self.scenario = scenario
@@ -39,22 +43,12 @@ class SimulationStudy(AbstractStudy):
 
     @property
     def simulations_path(self):
-        return op.join(ADAMONT_PATH, self.parameter, self.scenario)
-
-    @property
-    def parameter(self):
-        return self.variable_class_to_parameter[self.variable_class]
-
-    @classproperty
-    def variable_class_to_parameter(cls):
-        return {
-            SafranSnowfallSimulationVariable: 'Snow',
-            CrocusTotalSweVariable: 'SNOWSWE',
-        }
+        return op.join(ADAMONT_PATH, self.variable_class.keyword(), self.scenario)
 
     @property
     def nc_path(self):
         nc_file = os.listdir(self.simulations_path)[self.ensemble_idx]
+        print(" ", nc_file)
         return op.join(self.simulations_path, nc_file)
 
     @property
@@ -108,7 +102,7 @@ class SimulationStudy(AbstractStudy):
         year_to_data_list = {}
         for year in self.ordered_years:
             year_to_data_list[year] = []
-        data_list = self.dataset.variables[self.variable_class.keyword]
+        data_list = self.dataset.variables[self.variable_class.keyword()]
         data_year_list = self.winter_year_for_each_time_step
         assert len(data_list) == len(data_year_list)
         for year_data, data in zip(data_year_list, data_list):
@@ -124,7 +118,7 @@ class SimulationStudy(AbstractStudy):
         zs_list = [int(e) for e in np.array(self.dataset.variables['ZS'])]
         return np.array(zs_list) == self.altitude
 
-    @property
+    @cached_property
     def study_massif_names(self) -> List[str]:
         massif_ids = np.array(self.dataset.variables['MASSIF_NUMBER'])[self.flat_mask]
         return [self.massif_number_to_massif_name[massif_id] for massif_id in massif_ids]
diff --git a/extreme_data/meteo_france_data/adamont_data/snowfall_simulation.py b/extreme_data/meteo_france_data/adamont_data/snowfall_simulation.py
index 44e0d67b7bfcf2665be974fbfa345ca47c754e2b..e9d990ae26bfa84e7e3aafa8c544d9f4f130f838 100644
--- a/extreme_data/meteo_france_data/adamont_data/snowfall_simulation.py
+++ b/extreme_data/meteo_france_data/adamont_data/snowfall_simulation.py
@@ -7,6 +7,8 @@ from extreme_data.meteo_france_data.scm_models_data.utils import Season, FrenchR
 
 
 class SafranSnowfallSimulationVariable(AbstractVariable):
+    UNIT = 'kg $m^{-2}$'
+
 
     @property
     def daily_time_serie_array(self) -> np.ndarray:
@@ -21,10 +23,10 @@ class SafranSnowfallSimulationRCP85(SimulationStudy):
 
     def __init__(self, altitude: int = 1800, year_min=YEAR_MIN, year_max=YEAR_MAX,
                  multiprocessing=True, orientation=None, slope=20.0, season=Season.annual,
-                 french_region=FrenchRegion.alps, split_years=None):
+                 french_region=FrenchRegion.alps, split_years=None, ensemble_idx=0):
         super().__init__(SafranSnowfallSimulationVariable, altitude, year_min, year_max, multiprocessing, orientation,
                          slope,
-                         season, french_region, split_years, "RCP85")
+                         season, french_region, split_years, "RCP85", ensemble_idx)
 
 
 if __name__ == '__main__':
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 35d40760169157cc521560192e8eaa13a6d161d1..b48d95891efaa342096f4b990f627b9bd3cdebc1 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
@@ -24,7 +24,8 @@ from extreme_data.meteo_france_data.scm_models_data.abstract_variable import Abs
 from extreme_data.meteo_france_data.scm_models_data.utils import ALTITUDES, ZS_INT_23, ZS_INT_MASK, LONGITUDES, \
     LATITUDES, ORIENTATIONS, SLOPES, ORDERED_ALLSLOPES_ALTITUDES, ORDERED_ALLSLOPES_ORIENTATIONS, \
     ORDERED_ALLSLOPES_SLOPES, ORDERED_ALLSLOPES_MASSIFNUM, date_to_str, WP_PATTERN_MAX_YEAR, Season, \
-    first_day_and_last_day, FrenchRegion, ZS_INT_MASK_PYRENNES, alps_massif_order, ZS_INT_MASK_PYRENNES_LIST
+    first_day_and_last_day, FrenchRegion, ZS_INT_MASK_PYRENNES, alps_massif_order, ZS_INT_MASK_PYRENNES_LIST, \
+    season_to_str
 from extreme_data.meteo_france_data.scm_models_data.visualization.utils import get_km_formatter
 from extreme_fit.function.margin_function.abstract_margin_function import \
     AbstractMarginFunction
@@ -665,6 +666,10 @@ class AbstractStudy(object):
     def has_orientation(self):
         return self.orientation is not None
 
+    @property
+    def season_name(self):
+        return season_to_str(self.season)
+
     """  Spatial properties """
 
     @cached_property
diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py
index 271753085c518ba9a2aeb21e86cb751b4c8787ce..fb55246b692b1a3dc91b1f671a2532bc0e2162b6 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/main_study_visualizer.py
@@ -2,8 +2,10 @@ import time
 from collections import OrderedDict
 from typing import List
 
+from extreme_data.meteo_france_data.adamont_data.snowfall_simulation import SafranSnowfallSimulationRCP85
 from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
-from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusDepth, CrocusSweTotal, ExtendedCrocusDepth, \
+from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusDepth, CrocusSweTotal, \
+    ExtendedCrocusDepth, \
     ExtendedCrocusSweTotal, CrocusDaysWithSnowOnGround, CrocusSwe3Days, CrocusSnowLoad3Days, CrocusSnowLoadTotal, \
     CrocusSnowLoadEurocode, CrocusSnowLoad5Days, CrocusSnowLoad7Days
 from extreme_data.meteo_france_data.scm_models_data.crocus.crocus_variables import CrocusDensityVariable
@@ -48,10 +50,12 @@ SCM_STUDY_CLASS_TO_ABBREVIATION = {
     CrocusSnowLoad5Days: 'GSL5',
     CrocusSnowLoad7Days: 'GSL7',
     CrocusSnowDensityAtMaxofSwe: '{} when the max of GSL \nis reached'.format(snow_density_str),
-    CrocusDifferenceSnowLoadRescaledAndEurocodeToSeeSynchronization: 'max GSL rescaled - GSL from max HS \nboth with {}'.format(eurocode_snow_density),
+    CrocusDifferenceSnowLoadRescaledAndEurocodeToSeeSynchronization: 'max GSL rescaled - GSL from max HS \nboth with {}'.format(
+        eurocode_snow_density),
     CrocusDifferenceSnowLoad: ('max GSL - GSL from max HS \n with {}'.format(eurocode_snow_density)),
     CrocusSnowDepthDifference: 'max HS - HS at max of GSL',
     CrocusSnowDepthAtMaxofSwe: 'HS at max of GSL',
+    SafranSnowfallSimulationRCP85: 'SF1 RCP85 projections'
 }
 
 altitude_massif_name_and_study_class_for_poster = [
@@ -177,7 +181,6 @@ def case_study():
         print(y)
 
 
-
 def complete_analysis(only_first_one=False):
     """An overview of everything that is possible with study OR extended study"""
     for study_class, extended_study_class in list(SCM_STUDY_TO_EXTENDED_STUDY.items())[:]:
@@ -244,7 +247,6 @@ def main_run():
     max_graph_annual_maxima_poster()
 
 
-
 if __name__ == '__main__':
     start = time.time()
     main_run()
diff --git a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
index 91c427a3e6011b63a99693240dc4c397e21d0502..037dcbff233395c0349c43afc5ba9d14352dc6cc 100644
--- a/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
+++ b/extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py
@@ -13,7 +13,6 @@ import numpy as np
 import pandas as pd
 import seaborn as sns
 
-from extreme_data.meteo_france_data.scm_models_data.utils import season_to_str
 from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import load_plot
 from extreme_fit.estimator.margin_estimator.utils import fitted_stationary_gev
 from extreme_fit.model.margin_model.utils import fitmethod_to_str
@@ -555,7 +554,7 @@ class StudyVisualizer(VisualizationParameters):
             plt.show()
         if self.save_to_file:
             main_title, specific_title = '_'.join(self.study.title.split()).split('/')
-            main_title += season_to_str(self.study.season)
+            main_title += self.study.season_name
             if folder_for_variable:
                 filename = "{}/{}/".format(VERSION_TIME, main_title)
             else:
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
index 75671432d34f946e60cd169e96d73d1e117f8e2b..ea4282f823a7781f9516dcbd599beb920b8d14f8 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/utils.py
@@ -34,7 +34,8 @@ from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_mode
     NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation, \
     NonStationaryAltitudinalScaleQuadraticLocationLinearCrossTermForLocation, \
     NonStationaryAltitudinalLocationCubicScaleQuadraticCrossTermForLocation, \
-    NonStationaryAltitudinalLocationCubicScaleLinearCrossTermForLocation
+    NonStationaryAltitudinalLocationCubicScaleLinearCrossTermForLocation, \
+    NonStationaryAltitudinalLocationCubicScaleLinear
 from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_models import \
     StationaryGumbelAltitudinal, NonStationaryGumbelAltitudinalScaleLinear, \
     NonStationaryGumbelAltitudinalLocationLinear, NonStationaryGumbelAltitudinalLocationQuadratic, \
@@ -89,8 +90,12 @@ ALTITUDINAL_GEV_MODELS_LOCATION_QUADRATIC_MINIMUM = [
     NonStationaryAltitudinalLocationCubicCrossTermForLocation,
     # NonStationaryAltitudinalLocationOrder4CrossTermForLocation,
 
-]
+    NonStationaryAltitudinalLocationQuadraticScaleLinearCrossTermForLocation,
+    NonStationaryAltitudinalLocationQuadraticScaleLinear,
+    NonStationaryAltitudinalLocationCubicScaleLinear,
+    NonStationaryAltitudinalLocationCubicScaleLinearCrossTermForLocation,
 
+]
 
 ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES = [
     StationaryAltitudinalOnlyScale,
diff --git a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
index 940d4fa3f1bfd064ca00af9c507a17f49373717f..4c4fffa47612cbfbaa761bcdd6e42028ccfa2c64 100644
--- a/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py
@@ -46,6 +46,8 @@ class AltitudesStudies(object):
                                 s_split_temporal: pd.Series = None, top_n_values_to_remove=None):
         coordinate_values_to_maxima = {}
         massif_altitudes = self.massif_name_to_altitudes[massif_name]
+        if len(massif_altitudes) == 0:
+            print('{} has no data for these altitudes: {}'.format(massif_name, self.altitudes))
         for altitude in massif_altitudes:
             study = self.altitude_to_study[altitude]
             for year, maxima in zip(study.ordered_years, study.massif_name_to_annual_maxima[massif_name]):
@@ -71,6 +73,7 @@ class AltitudesStudies(object):
         if massif_altitudes is None or set(massif_altitudes) == set(self.altitudes):
             spatial_coordinates = self.spatial_coordinates
         else:
+            assert len(massif_altitudes) > 0
             spatial_coordinates = self.spatial_coordinates_for_altitudes(massif_altitudes)
         slicer_class = get_slicer_class_from_s_splits(s_split_spatial, s_split_temporal)
         return AbstractSpatioTemporalCoordinates(slicer_class=slicer_class,
@@ -107,10 +110,10 @@ class AltitudesStudies(object):
 
     # Some visualization
 
-    def show_or_save_to_file(self, plot_name, show=False, no_title=False):
+    def show_or_save_to_file(self, plot_name, show=False, no_title=False, tight_layout=None):
         study_visualizer = StudyVisualizer(study=self.study, show=show, save_to_file=not show)
         study_visualizer.plot_name = plot_name
-        study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500, no_title=no_title)
+        study_visualizer.show_or_save_to_file(add_classic_title=False, dpi=500, no_title=no_title, tight_layout=tight_layout)
 
     def run_for_each_massif(self, function, massif_names, **kwargs):
         massif_names = massif_names if massif_names is not None else self.study.all_massif_names()
@@ -150,16 +153,17 @@ class AltitudesStudies(object):
         if change is True or change is None:
             moment += ' change (between two block of 30 years) for'
         moment += 'mean' if not std else 'std'
-        if change is False:
-            moment += ' (for the 60 years of data)'
-        plot_name = '{} of annual maxima of {}'.format(moment, SCM_STUDY_CLASS_TO_ABBREVIATION[self.study_class])
+        # if change is False:
+            # moment += ' (for the 60 years of data)'
+        plot_name = '{} of {} maxima of {}'.format(moment, self.study.season_name,
+                                                   SCM_STUDY_CLASS_TO_ABBREVIATION[self.study_class])
         ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
         ax.set_xlabel('altitudes', fontsize=15)
         lim_down, lim_up = ax.get_ylim()
         lim_up += (lim_up - lim_down) / 3
         ax.set_ylim([lim_down, lim_up])
         ax.tick_params(axis='both', which='major', labelsize=13)
-        self.show_or_save_to_file(plot_name=plot_name, show=show)
+        self.show_or_save_to_file(plot_name=plot_name, show=show, no_title=True)
         ax.clear()
 
     def _plot_mean_maxima_against_altitude(self, massif_name, massif_id, ax=None, std=False, change=False):
diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index aee48ba6151de3e16ac75dd1a28db49608dc6060..87dd6ed59bae8fb51dc672248431a43f19277c63 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -2,6 +2,8 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+from extreme_data.meteo_france_data.adamont_data.abstract_simulation_study import SimulationStudy
+from extreme_data.meteo_france_data.adamont_data.snowfall_simulation import SafranSnowfallSimulationRCP85
 from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUDINAL_GEV_MODELS, \
     ALTITUDINAL_GEV_MODELS_LOCATION_QUADRATIC_MINIMUM, ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES, ALTITUDINAL_GEV_MODELS_LOCATION
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_total_aic, \
@@ -70,21 +72,32 @@ def main():
                      SafranPrecipitation7Days][:]
     study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranPrecipitation1Day
                      , SafranPrecipitation3Days][:1]
-    # seasons = [Season.automn, Season.winter, Season.spring][::-1]
-    seasons = [Season.winter]
-    # seasons = [Season.winter_extended]
 
+    # altitudes = [1500, 1800]
+    study_classes = [SafranSnowfall1Day, SafranSnowfallSimulationRCP85][:]
+
+    # Common parameters
+    altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600]
     massif_names = None
-    # massif_names = ['Aravis']
-    # massif_names = ['Chartreuse', 'Belledonne']
+    seasons = [Season.winter]
 
     for season in seasons:
         for study_class in study_classes:
-            studies = AltitudesStudies(study_class, altitudes, season=season)
-            print('inner loop', season, study_class)
-            # plot_time_series(studies, massif_names)
-            # plot_moments(studies, massif_names)
-            plot_altitudinal_fit(studies, massif_names)
+            if issubclass(study_class, SimulationStudy):
+                for ensemble_idx in list(range(14))[:1]:
+                    studies = AltitudesStudies(study_class, altitudes, season=season,
+                                               ensemble_idx=ensemble_idx)
+                    plot_studies(massif_names, season, studies, study_class)
+            else:
+                studies = AltitudesStudies(study_class, altitudes, season=season)
+                plot_studies(massif_names, season, studies, study_class)
+
+
+def plot_studies(massif_names, season, studies, study_class):
+    print('inner loop', season, study_class)
+    # plot_time_series(studies, massif_names)
+    # plot_moments(studies, massif_names)
+    plot_altitudinal_fit(studies, massif_names)
 
 
 if __name__ == '__main__':
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index 1b57bf0c7bc1d8a7c776c723c8ee7a9c0b6f20ea..3c7e4d70aa9b6f37755f0168c5da28bd280e2683 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -142,16 +142,15 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
 
     def plot_shape_map(self):
         self.plot_abstract_fast(self.massif_name_to_shape,
-                                label='Shape plot for {} {}'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)],
-                                                                    self.study.variable_unit),
+                                label='Shape parameter for {} maxima of {}'.format(self.study.season_name,
+                                                                                   SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)]),
                                 add_x_label=False, graduation=0.1, massif_name_to_text=self.massif_name_to_name)
 
     def plot_altitude_for_the_peak(self):
         pass
 
     def plot_year_for_the_peak(self, plot_mean=True):
-        t_list = 1959 + np.arange(20 + 41)
-        # t_list = 1800 + np.arange(400)
+        t_list = self.study.ordered_years
         return_period = 50
         for massif_name, one_fold_fit in self.massif_name_to_one_fold_fit.items():
             ax = plt.gca()
@@ -182,14 +181,14 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
             ax.set_ylim(ax_lim)
             ax.set_xlabel('Year')
             if plot_mean:
-                ylabel = 'Mean annual maxima'
+                ylabel = 'Mean {} maxima'.format(self.study.season_name)
             else:
                 ylabel = '{}-year return level'.format(return_period)
-            ax.set_ylabel('{} of {}'.format(ylabel, SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)]))
+            ax.set_ylabel('{} of {} in {} ({})'.format(ylabel, SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)], massif_name.replace('_', ' '), self.study.variable_unit))
             peak_year_folder = 'Peak year ' + ylabel
             plot_name = '{}/{}/Peak year for {}'.format(OneFoldFit.folder_for_plots, peak_year_folder,
                                                         massif_name.replace('_', ''))
-            self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show)
+            self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show, no_title=True, tight_layout=True)
             plt.close()
 
     # Plots "altitude switch" and "peak year"
@@ -198,8 +197,9 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
     def massif_name_to_is_decreasing_parabol(self):
         # For the test we only activate the Mont-Blanc massif
         d = {massif_name: False for massif_name in self.massif_name_to_one_fold_fit.keys()}
-        for massif_name in ['Mont-Blanc', 'Vanoise', 'Aravis', 'Beaufortain', 'Chablais']:
-            d[massif_name] = True
+        if max(self.study.ordered_years) < 2030:
+            for massif_name in ['Vanoise', 'Aravis', 'Beaufortain', 'Chablais']:
+                d[massif_name] = True
         return d
 
     @property
diff --git a/projects/ogorman/gorman_figures/temperature_of_snowfall_maxima.py b/projects/ogorman/gorman_figures/temperature_of_snowfall_maxima.py
index 913f51aa302203f01dcb825d03cc66112c9ef2fc..bd518ce869a5a9843d8a08860b2693a5cc87fdd8 100644
--- a/projects/ogorman/gorman_figures/temperature_of_snowfall_maxima.py
+++ b/projects/ogorman/gorman_figures/temperature_of_snowfall_maxima.py
@@ -2,7 +2,7 @@ import matplotlib as mpl
 mpl.rcParams['text.usetex'] = True
 mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
 
-from extreme_data.meteo_france_data.scm_models_data.utils import Season, season_to_str
+from extreme_data.meteo_france_data.scm_models_data.utils import Season
 
 from collections import OrderedDict
 
@@ -51,7 +51,7 @@ def distribution_temperature_of_maxima_of_snowfall(altitudes, temperature_at_max
     ax.boxplot([altitude_to_maxima_temperature[a] for a in altitudes], positions=altitudes, widths=width)
     ax.set_xlim([min(altitudes) - width, max(altitudes) + width])
     suffix = 'at maxima of snowfall' if temperature_at_maxima else ''
-    ylabel = 'Daily {} temperature {} ({})'.format(season_to_str(season), suffix, temperature_study.variable_class.UNIT)
+    ylabel = 'Daily {} temperature {} ({})'.format(temperature_study.season_name, suffix, temperature_study.variable_class.UNIT)
     print(ylabel)
     ax.set_ylabel(ylabel)
     ax.set_xlabel('Altitude (m)')
diff --git a/test/test_extreme_data/test_meteo_france_data/test_SCM_study.py b/test/test_extreme_data/test_meteo_france_data/test_SCM_study.py
index 6d0af7e6221e58fe46c269e1211698e49c4f6645..f098eadea20552278960b2de69f65f9572c54342 100644
--- a/test/test_extreme_data/test_meteo_france_data/test_SCM_study.py
+++ b/test/test_extreme_data/test_meteo_france_data/test_SCM_study.py
@@ -5,6 +5,7 @@ from random import sample
 
 import pandas as pd
 
+from extreme_data.meteo_france_data.adamont_data.abstract_simulation_study import SimulationStudy
 from extreme_data.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoad3Days
 from extreme_data.meteo_france_data.scm_models_data.safran.cumulated_study import NB_DAYS
 from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall, SafranTemperature, \
@@ -63,7 +64,9 @@ class TestSCMAllStudy(unittest.TestCase):
     def test_variables(self):
         for study_class in SCM_STUDY_CLASS_TO_ABBREVIATION.keys():
             study = study_class(year_max=1959)
-            _ = study.year_to_annual_maxima[1959]
+            if not issubclass(study_class, SimulationStudy):
+                print(study_class)
+                _ = study.year_to_annual_maxima[1959]
         self.assertTrue(True)
 
 
diff --git a/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py b/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py
index 4f72286be42c8d1f92c015aefb786106fc7d8fc9..eb59e3fb4ae1179548960a4bb9ee9c0c1bbf5a48 100644
--- a/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py
+++ b/test/test_extreme_data/test_meteo_france_data/test_adamont_study.py
@@ -6,6 +6,7 @@ from extreme_data.meteo_france_data.adamont_data.snowfall_simulation import Safr
 class TestAdamontStudy(unittest.TestCase):
 
     def test_year_to_date(self):
+
         study = SafranSnowfallSimulationRCP85(altitude=900)
         self.assertTrue(True)