diff --git a/src/midstream_model_testing.m b/src/midstream_model_testing.m
index 5fa15fd..6c9b041 100644
--- a/src/midstream_model_testing.m
+++ b/src/midstream_model_testing.m
@@ -1,7 +1,7 @@
 function[WarningMessages, Exit_ID] = midstream_model_testing(Training_Model_Samples, ...
     Testing_dataset_for_testing, model, Model_Selection, ir_units, ...
     iet_units, ir_threshold, iet_threshold, Dt_fcst, FcastTime_window, ...
-    Path_output_results, HighPrctile, LowPrctile, ScaleNumEventForecast, ...
+    path2SaveOutData, HighPrctile, LowPrctile, ScaleNumEventForecast, ...
     Script_Mode)
 
 
@@ -115,15 +115,6 @@ end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Plot results, save output data and plots
 
-%%
-%Create folder 
-path2SaveOutData = [Path_output_results model.folder{Model_Selection} '/Forecast_Nev_per_' cell2mat(FcastTime_window)];
-try
-    mkdir(path2SaveOutData)
-catch
-    disp('Forder for output data already exists')
-end
-
 %%
 figure(1)
 set(1, 'Position', [943   264   964   670])
diff --git a/src/midstream_model_testing_wrapper.m b/src/midstream_model_testing_wrapper.m
new file mode 100644
index 0000000..2f2b759
--- /dev/null
+++ b/src/midstream_model_testing_wrapper.m
@@ -0,0 +1,44 @@
+function midstream_model_testing_wrapper(Training_Model_Samples, testing_data_path, out_dir, Model_Selection, ir_units, ir_threshold, iet_units, ...
+  iet_threshold, time_window, HighPrctile, LowPrctile, ScaleNumEventForecast, Run_Mode)
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  %% Implemented models - general settings
+  model.name = {'Stationary', 'Power Law', 'Power Function'};
+  model.ID = {1, 2, 3};
+  model.Ord_mu = {1, 2, 2}; % Number of model parameters (Ord_mu)
+  model.function = {'mu = a', 'mu = a*(IR)^b', 'mu = a*10^(b*IR)'}; %to be written in output files, as reference
+  model.folder = {'Stationary', 'PowerLaw', 'PowerFunction'}; % names of folders to be crated (related to model names)
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  Script_Mode.name = {'Testing_Forecasting', 'Custom_Forecasting'};
+  Script_Mode.ID = Run_Mode;
+  Script_Mode.Description = {'Use this mode to test the model forecasting and compare with observations using the testing dataset', ...
+      'Use this mode to forecast seismcity rates using a custom time series of injection rate data'};
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  ForecastSetting.Dt_fcst_options = [1/(24*60) 1/24 1 7];
+  ForecastSetting.Name = {'minute', 'hour', 'day'};
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  %% Create folders for writing output data
+  mkdir(out_dir);
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  [out] = midstream_model_testing(Training_Model_Samples, ...
+      testing_data_path, model, Model_Selection, ...
+      ir_units, iet_units, ir_threshold, iet_threshold,  ...
+      ForecastSetting.Dt_fcst_options(time_window), ...
+      ForecastSetting.Name(time_window), out_dir, ...
+      HighPrctile, LowPrctile, ScaleNumEventForecast, Script_Mode);
+
+  disp('#################################################################')
+  disp('Output summary:')
+  disp('Eventual warning messages')
+  for i=1:length(out)
+      disp(cell2mat( out(i) ))
+  end
+  disp('#################################################################')
+
+end