diff --git a/probabilistic-analysis/probabilistic_assessment.py b/probabilistic-analysis/probabilistic_assessment.py
index c52fb91..98a97d1 100644
--- a/probabilistic-analysis/probabilistic_assessment.py
+++ b/probabilistic-analysis/probabilistic_assessment.py
@@ -169,6 +169,7 @@ def get_ongoing_recession(n_runs, start_year, end_year, sea_level_rise,
             'min'  (array_like): minimum value
             'mode' (array_like): most likely value
             'max'  (array_like): maximum value
+            'function'  (str): optional external function ('package.function')
         bruun_factor (dict):
             'min'  (float): minimum value
             'mode' (float): most likely value
@@ -180,6 +181,12 @@ def get_ongoing_recession(n_runs, start_year, end_year, sea_level_rise,
 
     Returns:
         the simulated recession distance (m)
+
+    Notes:
+    'sea_level_rise' is calculated with a triangular probability distribution
+    by default. Alternatively 'sea_level_rise' can be calculated using an
+    external function, to which the arguments 'n_runs', 'start_year', and
+    'end_year' are passed.
     """
 
     # Get time interval from input file
@@ -187,45 +194,54 @@ def get_ongoing_recession(n_runs, start_year, end_year, sea_level_rise,
     n_years = len(years)
 
     # Interpolate sea level rise projections (m)
-    slr_mode = np.interp(years,
-                         xp=sea_level_rise['year'],
-                         fp=sea_level_rise['mode'])[:, np.newaxis]
-
-    try:
-        slr_min = np.interp(years,
-                            xp=sea_level_rise['year'],
-                            fp=sea_level_rise['min'])[:, np.newaxis]
-    except ValueError:
-        # Use mode for deterministic beaches
-        slr_min = slr_mode
-
-    try:
-        slr_max = np.interp(years,
-                            xp=sea_level_rise['year'],
-                            fp=sea_level_rise['max'])[:, np.newaxis]
-    except ValueError:
-        # Use mode for deterministic beaches
-        slr_max = slr_mode
-
-    # Initialise sea level rise array
-    slr = np.zeros([n_years, n_runs])
-
-    for i in range(n_years):
-        # Use triangular distribution for SLR in each year (m)
+    if sea_level_rise['function']:  # Get slr from separate function
+        # Get names of package/script and function
+        pkg, func_name = sea_level_rise['function'].split('.')
+
+        # Import function from package
+        func = getattr(__import__(pkg), func_name)
+        slr = func(n_runs=n_runs, start_year=start_year, end_year=end_year)
+
+    else:  # Use triangular distribution
+        slr_mode = np.interp(years,
+                             xp=sea_level_rise['year'],
+                             fp=sea_level_rise['mode'])[:, np.newaxis]
+
         try:
-            slr[i, :] = np.random.triangular(left=slr_min[i],
-                                             mode=slr_mode[i],
-                                             right=slr_max[i],
-                                             size=n_runs)
+            slr_min = np.interp(years,
+                                xp=sea_level_rise['year'],
+                                fp=sea_level_rise['min'])[:, np.newaxis]
         except ValueError:
-            # Use constant value if slr_min == slr_max
-            slr[i, :] = np.ones([1, n_runs]) * slr_mode[i]
+            # Use mode for deterministic beaches
+            slr_min = slr_mode
 
-    # Sort each row, so SLR follows a smooth trajectory for each model run
-    slr = np.sort(slr, axis=1)
-
-    # Shuffle columns, so the order of model runs is randomised
-    slr = np.random.permutation(slr.T).T
+        try:
+            slr_max = np.interp(years,
+                                xp=sea_level_rise['year'],
+                                fp=sea_level_rise['max'])[:, np.newaxis]
+        except ValueError:
+            # Use mode for deterministic beaches
+            slr_max = slr_mode
+
+        # Initialise sea level rise array
+        slr = np.zeros([n_years, n_runs])
+
+        for i in range(n_years):
+            # Use triangular distribution for SLR in each year (m)
+            try:
+                slr[i, :] = np.random.triangular(left=slr_min[i],
+                                                 mode=slr_mode[i],
+                                                 right=slr_max[i],
+                                                 size=n_runs)
+            except ValueError:
+                # Use constant value if slr_min == slr_max
+                slr[i, :] = np.ones([1, n_runs]) * slr_mode[i]
+
+        # Sort each row, so SLR follows a smooth trajectory for each model run
+        slr = np.sort(slr, axis=1)
+
+        # Shuffle columns, so the order of model runs is randomised
+        slr = np.random.permutation(slr.T).T
 
     # Shift sea level so it is zero in the start year
     slr -= slr[0, :].mean(axis=0)