diff --git a/climix/etc/metadata.yml b/climix/etc/metadata.yml
index d2f20717010ddd5b37049f3ab092f6c6af7de269..214c9928509d8a3916cce38332e2918bfac8e7f0 100644
--- a/climix/etc/metadata.yml
+++ b/climix/etc/metadata.yml
@@ -20,6 +20,20 @@ index_functions:
       condition:
         kind: operator
 
+last_occurrence:
+    description: |
+      Calculates the last time some condition is met.
+      First, the threshold is transformed to the same standard_name and units as
+      the input data.
+      Then the thresholding is performed as condition(data, threshold), ie
+      if condition is <, data < threshold.
+      Finally, locate the last occurrence when condition is met.
+    parameters:
+      threshold:
+        kind: quantity
+      condition:
+        kind: operator
+
   spell_length:
     description: |
       Calculates statistics on lengths of spells.
diff --git a/climix/index_functions.py b/climix/index_functions.py
index 033ef41cb25bef3fc164c333b72306fdbd61160a..21bb71d0990e9894aaec4c41baa20f7413d98f06 100644
--- a/climix/index_functions.py
+++ b/climix/index_functions.py
@@ -82,6 +82,33 @@ class CountOccurrences:
         return res.astype('float32')
 
 
+class LastOccurence:
+    def __init__(self, threshold, condition, reducer):
+        self.threshold = threshold
+        self.condition = condition
+        self.reducer = reducer
+        self.invalid_val = -1  # preliminary
+
+    def prepare(self, input_cube):
+        change_units(self.threshold, input_cube.units, input_cube.standard_name)
+
+    def call_func(self, data, axis, **kwargs):
+        ndays = data.shape[axis]
+        # need to "flip" data and subtract from year_length (i.e. array length)
+        data1 = self.condition(np.flip(data, axis=axis), self.threshold.data)
+        res = ndays - np.where(data1.any(), data1.argmax(), self.invalid_val)
+        return float(res)
+
+    def lazy_func(self, data, axis, **kwargs):
+        ndays = data.shape[axis]
+        # need to "flip" data and subtract from year_length (i.e. array length)
+        data1 = self.condition(da.flip(data, axis=axis), self.threshold.data)
+        res = ndays - da.where(data1.any(), data1.argmax(), self.invalid_val)
+        return float(res)
+
+
+
+
 class SpellLength:
     def __init__(self, threshold, condition, reducer):
         self.threshold = threshold
diff --git a/setup.py b/setup.py
index 83194163ec95d6105bd1c1395ef09f81a0881de4..69b39f251d7557a81fa66dc2567eb5f3e4c08312 100755
--- a/setup.py
+++ b/setup.py
@@ -50,6 +50,7 @@ setuptools.setup(
         ],
         'climix.index_functions': [
             'count_occurrences=climix.index_functions:CountOccurrences',
+            'last_occurrence=climix.index_functions:LastOccurence',
             'spell_length=climix.index_functions:SpellLength',
             'thresholded_statistics=climix.index_functions:ThresholdedStatistics',
         ],