diff --git a/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py b/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py
index 8cd7e69d5f57d610b450c00f46809c1f8555d187..0748e83c02d846204cf3e67e785cb8d32fca8909 100644
--- a/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py
+++ b/extreme_data/meteo_france_data/adamont_data/cmip5/temperature_to_year.py
@@ -86,7 +86,7 @@ def plot_nb_data(is_temperature_interval, is_shift_interval):
     ]
     ax2.legend(handles=legend_elements, loc='upper center')
     ax2.set_yticks([])
-    # plt.show()
+    plt.show()
 
 
 def get_interval_limits(is_temperature_interval, is_shift_interval):
@@ -111,13 +111,13 @@ def get_ticks_labels_for_interval(is_temperature_interval, is_shift_interval):
     ticks_labels = [' +${}^o\mathrm{C}$ and +${}^o\mathrm{C}$'.format(left_limit, right_limit, **{'C': '{C}'})
                     if is_temperature_interval else '{} and {}'.format(left_limit, right_limit)
                     for left_limit, right_limit in zip(left_limits, right_limits)]
-    prefix = 'Maxima occured between \n'
+    prefix = 'Maxima between \n'
     ticks_labels = [prefix + l for l in ticks_labels]
     return ticks_labels
 
 
 if __name__ == '__main__':
     for shift_interval in [False, True]:
-        for temp_interval in [True, False]:
+        for temp_interval in [False, True]:
             print("shift = {}, temp_inteval = {}".format(shift_interval, temp_interval))
             plot_nb_data(is_temperature_interval=temp_interval, is_shift_interval=shift_interval)
diff --git a/extreme_trend/ensemble_fit/abstract_ensemble_fit.py b/extreme_trend/ensemble_fit/abstract_ensemble_fit.py
index 1a209af2e285f5a14c21a5730f1f704f7330acc6..5d352ea8e1d1a6b345b494fc77edfbf33ee07b59 100644
--- a/extreme_trend/ensemble_fit/abstract_ensemble_fit.py
+++ b/extreme_trend/ensemble_fit/abstract_ensemble_fit.py
@@ -2,9 +2,12 @@ from typing import Dict, Tuple
 
 from extreme_fit.model.margin_model.utils import MarginFitMethod
 from extreme_data.meteo_france_data.scm_models_data.altitudes_studies import AltitudesStudies
+from extreme_trend.one_fold_fit.one_fold_fit import OneFoldFit
 
 
 class AbstractEnsembleFit(object):
+    Median_merge = 'Median'
+    Mean_merge = 'Mean'
 
     def __init__(self, massif_names, gcm_rcm_couple_to_altitude_studies: Dict[Tuple[str, str], AltitudesStudies],
                  models_classes,
@@ -23,6 +26,14 @@ class AbstractEnsembleFit(object):
         self.only_models_that_pass_goodness_of_fit_test = only_models_that_pass_goodness_of_fit_test
         self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
 
+        # Set appropriate setting
+        OneFoldFit.last_year = 2100
+        OneFoldFit.nb_years = 95
+
     @property
     def altitudes(self):
+        raise self.visualizer_list.studies.altitudes
+
+    @property
+    def visualizer_list(self):
         raise NotImplementedError
\ No newline at end of file
diff --git a/extreme_trend/ensemble_fit/clustered_ensemble.py b/extreme_trend/ensemble_fit/clustered_ensemble.py
deleted file mode 100644
index 755905142c4fe40cdacc3326aca84b597cf51b68..0000000000000000000000000000000000000000
--- a/extreme_trend/ensemble_fit/clustered_ensemble.py
+++ /dev/null
@@ -1,19 +0,0 @@
-
-"""Instead of creating one big group, with all the ensemble together,
-and assuming a common temporal trend to all the group.
-
-We could create smaller groups, with an importance proportional to the number of GCM/RCM couples considered
-in the group.
-For instance, we could group them by GCM, or group them by RCM.
-Or we could try to find a metric to group them together.
-
-This is the idea of finding of sweet spot between:
--only independent fits with few assumptions
--one common fit with too much assumption
-
-
-it links with the idea of "climate model subset".
-
-Generally people try to find one model subset,
-the idea here, would be to find group of model subsets
-"""
\ No newline at end of file
diff --git a/extreme_trend/ensemble_fit/independent_ensemble_fit/independent_ensemble_fit.py b/extreme_trend/ensemble_fit/independent_ensemble_fit/independent_ensemble_fit.py
index c4c3aaca4989d7c45504ac6c3e63039d86bf7cae..42af498c8f343e217e5908b8f41b58aa63077896 100644
--- a/extreme_trend/ensemble_fit/independent_ensemble_fit/independent_ensemble_fit.py
+++ b/extreme_trend/ensemble_fit/independent_ensemble_fit/independent_ensemble_fit.py
@@ -9,15 +9,10 @@ from extreme_trend.one_fold_fit.one_fold_fit import OneFoldFit
 
 class IndependentEnsembleFit(AbstractEnsembleFit):
     """For each gcm_rcm_couple, we create a OneFoldFit"""
-    Median_merge = 'Median'
-    Mean_merge = 'Mean'
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
-        # Set appropriate setting
-        OneFoldFit.last_year = 2100
-        OneFoldFit.nb_years = 95
         # Load a classical visualizer
         self.gcm_rcm_couple_to_visualizer = {}
         for gcm_rcm_couple, studies in self.gcm_rcm_couple_to_altitude_studies.items():
@@ -47,9 +42,7 @@ class IndependentEnsembleFit(AbstractEnsembleFit):
         }
 
     @property
-    def visualizer(self):
-        return list(self.gcm_rcm_couple_to_visualizer.values())[0]
+    def visualizer_list(self):
+        return list(self.gcm_rcm_couple_to_visualizer.values()) \
+               + list(self.merge_function_name_to_visualizer.values())
 
-    @property
-    def altitudes(self):
-        return self.visualizer.studies.altitudes
diff --git a/extreme_trend/ensemble_fit/together_ensemble_fit.py b/extreme_trend/ensemble_fit/together_ensemble_fit/__init__.py
similarity index 100%
rename from extreme_trend/ensemble_fit/together_ensemble_fit.py
rename to extreme_trend/ensemble_fit/together_ensemble_fit/__init__.py
diff --git a/extreme_trend/ensemble_fit/together_ensemble_fit/together_ensemble_fit.py b/extreme_trend/ensemble_fit/together_ensemble_fit/together_ensemble_fit.py
new file mode 100644
index 0000000000000000000000000000000000000000..7534c8a2da7c705d5801696a016ab0694c1d062a
--- /dev/null
+++ b/extreme_trend/ensemble_fit/together_ensemble_fit/together_ensemble_fit.py
@@ -0,0 +1,24 @@
+from extreme_trend.ensemble_fit.abstract_ensemble_fit import AbstractEnsembleFit
+from extreme_trend.ensemble_fit.together_ensemble_fit.visualizer_non_stationary_ensemble import \
+    VisualizerNonStationaryEnsemble
+
+
+class TogetherEnsembleFit(AbstractEnsembleFit):
+    """We create a single OneFoldFit for all gcm_rcm_couples"""
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.visualizer = VisualizerNonStationaryEnsemble(self.gcm_rcm_couple_to_altitude_studies,
+                                                     self.models_classes,
+                                                     False,
+                                                     self.massif_names, self.fit_method,
+                                                     self.temporal_covariate_for_fit,
+                                                     self.only_models_that_pass_goodness_of_fit_test,
+                                                     self.confidence_interval_based_on_delta_method,
+                                                     self.remove_physically_implausible_models
+                                                     )
+
+    @property
+    def visualizer_list(self):
+        return [self.visualizer]
+
diff --git a/extreme_trend/ensemble_fit/together_ensemble_fit/visualizer_non_stationary_ensemble.py b/extreme_trend/ensemble_fit/together_ensemble_fit/visualizer_non_stationary_ensemble.py
new file mode 100644
index 0000000000000000000000000000000000000000..8781b583d6523d76c21799556d537f4f08a66779
--- /dev/null
+++ b/extreme_trend/ensemble_fit/together_ensemble_fit/visualizer_non_stationary_ensemble.py
@@ -0,0 +1,28 @@
+from typing import List, Dict, Tuple
+
+from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
+from extreme_data.meteo_france_data.scm_models_data.altitudes_studies import AltitudesStudies
+from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
+    AbstractSpatioTemporalPolynomialModel
+from extreme_fit.model.margin_model.utils import MarginFitMethod
+from extreme_trend.one_fold_fit.altitudes_studies_visualizer_for_non_stationary_models import \
+    AltitudesStudiesVisualizerForNonStationaryModels
+from extreme_trend.trend_test.visualizers.study_visualizer_for_non_stationary_trends import \
+    StudyVisualizerForNonStationaryTrends
+
+
+class VisualizerNonStationaryEnsemble(AltitudesStudiesVisualizerForNonStationaryModels):
+
+    def __init__(self, gcm_rcm_couple_to_studies: Dict[Tuple[str, str], AltitudesStudies], *args, **kwargs):
+        self.gcm_rcm_couple_to_studies = gcm_rcm_couple_to_studies
+        studies = list(self.gcm_rcm_couple_to_studies.values())[0]
+        super().__init__(studies, *args, **kwargs)
+
+    def get_massif_altitudes(self, massif_name):
+        # return self._get_massif_altitudes(massif_name, self.studies)
+        raise NotImplementedError
+
+    def get_dataset(self, massif_altitudes, massif_name):
+        raise NotImplementedError
+        # dataset = self.studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
+        # return dataset
diff --git a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
index 561387abc63859aae890198e9124f600c13c5031..665edb893d4cf38ea280b3d3b974c8e5441dd243 100644
--- a/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
+++ b/extreme_trend/ensemble_fit/visualizer_for_projection_ensemble.py
@@ -5,7 +5,9 @@ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_poly
     AbstractSpatioTemporalPolynomialModel
 from extreme_fit.model.margin_model.utils import MarginFitMethod
 from extreme_data.meteo_france_data.scm_models_data.altitudes_studies import AltitudesStudies
+from extreme_trend.ensemble_fit.abstract_ensemble_fit import AbstractEnsembleFit
 from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fit import IndependentEnsembleFit
+from extreme_trend.ensemble_fit.together_ensemble_fit.together_ensemble_fit import TogetherEnsembleFit
 from extreme_trend.one_fold_fit.altitude_group import get_altitude_class_from_altitudes
 from extreme_trend.one_fold_fit.plots.plot_histogram_altitude_studies import \
     plot_histogram_all_trends_against_altitudes, plot_shoe_plot_changes_against_altitude
@@ -75,29 +77,35 @@ class VisualizerForProjectionEnsemble(object):
                 ensemble_class_to_ensemble_fit[ensemble_fit_class] = ensemble_fit
             self.altitude_class_to_ensemble_class_to_ensemble_fit[altitude_class] = ensemble_class_to_ensemble_fit
 
+    def plot_for_visualizer_list(self, visualizer_list):
+        with_significance = False
+        for v in visualizer_list:
+            v.plot_moments()
+        plot_histogram_all_trends_against_altitudes(self.massif_names, visualizer_list,
+                                                    with_significance=with_significance)
+        for relative in [True, False]:
+            plot_shoe_plot_changes_against_altitude(self.massif_names, visualizer_list, relative=relative,
+                                                    with_significance=with_significance)
+
     def plot(self):
+        # Set limit for the plot
+        visualizer_list = []
+        for ensemble_fit_class in self.ensemble_fit_classes:
+            for ensemble_fit in self.ensemble_fits(ensemble_fit_class):
+                visualizer_list.extend(ensemble_fit.visualizer_list)
+        compute_and_assign_max_abs(visualizer_list)
+        # Plot
         if IndependentEnsembleFit in self.ensemble_fit_classes:
-            # Set max abs
-            visualizer_list = []
-            for ensemble_fit in self.ensemble_fits(IndependentEnsembleFit):
-                visualizer_list.extend(list(ensemble_fit.gcm_rcm_couple_to_visualizer.values()))
-            # Potentially I could add more visualizer here...
-            method_name_and_order_to_max_abs, max_abs_for_shape = compute_and_assign_max_abs(visualizer_list)
-            # Assign the same max abs for the 
-            for ensemble_fit in self.ensemble_fits(IndependentEnsembleFit):
-                for v in ensemble_fit.merge_function_name_to_visualizer.values():
-                    v._max_abs_for_shape = max_abs_for_shape
-                    v._method_name_and_order_to_max_abs = method_name_and_order_to_max_abs
-            # Plot
             self.plot_independent()
+        if TogetherEnsembleFit in self.ensemble_fit_classes:
+            self.plot_together()
 
     def plot_independent(self):
-        with_significance = False
         # Aggregated at gcm_rcm_level plots
-        merge_keys = [IndependentEnsembleFit.Median_merge, IndependentEnsembleFit.Mean_merge]
+        merge_keys = [AbstractEnsembleFit.Median_merge, AbstractEnsembleFit.Mean_merge]
         keys = self.gcm_rcm_couples + merge_keys
         # Only plot Mean for speed
-        keys = [IndependentEnsembleFit.Mean_merge]
+        keys = [AbstractEnsembleFit.Mean_merge]
         for key in keys:
             visualizer_list = [independent_ensemble_fit.gcm_rcm_couple_to_visualizer[key]
                                if key in self.gcm_rcm_couples
@@ -107,17 +115,15 @@ class VisualizerForProjectionEnsemble(object):
             if key in merge_keys:
                 for v in visualizer_list:
                     v.studies.study.gcm_rcm_couple = (key, "merge")
-            for v in visualizer_list:
-                v.plot_moments()
-            plot_histogram_all_trends_against_altitudes(self.massif_names, visualizer_list, with_significance=with_significance)
-            for relative in [True, False]:
-                plot_shoe_plot_changes_against_altitude(self.massif_names, visualizer_list, relative=relative, with_significance=with_significance)
+            self.plot_for_visualizer_list(visualizer_list)
+
+    def plot_together(self):
+        visualizer_list = [together_ensemble_fit.visualizer
+                           for together_ensemble_fit in self.ensemble_fits(TogetherEnsembleFit)]
+        self.plot_for_visualizer_list(visualizer_list)
 
     def ensemble_fits(self, ensemble_class):
         """Return the ordered ensemble fit for a given ensemble class (in the order of the altitudes)"""
         return [ensemble_class_to_ensemble_fit[ensemble_class]
                 for ensemble_class_to_ensemble_fit
                 in self.altitude_class_to_ensemble_class_to_ensemble_fit.values()]
-
-    def plot_together(self):
-        pass
diff --git a/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py b/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py
index 222fbc79e9e275c130d5315c0f3fe32385c41e15..8e8c7f10b83ca97eb919a0db307fecf74a00dc42 100644
--- a/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py
+++ b/extreme_trend/ensemble_fit/visualizer_for_sensitivity.py
@@ -8,6 +8,7 @@ from extreme_data.meteo_france_data.scm_models_data.utils import Season
 from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
     AbstractSpatioTemporalPolynomialModel
 from extreme_fit.model.margin_model.utils import MarginFitMethod
+from extreme_trend.ensemble_fit.abstract_ensemble_fit import AbstractEnsembleFit
 from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fit import IndependentEnsembleFit
 from extreme_trend.ensemble_fit.visualizer_for_projection_ensemble import VisualizerForProjectionEnsemble
 from extreme_trend.one_fold_fit.altitude_group import get_altitude_class_from_altitudes, \
@@ -25,7 +26,7 @@ class VisualizerForSensivity(object):
                  display_only_model_that_pass_gof_test=False,
                  confidence_interval_based_on_delta_method=False,
                  remove_physically_implausible_models=False,
-                 merge_visualizer_str=IndependentEnsembleFit.Median_merge,  # if we choose the Mean merge, then it is almost impossible to obtain stationary trends
+                 merge_visualizer_str=AbstractEnsembleFit.Median_merge,  # if we choose the Mean merge, then it is almost impossible to obtain stationary trends
                  is_temperature_interval=False,
                  is_shift_interval=False,
                  ):
diff --git a/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py b/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
index 6263daed5f8e63eb8c4ffed5acc31ba78a989d2f..c555a33cb111435b5ad80f0b1b965a0322b6e5cb 100644
--- a/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/extreme_trend/one_fold_fit/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -79,7 +79,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
             # Save the massif altitudes only for those who pass the condition
             self.massif_name_to_massif_altitudes[massif_name] = massif_altitudes
             # Load dataset
-            dataset = self.studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
+            dataset = self.get_dataset(massif_altitudes, massif_name)
             old_fold_fit = OneFoldFit(massif_name, dataset, self.model_classes, self.fit_method,
                                       self.temporal_covariate_for_fit,
                                       type(self.altitude_group),
@@ -90,15 +90,22 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         else:
             return None
 
+    def get_dataset(self, massif_altitudes, massif_name):
+        dataset = self.studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
+        return dataset
+
     moment_names = ['moment', 'changes_of_moment', 'relative_changes_of_moment'][:]
     orders = [1, 2, None][2:]
 
     def get_massif_altitudes(self, massif_name):
-        valid_altitudes = [altitude for altitude, study in self.studies.altitude_to_study.items()
+        return self._get_massif_altitudes(massif_name, self.studies)
+
+    def _get_massif_altitudes(self, massif_name, studies):
+        valid_altitudes = [altitude for altitude, study in studies.altitude_to_study.items()
                            if massif_name in study.study_massif_names]
         massif_altitudes = []
         for altitude in valid_altitudes:
-            study = self.studies.altitude_to_study[altitude]
+            study = studies.altitude_to_study[altitude]
             annual_maxima = study.massif_name_to_annual_maxima[massif_name]
             percentage_of_non_zeros = 100 * np.count_nonzero(annual_maxima) / len(annual_maxima)
             if percentage_of_non_zeros > 90:
diff --git a/projects/projected_extreme_snowfall/evaluation/main_comparison_on_quantile_period.py b/projects/projected_extreme_snowfall/evaluation/main_comparison_on_quantile_period.py
index ff21e6353c9ca157b108a7bb1e6851d1d863ce18..8738cbc315ab6f308c74fafbc7583a114ba9ff4b 100644
--- a/projects/projected_extreme_snowfall/evaluation/main_comparison_on_quantile_period.py
+++ b/projects/projected_extreme_snowfall/evaluation/main_comparison_on_quantile_period.py
@@ -135,7 +135,7 @@ def main_comparaison_plot():
                     ax.set_ylabel('Altitude (m)', fontsize=10)
                     massif_str = 'all massifs' if massif_names is None else 'the {} massif'.format(massif_names[0])
                     unit = '%' if relative_bias else study.variable_unit
-                    bias_name = 'Relative bias' if relative_bias else 'Bias'
+                    bias_name = 'Relative difference' if relative_bias else 'Difference'
                     mean_str = 'mean' if mean else 'std'
                     title = '{} in the {} annual maxima of {} of {}\n' \
                                      'for ADAMONT v{}' \
@@ -144,7 +144,7 @@ def main_comparaison_plot():
                                                                                               adamont_version,
                                                                                               comparaison_study_class,
                                                                                               unit)
-                    folder = 'relative bias' if relative_bias else 'absolute bias'
+                    folder = 'relative difference' if relative_bias else 'difference'
                     plot_name = op.join(folder, title)
                     ax.set_xlabel(title, fontsize=10)
                     reanalysis_altitude_studies.show_or_save_to_file(plot_name=plot_name, no_title=True)
diff --git a/projects/projected_extreme_snowfall/main_elevation_temporal_for_projections_ensemble.py b/projects/projected_extreme_snowfall/main_projections_ensemble.py
similarity index 69%
rename from projects/projected_extreme_snowfall/main_elevation_temporal_for_projections_ensemble.py
rename to projects/projected_extreme_snowfall/main_projections_ensemble.py
index 5e1eac74b84f6e07368d55afd041b62f0a4cf27c..9cb4e279116f88220c3e06464ada16401c023565 100644
--- a/projects/projected_extreme_snowfall/main_elevation_temporal_for_projections_ensemble.py
+++ b/projects/projected_extreme_snowfall/main_projections_ensemble.py
@@ -1,9 +1,15 @@
 import datetime
 import time
 from typing import List
+import matplotlib
+
+from extreme_trend.ensemble_fit.together_ensemble_fit.together_ensemble_fit import TogetherEnsembleFit
 
+matplotlib.use('Agg')
 import matplotlib as mpl
 
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
 
 from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fit import IndependentEnsembleFit
 from extreme_trend.ensemble_fit.visualizer_for_projection_ensemble import VisualizerForProjectionEnsemble
@@ -11,20 +17,16 @@ from extreme_trend.ensemble_fit.visualizer_for_sensitivity import VisualizerForS
 from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
     AnomalyTemperatureWithSplineTemporalCovariate
 
-mpl.rcParams['text.usetex'] = True
-mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
-
 from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
     ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
-import matplotlib
+
 from extreme_fit.model.utils import set_seed_for_test
 
 from extreme_data.meteo_france_data.adamont_data.adamont.adamont_safran import AdamontSnowfall
 from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \
     rcp_scenarios
-from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import TimeTemporalCovariate
-
-matplotlib.use('Agg')
+from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
+    TimeTemporalCovariate
 
 from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
     AbstractExtractEurocodeReturnLevel
@@ -37,21 +39,20 @@ from extreme_data.meteo_france_data.scm_models_data.utils import Season
 def main():
     start = time.time()
     study_class = AdamontSnowfall
-    ensemble_fit_class = [IndependentEnsembleFit]
+    ensemble_fit_classes = [IndependentEnsembleFit, TogetherEnsembleFit][1:]
     temporal_covariate_for_fit = [TimeTemporalCovariate,
                                   AnomalyTemperatureWithSplineTemporalCovariate][0]
     set_seed_for_test()
     AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
 
     fast = True
-    sensitivity_plot = True
     scenarios = rcp_scenarios[::-1] if fast is False else [AdamontScenario.rcp85]
 
     for scenario in scenarios:
         gcm_rcm_couples = get_gcm_rcm_couples(scenario)
         if fast is None:
             massif_names = None
-            gcm_rcm_couples = None
+            gcm_rcm_couples = gcm_rcm_couples
             AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
             altitudes_list = altitudes_for_groups[3:]
         elif fast:
@@ -65,35 +66,15 @@ def main():
 
         assert isinstance(gcm_rcm_couples, list)
 
-        main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_class, ensemble_fit_class, scenario,
-                  temporal_covariate_for_fit, sensitivity_plot=sensitivity_plot)
+        assert isinstance(altitudes_list, List)
+        assert isinstance(altitudes_list[0], List)
+        print('Scenario is', scenario)
+        print('Covariate is {}'.format(temporal_covariate_for_fit))
 
-    end = time.time()
-    duration = str(datetime.timedelta(seconds=end - start))
-    print('Total duration', duration)
+        model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
+        assert scenario in rcp_scenarios
+        remove_physically_implausible_models = True
 
-
-def main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_class, ensemble_fit_classes, scenario,
-              temporal_covariate_for_fit, sensitivity_plot=False):
-    assert isinstance(altitudes_list, List)
-    assert isinstance(altitudes_list[0], List)
-    print('Scenario is', scenario)
-    print('Covariate is {}'.format(temporal_covariate_for_fit))
-
-    model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
-    assert scenario in rcp_scenarios
-    remove_physically_implausible_models = True
-
-    if sensitivity_plot:
-        visualizer = VisualizerForSensivity(
-            altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,
-            model_classes=model_classes,
-            ensemble_fit_classes=ensemble_fit_classes,
-            massif_names=massif_names,
-            temporal_covariate_for_fit=temporal_covariate_for_fit,
-            remove_physically_implausible_models=remove_physically_implausible_models,
-        )
-    else:
         visualizer = VisualizerForProjectionEnsemble(
             altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,
             model_classes=model_classes,
@@ -103,7 +84,11 @@ def main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_class, ensemb
             remove_physically_implausible_models=remove_physically_implausible_models,
             gcm_to_year_min_and_year_max=None,
         )
-    visualizer.plot()
+        visualizer.plot()
+
+    end = time.time()
+    duration = str(datetime.timedelta(seconds=end - start))
+    print('Total duration', duration)
 
 
 if __name__ == '__main__':
diff --git a/projects/projected_extreme_snowfall/main_sensitivity.py b/projects/projected_extreme_snowfall/main_sensitivity.py
new file mode 100644
index 0000000000000000000000000000000000000000..4ec37872e7f3f5b305c90d5511e7b8af3658c77b
--- /dev/null
+++ b/projects/projected_extreme_snowfall/main_sensitivity.py
@@ -0,0 +1,99 @@
+import datetime
+import time
+from typing import List
+import matplotlib
+
+from extreme_trend.ensemble_fit.abstract_ensemble_fit import AbstractEnsembleFit
+
+matplotlib.use('Agg')
+import matplotlib as mpl
+mpl.rcParams['text.usetex'] = True
+mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
+
+
+from extreme_trend.ensemble_fit.independent_ensemble_fit.independent_ensemble_fit import IndependentEnsembleFit
+from extreme_trend.ensemble_fit.visualizer_for_projection_ensemble import VisualizerForProjectionEnsemble
+from extreme_trend.ensemble_fit.visualizer_for_sensitivity import VisualizerForSensivity
+from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
+    AnomalyTemperatureWithSplineTemporalCovariate
+
+
+from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
+    ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
+
+from extreme_fit.model.utils import set_seed_for_test
+
+from extreme_data.meteo_france_data.adamont_data.adamont.adamont_safran import AdamontSnowfall
+from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \
+    rcp_scenarios
+from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import TimeTemporalCovariate
+
+
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
+    AbstractExtractEurocodeReturnLevel
+
+from extreme_trend.one_fold_fit.altitude_group import altitudes_for_groups
+
+from extreme_data.meteo_france_data.scm_models_data.utils import Season
+
+
+def main():
+    start = time.time()
+    study_class = AdamontSnowfall
+    ensemble_fit_classes = [IndependentEnsembleFit]
+    temporal_covariate_for_fit = [TimeTemporalCovariate,
+                                  AnomalyTemperatureWithSplineTemporalCovariate][0]
+    set_seed_for_test()
+    AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
+
+    fast = False
+    scenarios = rcp_scenarios[-1:] if fast is False else [AdamontScenario.rcp85]
+
+    for scenario in scenarios:
+        gcm_rcm_couples = get_gcm_rcm_couples(scenario)
+        if fast is None:
+            massif_names = None
+            gcm_rcm_couples = gcm_rcm_couples
+            AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
+            altitudes_list = altitudes_for_groups[1:2]
+        elif fast:
+            massif_names = ['Vanoise', 'Haute-Maurienne']
+            gcm_rcm_couples = gcm_rcm_couples[4:6]
+            AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
+            altitudes_list = altitudes_for_groups[:1]
+        else:
+            massif_names = None
+            altitudes_list = altitudes_for_groups[:]
+
+        assert isinstance(gcm_rcm_couples, list)
+
+        assert isinstance(altitudes_list, List)
+        assert isinstance(altitudes_list[0], List)
+        print('Scenario is', scenario)
+        print('Covariate is {}'.format(temporal_covariate_for_fit))
+
+        model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
+        assert scenario in rcp_scenarios
+        remove_physically_implausible_models = True
+
+        visualizer = VisualizerForSensivity(
+            altitudes_list, gcm_rcm_couples, study_class, Season.annual, scenario,
+            model_classes=model_classes,
+            ensemble_fit_classes=ensemble_fit_classes,
+            massif_names=massif_names,
+            temporal_covariate_for_fit=temporal_covariate_for_fit,
+            remove_physically_implausible_models=remove_physically_implausible_models,
+            merge_visualizer_str=AbstractEnsembleFit.Median_merge,
+            is_temperature_interval=False,
+            is_shift_interval=True,
+        )
+        visualizer.plot()
+
+    end = time.time()
+    duration = str(datetime.timedelta(seconds=end - start))
+    print('Total duration', duration)
+
+
+
+if __name__ == '__main__':
+    main()