Repository: Ouranosinc/xclim

Issue:
PCC indicators
### Description

Indicators used by PCC. Should we plan their inclusion ? 


Variable | Description | Importance
-- | -- | --
Coldest Minimum Temperature | The very coldest temperature of the year. | Affects safety, recreation,  buildings, transportation, energy use, etc. 
Cooling Degree Days | Annual sum of the number of degrees Celsius that each day’s mean temperature is above 18 °C. | Often used to estimate how much air-conditioning is required in a year. An increase implies hotter or longer summers
Corn Heat Units | Temperature index used by farmers for when it’s ideal to grow corn. | Used by farmers and agricultural researchers to make planting decisions and maximize production.
Date of First Fall Frost | The earliest date in the autumn when temperatures drop below freezing. | Marks the approximate end of the growing season and indicates the return of cold weather.
Date of Last Spring Frost | The latest date in the spring when temperatures drop below freezing. | Marks the approximate beginning of growing season and indicates the return of warm weather.
Freeze-Thaw Cycles | Total number of days per year when temperatures fluctuate between freezing and non-freezing. | Can have major structural impact on infrastructure such as roadways, sidewalks, etc.
Freezing Degree Days | Annual sum of the number of degrees Celsius that each day’s mean temperature is below 0 °C. | High values are associated with relatively cold conditions, and indicate many days with temperatures significantly below freezing.
Frost Days | Total number of days per year when the coldest temperature of the day is less than 0°C. | These conditions indicate frost might form at ground level or on cold services.
Frost-Free Season | Interval between the first frost of the fall and the final frost of the spring. | Indicates when no freezing temperatures occur that might kill or damage plants or animals.
Growing Degree Days (5, 10, 15 °C) | Index of heat available for organisms (such as plants and insects) for growth. | Helpful for crop selection, insect control, and other management.
Heating Degree Days | Annual sum of the number of degrees Celsius a given day's mean temperature is below 18 °C. | Often used as a measure of how much heating is required in a year.
Heavy Precipitation Days (10mm / 20mm) | Total number of days per year when at least 10 mm or 20 mm of rain or frozen precipitation falls. | Heavy rainfall can impact storm drains, flooding, erosion, roads, etc
Icing Days | Total number of days per year when the warmest temperature of the day does not go above freezing (0 °C). | Indicator of the length and/or severity of the winter season, which affects human health and safety, recreation, heating costs, and more.
Maximum Temperature | The highest temperature of the day. | The temperature range we expect within a season or year is a very important aspect of climate. Changes in average and extreme temperatures can dramatically affect our everyday lives as well as a wide range of planning and policy decisions.
Mean Temperature | The average temperature of the day. | idem
Minimum Temperatures | The lowest temperature of the day. | idem
P-Days | Temperature index used by farmers for when it’s ideal to grow potatoes. | Often used by farmers and agricultural researchers to maximize production.
Precipitation | Total annual amount of rain, drizzle, snow, sleet, etc. | Affects water availability, crop viability, hydroelectric power generation, fire risk, snow accumulation, flooding, and drought.
Summer Days | Total number of days per year when the temperature rises to at least 25°C. | Increased health risks to vulnerable people, limits outdoor recreation, increased energy use, increased drought and fire risk, changes to ecosystems and biodiversity, etc. 
Tropical Nights | Total number of days per year when the lowest temperature of the day does not go below 20°C. | idem
Very Cold Days (-30 °C) | Total number of days per year when the temperature drops to -30°C or below. | Indicator of winter severity that affects health, safety, recreation, ecosystems, etc.
Very Hot Days (+30 °C) | Total number of days per year when the temperature rises to at least 30°C. | Increased health risks to vulnerable people, limits outdoor recreation, increased energy use, increased drought and fire risk, changes to ecosystems and biodiversity, etc.
Warmest Maximum Temperature | The highest temperature of the year | idem


Discussion:
coveralls: ## Pull Request Test Coverage Report for [Build 1622](https://coveralls.io/builds/29006985)

* **27** of **75**   **(36.0%)**  changed or added relevant lines in **2** files are covered.
* **6** unchanged lines in **1** file lost coverage.
* Overall coverage decreased (**-1.8%**) to **86.018%**

---

|  Changes Missing Coverage | Covered Lines | Changed/Added Lines | % |
| :-----|--------------|--------|---: |
| [xclim/indices/_threshold.py](https://coveralls.io/builds/29006985/source?filename=xclim%2Findices%2F_threshold.py#L181) | 18 | 66 | 27.27%
<!-- | **Total:** | **27** | **75** | **36.0%** | -->

|  Files with Coverage Reduction | New Missed Lines | % |
| :-----|--------------|--: |
| [xclim/checks.py](https://coveralls.io/builds/29006985/source?filename=xclim%2Fchecks.py#L80) | 6 | 56.44% |
<!-- | **Total:** | **6** |  | -->

|  Totals | [![Coverage Status](https://coveralls.io/builds/29006985/badge)](https://coveralls.io/builds/29006985) |
| :-- | --: |
| Change from base [Build 1621](https://coveralls.io/builds/29002044): |  -1.8% |
| Covered Lines: | 2024 |
| Relevant Lines: | 2353 |

---
##### 💛  - [Coveralls](https://coveralls.io)

Zeitsperre: @aulemahal @huard 

This is good to finally be reviewed. Thanks!
Zeitsperre: @huard we ready to close this?
Zeitsperre: This PR will now also close #391 

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 8db14c96..913c6a60 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -4,6 +4,8 @@ History
 
 0.15.x
 ------
+* Added example notebook for creating a catalog of selected indices
+* Added `growing_season_end`, `last_spring_frost`, `dry_days`,  `hot_spell_frequency`, `hot_spell_max_length`, and `maximum_consecutive_frost_free_days` indices.
 * Dropped use of `fiona.crs` class in lieu of the newer pyproj CRS handler for `subset_shape` operations.
 * Complete internal reorganization of xclim.
 * Internationalization of xclim : add `locales` submodule for localized metadata.
diff --git a/docs/notebooks/index.rst b/docs/notebooks/index.rst
index bdf86990..37c8fd00 100644
--- a/docs/notebooks/index.rst
+++ b/docs/notebooks/index.rst
@@ -8,3 +8,4 @@ Examples
    example
    subset
    ensembles
+   mappings
diff --git a/docs/notebooks/mappings.ipynb b/docs/notebooks/mappings.ipynb
new file mode 100644
index 00000000..a998c357
--- /dev/null
+++ b/docs/notebooks/mappings.ipynb
@@ -0,0 +1,351 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Mapped modules\n",
+    "\n",
+    "`xclim` gives users the ability to generate their own modules as one-to-one mappings with the existing indices library. These mappings can help in emulating existing libraries (such as ICCLIM), with the added benefit of CF-compliant metadata, multilingual metadata support, and optimized calculations using federated resources (using Dask). This can be used for example to tailor existing indices with predefined thresholds without having to rewrite indices.\n",
+    "\n",
+    "Presently, xclim is capable of approximating the indices developed in ICCLIM (https://icclim.readthedocs.io/en/latest/intro.html) and is open to contributions of new indices and library mappings.\n",
+    "\n",
+    "This notebook serves as an example of how one might go about creating their own library of mapped indices."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We'll begin by examining how the mapped modules behave"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# importing libraries\n",
+    "\n",
+    "import xarray as xr\n",
+    "xr.set_options(display_style='html')\n",
+    "\n",
+    "import xclim\n",
+    "from xclim import indices\n",
+    "from xclim import icclim\n",
+    "from xclim import subset\n",
+    "import numpy as np\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = xr.tutorial.open_dataset('air_temperature')\n",
+    "daily_tas = ds.air.resample(time='D').mean(keep_attrs=True)\n",
+    "\n",
+    "fig, axes = plt.subplots(ncols=2, figsize=(12,5))\n",
+    "\n",
+    "daily_tas.sel(lat=45, lon=-71, method=\"nearest\").plot(ax=axes[0], color=\"purple\")\n",
+    "daily_tas.sel(lat=45, lon=-71, method=\"nearest\").plot.hist(ax=axes[1], color=\"green\")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# The following comparisons for Growing Degree Days from xclim.indices and xclim.icclim should be the same.\n",
+    "\n",
+    "gdd_xclim = indices.growing_degree_days(daily_tas, freq=\"MS\").compute()\n",
+    "gdd_icclim = icclim.GD4(daily_tas, freq=\"MS\").compute()\n",
+    "\n",
+    "fig, axes = plt.subplots(ncols=2, figsize=(12,5))\n",
+    "\n",
+    "gdd_icclim.sel(lat=45, lon=-71, method=\"nearest\").plot(ax=axes[0], color=\"cyan\")\n",
+    "gdd_xclim.sel(lat=45, lon=-71, method=\"nearest\").plot(ax=axes[1], color=\"orange\")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Assert that the values are the same between both the original and the mapping of the indice.\n",
+    "\n",
+    "if np.alltrue(gdd_icclim.values == gdd_xclim.values):\n",
+    "    print(\"Datasets hold the same values!\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Let's look more closely at the ICCLIM module\n",
+    "from xclim import icclim\n",
+    "import inspect\n",
+    "\n",
+    "print(icclim.__doc__)\n",
+    "for i in inspect.getmembers(icclim):\n",
+    "    if i[0][0] != '_':\n",
+    "        print(f\"{i[0]} : {i[1]}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Creating your own custom mapped module\n",
+    "\n",
+    "As we can see from the layout of ICCLIM, the module is effectively several wrapped instances of existing indices, with a majority comprising unchanged indices and several composed of indices with different threshold values.\n",
+    "\n",
+    "These mappings can be found in `xclim/xclim/__init__.py` and can be create using the `build_module` function. A mapping requires two types of information:\n",
+    "\n",
+    "* A dictionary of mapped indices\n",
+    "* A docstring, usually written in ReStructuredText format\n",
+    "\n",
+    "Let's try creating our own custom indices for our specific purposes:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from xclim.indicators import build_module\n",
+    "from xclim.core.utils import wrapped_partial\n",
+    "\n",
+    "def awesome_module(mode: str = \"ignore\"):\n",
+    "    mapping = dict(\n",
+    "        egg_cooking_season=wrapped_partial(indices.maximum_consecutive_tx_days, thresh=\"35 degC\"),\n",
+    "        precipitation_above_14mm=wrapped_partial(indices.wetdays, thresh=\"14.0 mm/day\"),\n",
+    "        sweater_weather=indices.tg_min\n",
+    "    )\n",
+    "    mod = build_module(\n",
+    "        \"xclim.awesome\",\n",
+    "        mapping,\n",
+    "        doc=\"\"\"\n",
+    "            =========================\n",
+    "            My Awesome Custom indices\n",
+    "            =========================\n",
+    "            There are only 3 indices that really matter when you come down to brass tacks.\n",
+    "            This mapping library exposes them to users who want to perform real deal \n",
+    "            climate science. Note that none of the checks performed by the `xclim.utils.Indicator` \n",
+    "            class (like with `xclim.atmos` indicators) are performed in this module.\n",
+    "\n",
+    "            \"\"\",\n",
+    "        mode=mode,\n",
+    "    )\n",
+    "    return mod\n",
+    "\n",
+    "AWESOME = awesome_module()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Let's look at our new awesome module\n",
+    "from xclim import awesome\n",
+    "import inspect\n",
+    "\n",
+    "print(awesome.__doc__)\n",
+    "for i in inspect.getmembers(awesome):\n",
+    "    if i[0][0] != '_':\n",
+    "        print(f\"{i[0]} : {i[1]}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sweater_weather_day_temperature = awesome.sweater_weather(daily_tas).compute()\n",
+    "\n",
+    "fig, axes = plt.subplots(ncols=2, figsize=(12,5))\n",
+    "\n",
+    "sweater_weather_day_temperature.isel(time=0).plot(ax=axes[0])\n",
+    "sweater_weather_day_temperature.isel(time=0).plot.hist(ax=axes[1], color=\"red\")\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can also use this to recreate existing indice libraries such as the indicators used by the Prairie Climate Centre's Climate Atlas of Canada (https://climateatlas.ca/variables). \n",
+    "\n",
+    "Here's how we could go about creating a good deal of them:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def __build_pcc(mode: str = \"warn\"):\n",
+    "    \n",
+    "    # The following indices are not implemented here:\n",
+    "        # average_length_of_heat_waves\n",
+    "        # corn_heat_units\n",
+    "        # maximum_temperature\n",
+    "        # minimum_temperature\n",
+    "\n",
+    "    mapping = dict(\n",
+    "        coldest_minimum_temperature=indices.tn_min,\n",
+    "        cooling_degree_days=indices.cooling_degree_days,\n",
+    "        date_of_first_fall_frost=wrapped_partial(\n",
+    "            indices.growing_season_end, thresh=\"0 degC\", mid_date=\"07-15\", window=1\n",
+    "        ),\n",
+    "        date_of_last_spring_frost=indices.last_spring_frost,\n",
+    "        days_above_32_celsius=wrapped_partial(\n",
+    "            indices.tx_days_above, thresh=\"32.0 degC\"\n",
+    "        ),\n",
+    "        days_above_34_celsius=wrapped_partial(\n",
+    "            indices.tx_days_above, thresh=\"34.0 degC\"\n",
+    "        ),\n",
+    "        dry_days=indices.dry_days,\n",
+    "        freeze_thaw_cycles=wrapped_partial(\n",
+    "            indices.daily_freezethaw_cycles,\n",
+    "            thresh_tasmax=\"0 degC\",\n",
+    "            thresh_tasmin=\"-1.0 degC\",\n",
+    "        ),\n",
+    "        freezing_degree_days=wrapped_partial(\n",
+    "            indices.heating_degree_days, thresh=\"0 degC\"\n",
+    "        ),\n",
+    "        frost_days=indices.frost_days,\n",
+    "        frost_free_season=indices.maximum_consecutive_frost_free_days,\n",
+    "        growing_degree_days_base_4_celsius=wrapped_partial(\n",
+    "            indices.growing_degree_days, thresh=\"4.0 degC\"\n",
+    "        ),\n",
+    "        growing_degree_days_base_5_celsius=wrapped_partial(\n",
+    "            indices.growing_degree_days, thresh=\"5.0 degC\"\n",
+    "        ),\n",
+    "        growing_degree_days_base_10_celsius=wrapped_partial(\n",
+    "            indices.growing_degree_days, thresh=\"10.0 degC\"\n",
+    "        ),\n",
+    "        growing_degree_days_base_15_celsius=wrapped_partial(\n",
+    "            indices.growing_degree_days, thresh=\"15.0 degC\"\n",
+    "        ),\n",
+    "        heating_degree_days=wrapped_partial(\n",
+    "            indices.heating_degree_days, thresh=\"18.0 degC\"\n",
+    "        ),\n",
+    "        heavy_precipitation_days_10mm=wrapped_partial(\n",
+    "            indices.wetdays, thresh=\"10.0 mm/day\"\n",
+    "        ),\n",
+    "        heavy_precipitation_days_20mm=wrapped_partial(\n",
+    "            indices.wetdays, thresh=\"20.0 mm/day\"\n",
+    "        ),\n",
+    "        hot_season=wrapped_partial(\n",
+    "            indices.maximum_consecutive_tx_days, thresh=\"30.0 degC\"\n",
+    "        ),\n",
+    "        icing_days=indices.ice_days,\n",
+    "        longest_spell_of_30_celsius_days=indices.hot_spell_max_length,\n",
+    "        max_1_day_precipitation=indices.max_1day_precipitation_amount,\n",
+    "        max_3_day_precipitation=wrapped_partial(\n",
+    "            indices.max_n_day_precipitation_amount, window=3\n",
+    "        ),\n",
+    "        max_5_day_precipitation=wrapped_partial(\n",
+    "            indices.max_n_day_precipitation_amount, window=5\n",
+    "        ),\n",
+    "        mean_temperature=indices.tas,\n",
+    "        mild_winter_days=wrapped_partial(indices.tn_days_below, thresh=\"-5.0 degC\"),\n",
+    "        number_of_heat_waves=indices.hot_spell_frequency,\n",
+    "        summer_days=indices.tx_days_above,\n",
+    "        tropical_nights=wrapped_partial(\n",
+    "            indices.tropical_nights\n",
+    "        ),  \n",
+    "        very_cold_days=wrapped_partial(indices.tn_days_below, thresh=\"-30.0 degC\"),\n",
+    "        very_hot_days=wrapped_partial(indices.tx_days_above, thresh=\"-30.0 degC\"),\n",
+    "        warmest_maximum_temperature=indices.tx_max,\n",
+    "        wet_days=wrapped_partial(indices.wetdays, thresh=\"0.2 mm/day\"),\n",
+    "        winter_days=wrapped_partial(indices.tn_days_below, thresh=\"-15.0 degC\"),\n",
+    "    )\n",
+    "    mod = build_module(\n",
+    "        \"xclim.pcc\",\n",
+    "        mapping,\n",
+    "        doc=\"\"\"\n",
+    "            ==============================\n",
+    "            Prairie Climate Centre indices\n",
+    "            ==============================\n",
+    "            The Prairie Climate Centre (`PCC`_) defines a set of 23 core climate `indices`_.\n",
+    "            They have been made accessible directly in xclim by their climate atlas names.\n",
+    "            However, the methods in this module are only wrappers around the corresponding\n",
+    "            methods of `xclim.indices`. Note that none of the checks performed by\n",
+    "            the `xclim.utils.Indicator` class (like with `xclim.atmos` indicators)\n",
+    "            are performed in this module.\n",
+    "\n",
+    "            .. _PCC: http://prairieclimatecentre.ca/\n",
+    "            .. _indices: https://climateatlas.ca/variables\n",
+    "            \"\"\",\n",
+    "        mode=mode,\n",
+    "    )\n",
+    "    return mod\n",
+    "\n",
+    "PCC = __build_pcc()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from xclim import pcc\n",
+    "import inspect\n",
+    "\n",
+    "print(pcc.__doc__)\n",
+    "for i in inspect.getmembers(pcc):\n",
+    "    if i[0][0] != '_':\n",
+    "        print(f\"{i[0]} : {i[1]}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Let's see what the documentation will look like\n",
+    "\n",
+    "help(pcc)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/notebooks/subset.ipynb b/docs/notebooks/subset.ipynb
index f29b7311..25801a32 100644
--- a/docs/notebooks/subset.ipynb
+++ b/docs/notebooks/subset.ipynb
@@ -593,7 +593,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.5"
+   "version": "3.7.3"
   },
   "nbsphinx": {
    "execute": "never"
diff --git a/docs/notebooks/usage.ipynb b/docs/notebooks/usage.ipynb
index d8670592..2d53ea97 100644
--- a/docs/notebooks/usage.ipynb
+++ b/docs/notebooks/usage.ipynb
@@ -208,7 +208,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.5"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,
diff --git a/setup.cfg b/setup.cfg
index eaaedc74..3ee31e5c 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.14.3-beta
+current_version = 0.14.4-beta
 commit = False
 tag = False
 parse = (?P<major>\d+)\.(?P<minor>\d+).(?P<patch>\d+)(\-(?P<release>[a-z]+))?
diff --git a/setup.py b/setup.py
index 7f316759..a45451eb 100644
--- a/setup.py
+++ b/setup.py
@@ -12,7 +12,7 @@ URL = "https://github.com/Ouranosinc/xclim"
 AUTHOR = "Travis Logan"
 AUTHOR_EMAIL = "logan.travis@ouranos.ca"
 REQUIRES_PYTHON = ">=3.6.0"
-VERSION = "0.14.3-beta"
+VERSION = "0.14.4-beta"
 LICENSE = "Apache Software License 2.0"
 
 with open("README.rst") as readme_file:
@@ -38,6 +38,7 @@ requirements = [
     "bottleneck~=1.3.1",
     "boltons>=19.0",
     "scikit-learn>=0.21.3",
+    "packaging>=20.0",
 ]
 
 if not os.getenv("READTHEDOCS"):
diff --git a/xclim/__init__.py b/xclim/__init__.py
index a7aa01a5..d58d9817 100644
--- a/xclim/__init__.py
+++ b/xclim/__init__.py
@@ -8,4 +8,4 @@ from xclim.indicators import seaIce
 
 __author__ = """Travis Logan"""
 __email__ = "logan.travis@ouranos.ca"
-__version__ = "0.14.3-beta"
+__version__ = "0.14.4-beta"
diff --git a/xclim/indicators/atmos/_precip.py b/xclim/indicators/atmos/_precip.py
index 23b2bc6f..e1febb4b 100644
--- a/xclim/indicators/atmos/_precip.py
+++ b/xclim/indicators/atmos/_precip.py
@@ -9,6 +9,7 @@ __all__ = [
     "max_1day_precipitation_amount",
     "max_n_day_precipitation_amount",
     "wetdays",
+    "dry_days",
     "maximum_consecutive_dry_days",
     "maximum_consecutive_wet_days",
     "daily_pr_intensity",
@@ -73,6 +74,16 @@ wetdays = Pr(
     compute=indices.wetdays,
 )
 
+dry_days = Pr(
+    identifier="dry_days",
+    units="days",
+    standard_name="number_of_days_with_lwe_thickness_of_precipitation_amount_below_threshold",
+    long_name="Number of dry days (precip < {thresh})",
+    description="{freq} number of days with daily precipitation under {thresh}.",
+    cell_methods="time: sum within days time: sum over days",
+    compute=indices.dry_days,
+)
+
 maximum_consecutive_wet_days = Pr(
     identifier="cwd",
     units="days",
diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py
index 3c6a1f01..527d04e0 100644
--- a/xclim/indicators/atmos/_temperature.py
+++ b/xclim/indicators/atmos/_temperature.py
@@ -14,6 +14,8 @@ __all__ = [
     "heat_wave_max_length",
     "heat_wave_total_length",
     "heat_wave_index",
+    "hot_spell_frequency",
+    "hot_spell_max_length",
     "tg",
     "tg_mean",
     "tg10p",
@@ -39,9 +41,12 @@ __all__ = [
     "growing_degree_days",
     "freshet_start",
     "frost_days",
+    "last_spring_frost",
     "ice_days",
     "consecutive_frost_days",
+    "maximum_consecutive_frost_free_days",
     "growing_season_length",
+    "growing_season_end",
     "tropical_nights",
 ]
 
@@ -124,8 +129,8 @@ tx_tn_days_above = TasminTasmax(
     units="days",
     standard_name="number_of_days_with_air_temperature_above_threshold",
     long_name="Number of days with Tmax > {thresh_tasmax} and Tmin > {thresh_tasmin}",
-    description="{freq} number of days where daily maximum temperature exceeds"
-    " {thresh_tasmax} and minimum temperature exceeds {thresh_tasmin}",
+    description="{freq} number of days where daily maximum temperature exceeds "
+    "{thresh_tasmax} and minimum temperature exceeds {thresh_tasmin}",
     cell_methods="",
     compute=indices.tx_tn_days_above,
 )
@@ -134,7 +139,7 @@ heat_wave_frequency = TasminTasmax(
     identifier="heat_wave_frequency",
     units="",
     standard_name="heat_wave_events",
-    long_name="Number of heat wave events (Tmin > {thresh_tasmin}"
+    long_name="Number of heat wave events (Tmin > {thresh_tasmin} "
     "and Tmax > {thresh_tasmax} for >= {window} days)",
     description="{freq} number of heat wave events over a given period. "
     "An event occurs when the minimum and maximum daily "
@@ -152,7 +157,7 @@ heat_wave_max_length = TasminTasmax(
     standard_name="spell_length_of_days_with_air_temperature_above_threshold",
     long_name="Maximum length of heat wave events (Tmin > {thresh_tasmin}"
     "and Tmax > {thresh_tasmax} for >= {window} days)",
-    description="{freq} maximum length of heat wave events occuring in a given period."
+    description="{freq} maximum length of heat wave events occurring in a given period. "
     "An event occurs when the minimum and maximum daily "
     "temperature both exceeds specific thresholds "
     "(Tmin > {thresh_tasmin} and Tmax > {thresh_tasmax}) over "
@@ -166,9 +171,9 @@ heat_wave_total_length = TasminTasmax(
     identifier="heat_wave_total_length",
     units="days",
     standard_name="spell_length_of_days_with_air_temperature_above_threshold",
-    long_name="Total length of heat wave events (Tmin > {thresh_tasmin}"
+    long_name="Total length of heat wave events (Tmin > {thresh_tasmin} "
     "and Tmax > {thresh_tasmax} for >= {window} days)",
-    description="{freq} total length of heat wave events occuring in a given period."
+    description="{freq} total length of heat wave events occurring in a given period. "
     "An event occurs when the minimum and maximum daily "
     "temperature both exceeds specific thresholds "
     "(Tmin > {thresh_tasmin} and Tmax > {thresh_tasmax}) over "
@@ -190,6 +195,32 @@ heat_wave_index = Tasmax(
     compute=indices.heat_wave_index,
 )
 
+hot_spell_frequency = Tasmax(
+    identifier="hot_spell_frequency",
+    units="",
+    standard_name="hot_spell_events",
+    long_name="Number of hot spell events (Tmax > {thresh_tasmax} for >= {window} days)",
+    description="{freq} number of hot spell events over a given period. "
+    "An event occurs when the maximum daily temperature exceeds a specific threshold: (Tmax > {thresh_tasmax}) "
+    "over a minimum number of days ({window}).",
+    cell_methods="",
+    keywords="health,",
+    compute=indices.hot_spell_frequency,
+)
+
+hot_spell_max_length = Tasmax(
+    identifier="heat_wave_max_length",
+    units="days",
+    standard_name="spell_length_of_days_with_air_temperature_above_threshold",
+    long_name="Maximum length of hot spell events (Tmax > {thresh_tasmax} for >= {window} days)",
+    description="{freq} maximum length of hot spell events occurring in a given period. "
+    "An event occurs when the maximum daily temperature exceeds a specific threshold: (Tmax > {thresh_tasmax}) "
+    "over a minimum number of days ({window}).",
+    cell_methods="",
+    keywords="health,",
+    compute=indices.hot_spell_max_length,
+)
+
 tg = TasminTasmax(
     identifier="tg",
     units="K",
@@ -285,8 +316,8 @@ daily_temperature_range_variability = TasminTasmax(
     units="K",
     standard_name="air_temperature",
     long_name="Mean Diurnal Temperature Range Variability",
-    description="{freq} mean diurnal temparature range variability ("
-    "defined as the average day-to-day variation "
+    description="{freq} mean diurnal temparature range variability "
+    "(defined as the average day-to-day variation "
     "in daily temperature range "
     "for the given time period)",
     cell_methods="time range within days time: difference "
@@ -311,8 +342,8 @@ cold_spell_duration_index = Tasmin(
     standard_name="cold_spell_duration_index",
     long_name="Cold Spell Duration Index, count of days with at "
     "least {window} consecutive days when Tmin < 10th percentile",
-    description="{freq} number of days with at least {window} consecutive days"
-    " where the daily minimum temperature is below the 10th "
+    description="{freq} number of days with at least {window} consecutive days "
+    "where the daily minimum temperature is below the 10th "
     "percentile. The 10th percentile should be computed for "
     "a 5-day window centred on each calendar day in the  1961-1990 period",
     cell_methods="",
@@ -387,17 +418,27 @@ frost_days = Tasmin(
     units="days",
     standard_name="days_with_air_temperature_below_threshold",
     long_name="Number of Frost Days (Tmin < 0C)",
-    description="{freq} number of days with minimum daily " "temperature below 0℃.",
+    description="{freq} number of days with minimum daily temperature below 0℃.",
     cell_methods="time: minimum within days time: sum over days",
     compute=indices.frost_days,
 )
 
+last_spring_frost = Tasmin(
+    identifier="last_spring_frost",
+    units="",
+    standard_name="day_of_year",
+    long_name="Day of year of last spring frost",
+    description="Day of year of last spring frost, defined as the last day a minimum temperature "
+    "threshold of {thresh} is not exceeded before a given date.",
+    compute=indices.last_spring_frost,
+)
+
 ice_days = Tasmax(
     identifier="ice_days",
     standard_name="days_with_air_temperature_below_threshold",
     units="days",
     long_name="Number of Ice Days (Tmax < 0℃)",
-    description="{freq} number of days with maximum daily " "temperature below 0℃",
+    description="{freq} number of days with maximum daily temperature below 0℃",
     cell_methods="time: maximum within days time: sum over days",
     compute=indices.ice_days,
 )
@@ -413,6 +454,17 @@ consecutive_frost_days = Tasmin(
     compute=indices.consecutive_frost_days,
 )
 
+maximum_consecutive_frost_free_days = Tasmin(
+    identifier="consecutive_frost_free_days",
+    units="days",
+    standard_name="spell_length_of_days_with_air_temperature_above_threshold",
+    long_name="Maximum number of consecutive days with Tmin > {thresh}",
+    description="{freq} maximum number of consecutive days with "
+    "minimum daily temperature above {thresh}.",
+    cell_methods="time: min within days time: maximum over days",
+    compute=indices.maximum_consecutive_frost_free_days,
+)
+
 growing_season_length = Tas(
     identifier="growing_season_length",
     units="days",
@@ -426,6 +478,18 @@ growing_season_length = Tas(
     compute=indices.growing_season_length,
 )
 
+growing_season_end = Tas(
+    identifier="growing_season_end",
+    units="",
+    standard_name="day_of_year",
+    long_name="Day of year of growing season end",
+    description="Day of year of end of growing season, defined as the first day of "
+    "consistent inferior threshold temperature of {thresh} after a run of "
+    "{window} days superior to threshold temperature.",
+    cell_methods="",
+    compute=indices.growing_season_end,
+)
+
 tropical_nights = Tasmin(
     identifier="tropical_nights",
     units="days",
diff --git a/xclim/indices/_multivariate.py b/xclim/indices/_multivariate.py
index 4e05f546..d2a252ca 100644
--- a/xclim/indices/_multivariate.py
+++ b/xclim/indices/_multivariate.py
@@ -445,7 +445,8 @@ def drought_code(
 ):
     r"""The daily drought code (FWI component)
 
-    The drought code is part of the Canadian Forest Fire Weather Index System. It is a numeric rating of the average moisture content of organic layers.
+    The drought code is part of the Canadian Forest Fire Weather Index System.
+    It is a numeric rating of the average moisture content of organic layers.
 
     Parameters
     ----------
@@ -461,6 +462,8 @@ def drought_code(
       Initial values of the drought code.
     start_date : str, datetime.datetime
       Date at which to start the computation, dc0/dmc0/ffcm0 should be given at the day before.
+    start_up_mode : str
+      #TODO: Add a docstring for this signature.
     params :
       Any other keyword parameters as defined in `xclim.indices.fwi.fire_weather_ufunc`.
 
@@ -1346,7 +1349,7 @@ def warm_spell_duration_index(
     tx90 : float
       90th percentile of daily maximum temperature [℃] or [K]
     window : int
-      Minimum number of days with temperature below threshold to qualify as a warm spell.
+      Minimum number of days with temperature above threshold to qualify as a warm spell.
     freq : str
       Resampling frequency; Defaults to "YS".
 
diff --git a/xclim/indices/_threshold.py b/xclim/indices/_threshold.py
index eda2fa0c..2e3b1be6 100644
--- a/xclim/indices/_threshold.py
+++ b/xclim/indices/_threshold.py
@@ -22,20 +22,26 @@ xarray.set_options(enable_cftimeindex=True)  # Set xarray to use cftimeindex
 __all__ = [
     "cold_spell_days",
     "daily_pr_intensity",
-    "maximum_consecutive_wet_days",
     "cooling_degree_days",
     "freshet_start",
     "growing_degree_days",
+    "growing_season_end",
     "growing_season_length",
+    "last_spring_frost",
     "heat_wave_index",
     "heating_degree_days",
+    "hot_spell_frequency",
+    "hot_spell_max_length",
     "tn_days_below",
     "tx_days_above",
     "warm_day_frequency",
     "warm_night_frequency",
     "wetdays",
+    "dry_days",
     "maximum_consecutive_dry_days",
+    "maximum_consecutive_frost_free_days",
     "maximum_consecutive_tx_days",
+    "maximum_consecutive_wet_days",
     "sea_ice_area",
     "sea_ice_extent",
     "tropical_nights",
@@ -43,7 +49,9 @@ __all__ = [
 
 
 @declare_units("days", tas="[temperature]", thresh="[temperature]")
-def cold_spell_days(tas, thresh="-10 degC", window: int = 5, freq="AS-JUL"):
+def cold_spell_days(
+    tas, thresh: str = "-10 degC", window: int = 5, freq: str = "AS-JUL"
+):
     r"""Cold spell days
 
     The number of days that are part of a cold spell, defined as five or more consecutive days with mean daily
@@ -85,7 +93,7 @@ def cold_spell_days(tas, thresh="-10 degC", window: int = 5, freq="AS-JUL"):
 
 
 @declare_units("mm/day", pr="[precipitation]", thresh="[precipitation]")
-def daily_pr_intensity(pr, thresh="1 mm/day", freq="YS"):
+def daily_pr_intensity(pr, thresh: str = "1 mm/day", freq: str = "YS"):
     r"""Average daily precipitation intensity
 
     Return the average precipitation over wet days.
@@ -142,6 +150,39 @@ def daily_pr_intensity(pr, thresh="1 mm/day", freq="YS"):
     return sd / wd
 
 
+@declare_units("days", pr="[precipitation]", thresh="[precipitation]")
+def dry_days(pr: xarray.DataArray, thresh: str = "0.2 mm/d", freq: str = "YS"):
+    r"""Dry days
+
+    The number of days with daily precipitation below threshold.
+
+    Parameters
+    ----------
+    pr : xarray.DataArray
+      Daily precipitation [mm/d or kg/m²/s]
+    thresh : str
+      Threshold temperature on which to base evaluation. Default: '0.2 mm/d'.
+    freq : str
+      Resampling frequency; Defaults to "YS".
+
+    Returns
+    -------
+    xarray.DataArray
+      Number of days with daily precipitation below threshold.
+
+    Notes
+    -----
+    Let :math:`PR_{ij}` be the daily precipitation at day :math:`i` of period :math:`j`. Then
+    counted is the number of days where:
+
+    .. math::
+
+        \sum PR_{ij} < Threshold [mm/day]
+    """
+    thresh = convert_units_to(thresh, pr)
+    return pr.pipe(lambda x: (pr < thresh) * 1).resample(time=freq).sum(dim="time")
+
+
 @declare_units("days", pr="[precipitation]", thresh="[precipitation]")
 def maximum_consecutive_wet_days(
     pr: xarray.DataArray, thresh: str = "1 mm/day", freq: str = "YS"
@@ -237,18 +278,18 @@ def freshet_start(
     Parameters
     ----------
     tas : xarray.DataArray
-      Mean daily temperature [℃] or [K]
+      Mean daily temperature [℃] or [K].
     thresh : str
-      Threshold temperature on which to base evaluation [℃] or [K]. Default '0 degC'
+      Threshold temperature on which to base evaluation [℃] or [K]. Default '0 degC'.
     window : int
-      Minimum number of days with temperature above threshold needed for evaluation
+      Minimum number of days with temperature above threshold needed for evaluation.
     freq : str
       Resampling frequency; Defaults to "YS".
 
     Returns
     -------
     xarray.DataArray
-      Day of the year when temperature exceeds threshold over a given number of days for the first time. If there are
+      Day of the year when temperature exceeds threshold over a given number of days for the first time. If there is
       no such day, return np.nan.
 
     Notes
@@ -260,7 +301,7 @@ def freshet_start(
 
        \prod_{j=i}^{i+w} [x_j > thresh]
 
-    is true, where :math:`w` is the number of days the temperature threshold should be exceeded,  and :math:`[P]` is
+    is true, where :math:`w` is the number of days the temperature threshold should be exceeded, and :math:`[P]` is
     1 if :math:`P` is true, and 0 if false.
     """
     thresh = convert_units_to(thresh, tas)
@@ -281,16 +322,16 @@ def growing_degree_days(
     Parameters
     ---------
     tas : xarray.DataArray
-      Mean daily temperature [℃] or [K]
+      Mean daily temperature [℃] or [K].
     thresh : str
       Threshold temperature on which to base evaluation [℃] or [K]. Default: '4.0 degC'.
     freq : str
-      Resampling frequency; Defaults to "YS".
+      Resampling frequency. Default: "YS".
 
     Returns
     -------
     xarray.DataArray
-      The sum of growing degree-days above 4℃
+      The sum of growing degree-days above a given threshold.
 
     Notes
     -----
@@ -307,6 +348,51 @@ def growing_degree_days(
     )
 
 
+@declare_units("", tas="[temperature]", thresh="[temperature]")
+def growing_season_end(
+    tas: xarray.DataArray,
+    thresh: str = "5.0 degC",
+    mid_date: str = "07-01",
+    window: int = 5,
+    freq: str = "YS",
+):
+    r"""Day of the year of the start of a sequence of days with a temperature consistently below a threshold, after a period with temperatures consistently above the same threshold.
+
+    Returns the first day of period where a temperature is inferior to a threshold
+    over a given run of days.
+
+    Parameters
+    ----------
+    tas : xarray.DataArray
+      Mean daily temperature [℃] or [K].
+    thresh : str
+      Threshold temperature on which to base evaluation [℃] or [K]. Default '5.0 degC'.
+    mid_date : str
+      Date of the year after which to look for the end of the season. Should have the format '%m-%d'.
+    window : int
+      Minimum number of days with temperature below threshold needed for evaluation.
+    freq : str
+      Resampling frequency. Default: "YS".
+
+    Returns
+    -------
+    xarray.DataArray
+      Day of the year when temperature is inferior to a threshold over a given number of days for the first time.
+      If there is no such day or if a growing season is not detected, returns np.nan.
+      If the growing season does not end within the time period, returns the last day of the period.
+    """
+    thresh = convert_units_to(thresh, tas)
+    cond = tas >= thresh
+
+    return cond.resample(time=freq).map(
+        rl.run_end_after_date,
+        window=window,
+        date=mid_date,
+        dim="time",
+        coord="dayofyear",
+    )
+
+
 @declare_units("days", tas="[temperature]", thresh="[temperature]")
 def growing_season_length(
     tas: xarray.DataArray,
@@ -318,17 +404,17 @@ def growing_season_length(
     r"""Growing season length.
 
     The number of days between the first occurrence of at least
-    six consecutive days with mean daily temperature over 5℃ and
-    the first occurrence of at least six consecutive days with
-    mean daily temperature below 5℃ after a certain date. (Usually
-    July 1st in the northern hemisphere and January 1st in the southern hemisphere.)
+    six consecutive days with mean daily temperature over a threshold (default: 5℃) and
+    the first occurrence of at least six consecutive days with mean daily temperature
+    below the same threshold after a certain date. (Usually July 1st in the northern
+    hemisphere and January 1st in the southern hemisphere.)
 
-    WARNING: The default values are only valid for the northern hemisphere.
+    WARNING: The default calendar values are only valid for the northern hemisphere.
 
     Parameters
     ---------
     tas : xarray.DataArray
-      Mean daily temperature [℃] or [K]
+      Mean daily temperature [℃] or [K].
     thresh : str
       Threshold temperature on which to base evaluation [℃] or [K]. Default: '5.0 degC'.
     window : int
@@ -336,7 +422,7 @@ def growing_season_length(
     mid_date : str
       Date of the year after which to look for the end of the season. Should have the format '%m-%d'.
     freq : str
-      Resampling frequency; Defaults to "YS".
+      Resampling frequency. Default: "YS".
 
     Returns
     -------
@@ -365,39 +451,55 @@ def growing_season_length(
     >>> gsl = growing_season_length(tas, mid_date='01-01', freq='AS-Jul')
     """
     thresh = convert_units_to(thresh, tas)
+    cond = tas >= thresh
+
+    return cond.resample(time=freq).map(
+        rl.run_length_with_date, window=window, date=mid_date, dim="time",
+    )
+
+
+@declare_units("days", tas="[temperature]", thresh="[temperature]")
+def last_spring_frost(
+    tas: xarray.DataArray,
+    thresh: str = "0 degC",
+    before_date: str = "07-01",
+    window: int = 1,
+    freq: str = "YS",
+):
+    r"""Last day of temperatures inferior to a threshold temperature.
 
-    mid_doy = datetime.datetime.strptime(mid_date, "%m-%d").timetuple().tm_yday
-
-    def compute_gsl(yrdata):
-        if (
-            yrdata.chunks is not None
-            and len(yrdata.chunks[yrdata.dims.index("time")]) > 1
-        ):
-            yrdata = yrdata.chunk({"time": -1})
-        mid_idx = np.where(yrdata.time.dt.dayofyear == mid_doy)[0]
-        if (
-            mid_idx.size == 0
-        ):  # The mid date is not in the group. Happens at boundaries.
-            allNans = xarray.full_like(yrdata.isel(time=0), np.nan)
-            allNans.attrs = {}
-            return allNans
-        end = rl.first_run(
-            yrdata.where(yrdata.time >= yrdata.time[mid_idx][0]) < thresh,
-            window,
-            "time",
-        )
-        beg = rl.first_run(yrdata > thresh, window, "time")
-        sl = end - beg
-        sl = xarray.where(
-            beg.isnull() & end.notnull(), 0, sl
-        )  # If everything is under thresh
-        sl = xarray.where(
-            beg.notnull() & end.isnull(), yrdata.time.size - beg, sl
-        )  # If gs is not ended at end of year
-        return sl.where(sl >= 0)  # When end happens before beg.
-
-    gsl = tas.resample(time=freq).map(compute_gsl)
-    return gsl
+    Returns last day of period where a temperature is inferior to a threshold
+    over a given number of days and limited to a final calendar date.
+
+    Parameters
+    ----------
+    tas : xarray.DataArray
+      Mean daily temperature [℃] or [K].
+    thresh : str
+      Threshold temperature on which to base evaluation [℃] or [K]. Default '0 degC'.
+    before_date : str
+      Date of the year before which to look for the final frost event. Should have the format '%m-%d'.
+    window : int
+      Minimum number of days with temperature below threshold needed for evaluation.
+    freq : str
+      Resampling frequency; Defaults to "YS".
+
+    Returns
+    -------
+    xarray.DataArray
+      Day of the year when temperature is inferior to a threshold over a given number of days for the first time.
+      If there is no such day, return np.nan.
+    """
+    thresh = convert_units_to(thresh, tas)
+    cond = tas < thresh
+
+    return cond.resample(time=freq).map(
+        rl.last_run_before_date,
+        window=window,
+        date=before_date,
+        dim="time",
+        coord="dayofyear",
+    )
 
 
 @declare_units("days", tasmax="[temperature]", thresh="[temperature]")
@@ -470,11 +572,128 @@ def heating_degree_days(
     return tas.pipe(lambda x: thresh - x).clip(0).resample(time=freq).sum(dim="time")
 
 
+@declare_units(
+    "days", tasmax="[temperature]", thresh_tasmax="[temperature]",
+)
+def hot_spell_max_length(
+    tasmax: xarray.DataArray,
+    thresh_tasmax: str = "30 degC",
+    window: int = 1,
+    freq: str = "YS",
+) -> xarray.DataArray:
+    # Dev note : we should decide if it is deg K or C
+    r"""Longest hot spell
+
+    Longest spell of high temperatures over a given period.
+
+    The longest series of consecutive days with tasmax ≥ 30 °C. Here, there is no minimum threshold for number of
+    days in a row that must be reached or exceeded to count as a spell. A year with zero +30 °C days will return a
+    longest spell value of zero.
+
+    Parameters
+    ----------
+    tasmax : xarray.DataArray
+      Maximum daily temperature [℃] or [K]
+    thresh_tasmax : str
+      The maximum temperature threshold needed to trigger a heatwave event [℃] or [K]. Default : '30 degC'
+    window : int
+      Minimum number of days with temperatures above thresholds to qualify as a heatwave.
+    freq : str
+      Resampling frequency; Defaults to "YS".
+
+    Returns
+    -------
+    xarray.DataArray
+      Maximum length of continuous hot days at the wanted frequency
+
+    Notes
+    -----
+    The thresholds of 22° and 25°C for night temperatures and 30° and 35°C for day temperatures were selected by
+    Health Canada professionals, following a temperature–mortality analysis. These absolute temperature thresholds
+    characterize the occurrence of hot weather events that can result in adverse health outcomes for Canadian
+    communities (Casati et al., 2013).
+
+    In Robinson (2001), the parameters would be `thresh_tasmin=27.22, thresh_tasmax=39.44, window=2` (81F, 103F).
+
+    References
+    ----------
+    Casati, B., A. Yagouti, and D. Chaumont, 2013: Regional Climate Projections of Extreme Heat Events in Nine Pilot
+    Canadian Communities for Public Health Planning. J. Appl. Meteor. Climatol., 52, 2669–2698,
+    https://doi.org/10.1175/JAMC-D-12-0341.1
+
+    Robinson, P.J., 2001: On the Definition of a Heat Wave. J. Appl. Meteor., 40, 762–775,
+    https://doi.org/10.1175/1520-0450(2001)040<0762:OTDOAH>2.0.CO;2
+    """
+    thresh_tasmax = convert_units_to(thresh_tasmax, tasmax)
+
+    cond = tasmax > thresh_tasmax
+    group = cond.resample(time=freq)
+    max_l = group.apply(rl.longest_run, dim="time")
+    return max_l.where(max_l >= window, 0)
+
+
+@declare_units(
+    "", tasmax="[temperature]", thresh_tasmax="[temperature]",
+)
+def hot_spell_frequency(
+    tasmax: xarray.DataArray,
+    thresh_tasmax: str = "30 degC",
+    window: int = 3,
+    freq: str = "YS",
+) -> xarray.DataArray:
+    # Dev note : we should decide if it is deg K or C
+    r"""Hot spell frequency
+
+    Number of hot spells over a given period. A hot spell is defined as an event
+    where the maximum daily temperature exceeds a specific threshold
+    over a minimum number of days.
+
+    Parameters
+    ----------
+    tasmax : xarray.DataArray
+      Maximum daily temperature [℃] or [K]
+    thresh_tasmax : str
+      The maximum temperature threshold needed to trigger a heatwave event [℃] or [K]. Default : '30 degC'
+    window : int
+      Minimum number of days with temperatures above thresholds to qualify as a heatwave.
+    freq : str
+      Resampling frequency; Defaults to "YS".
+
+    Returns
+    -------
+    xarray.DataArray
+      Number of heatwave at the wanted frequency
+
+    Notes
+    -----
+    The thresholds of 22° and 25°C for night temperatures and 30° and 35°C for day temperatures were selected by
+    Health Canada professionals, following a temperature–mortality analysis. These absolute temperature thresholds
+    characterize the occurrence of hot weather events that can result in adverse health outcomes for Canadian
+    communities (Casati et al., 2013).
+
+    In Robinson (2001), the parameters would be `thresh_tasmin=27.22, thresh_tasmax=39.44, window=2` (81F, 103F).
+
+    References
+    ----------
+    Casati, B., A. Yagouti, and D. Chaumont, 2013: Regional Climate Projections of Extreme Heat Events in Nine Pilot
+    Canadian Communities for Public Health Planning. J. Appl. Meteor. Climatol., 52, 2669–2698,
+    https://doi.org/10.1175/JAMC-D-12-0341.1
+
+    Robinson, P.J., 2001: On the Definition of a Heat Wave. J. Appl. Meteor., 40, 762–775,
+    https://doi.org/10.1175/1520-0450(2001)040<0762:OTDOAH>2.0.CO;2
+    """
+    thresh_tasmax = convert_units_to(thresh_tasmax, tasmax)
+
+    cond = tasmax > thresh_tasmax
+    group = cond.resample(time=freq)
+    return group.apply(rl.windowed_run_events, window=window, dim="time")
+
+
 @declare_units("days", tasmin="[temperature]", thresh="[temperature]")
 def tn_days_below(
     tasmin: xarray.DataArray, thresh: str = "-10.0 degC", freq: str = "YS"
 ):
-    r"""Number of days with tmin below a threshold in
+    r"""Number of days with tmin below a threshold
 
     Number of days where daily minimum temperature is below a threshold.
 
@@ -685,13 +904,57 @@ def maximum_consecutive_dry_days(
     return group.apply(rl.longest_run, dim="time")
 
 
+@declare_units("days", tasmin="[temperature]", thresh="[temperature]")
+def maximum_consecutive_frost_free_days(
+    tasmin: xarray.DataArray, thresh: str = "0 degC", freq: str = "YS"
+):
+    r"""Maximum number of consecutive frost free days (Tn > 0℃)
+
+    Return the maximum number of consecutive days within the period where the
+    minimum temperature is above a certain threshold.
+
+    Parameters
+    ----------
+    tasmin : xarray.DataArray
+      Max daily temperature [K]
+    thresh : str
+      Threshold temperature [K].
+    freq : str
+      Resampling frequency; Defaults to "YS".
+
+    Returns
+    -------
+    xarray.DataArray
+      The maximum number of consecutive frost free days.
+
+    Notes
+    -----
+    Let :math:`\mathbf{t}=t_0, t_1, \ldots, t_n` be a daily minimum temperature series and :math:`thresh` the threshold
+    above which a day is considered a frost free day. Let :math:`\mathbf{s}` be the sorted vector of indices :math:`i`
+    where :math:`[t_i < thresh] \neq [t_{i+1} < thresh]`, that is, the days when the temperature crosses the threshold.
+    Then the maximum number of consecutive frost free days is given by
+
+    .. math::
+
+       \max(\mathbf{d}) \quad \mathrm{where} \quad d_j = (s_j - s_{j-1}) [t_{s_j} > thresh]
+
+    where :math:`[P]` is 1 if :math:`P` is true, and 0 if false. Note that this formula does not handle sequences at
+    the start and end of the series, but the numerical algorithm does.
+    """
+    t = convert_units_to(thresh, tasmin)
+    group = (tasmin > t).resample(time=freq)
+
+    return group.apply(rl.longest_run, dim="time")
+
+
 @declare_units("days", tasmax="[temperature]", thresh="[temperature]")
 def maximum_consecutive_tx_days(
     tasmax: xarray.DataArray, thresh: str = "25 degC", freq: str = "YS"
 ):
     r"""Maximum number of consecutive summer days (Tx > 25℃)
 
-    Return the maximum number of consecutive days within the period where temperature is above a certain threshold.
+    Return the maximum number of consecutive days within the period where
+    the maximum temperature is above a certain threshold.
 
     Parameters
     ----------
@@ -728,7 +991,7 @@ def maximum_consecutive_tx_days(
 
 
 @declare_units("[area]", sic="[]", area="[area]", thresh="[]")
-def sea_ice_area(sic, area, thresh="15 pct"):
+def sea_ice_area(sic, area, thresh: str = "15 pct"):
     """Return the total sea ice area.
 
     Sea ice area measures the total sea ice covered area where sea ice concentration is above a threshold,
@@ -765,7 +1028,7 @@ def sea_ice_area(sic, area, thresh="15 pct"):
 
 
 @declare_units("[area]", sic="[]", area="[area]", thresh="[]")
-def sea_ice_extent(sic, area, thresh="15 pct"):
+def sea_ice_extent(sic, area, thresh: str = "15 pct"):
     """Return the total sea ice extent.
 
     Sea ice extent measures the *ice-covered* area, where a region is considered ice-covered if its sea ice
@@ -801,7 +1064,7 @@ def sea_ice_extent(sic, area, thresh="15 pct"):
 
 @declare_units("days", tasmin="[temperature]", thresh="[temperature]")
 def tropical_nights(
-    tasmin: xarray.DataArray, thresh: str = "20.0 degC", freq: str = "YS"
+    tasmin: xarray.DataArray, thresh: str = "20.0 degC", freq: str = "YS",
 ):
     r"""Tropical nights
 
diff --git a/xclim/indices/run_length.py b/xclim/indices/run_length.py
index c769fc7d..e5c9b55e 100644
--- a/xclim/indices/run_length.py
+++ b/xclim/indices/run_length.py
@@ -6,6 +6,7 @@ Run length algorithms submodule
 Computation of statistics on runs of True values in boolean arrays.
 """
 import logging
+from datetime import datetime
 from functools import partial
 from typing import Optional
 from typing import Sequence
@@ -25,17 +26,16 @@ npts_opt = 9000
 def get_npts(da: xr.DataArray) -> int:
     """Return the number of gridpoints in a DataArray.
 
-        Parameters
-        ----------
-        da : xarray.DataArray
-          N-dimensional input array
-
-        Returns
-        -------
-        int
-          Product of input DataArray coordinate sizes excluding the dimension 'time'
-        """
+    Parameters
+    ----------
+    da : xarray.DataArray
+      N-dimensional input array
 
+    Returns
+    -------
+    int
+      Product of input DataArray coordinate sizes excluding the dimension 'time'
+    """
     coords = list(da.coords)
     coords.remove("time")
     npts = 1
@@ -86,21 +86,21 @@ def longest_run(
 ):
     """Return the length of the longest consecutive run of True values.
 
-        Parameters
-        ----------
-        da : xr.DataArray
-          N-dimensional array (boolean)
-        dim : str
-          Dimension along which to calculate consecutive run; Default: 'time'.
-        ufunc_1dim : Union[str, bool]
-          Use the 1d 'ufunc' version of this function : default (auto) will attempt to select optimal
-          usage based on number of data points.  Using 1D_ufunc=True is typically more efficient
-          for dataarray with a small number of gridpoints.
-        Returns
-        -------
-        N-dimensional array (int)
-          Length of longest run of True values along dimension
-        """
+    Parameters
+    ----------
+    da : xr.DataArray
+      N-dimensional array (boolean)
+    dim : str
+      Dimension along which to calculate consecutive run; Default: 'time'.
+    ufunc_1dim : Union[str, bool]
+      Use the 1d 'ufunc' version of this function : default (auto) will attempt to select optimal
+      usage based on number of data points.  Using 1D_ufunc=True is typically more efficient
+      for dataarray with a small number of gridpoints.
+    Returns
+    -------
+    N-dimensional array (int)
+      Length of longest run of True values along dimension
+    """
     if ufunc_1dim == "auto":
         npts = get_npts(da)
         ufunc_1dim = npts <= npts_opt
@@ -122,23 +122,23 @@ def windowed_run_events(
 ) -> xr.DataArray:
     """Return the number of runs of a minimum length.
 
-        Parameters
-        ----------
-        da: xr.DataArray
-          Input N-dimensional DataArray (boolean)
-        window : int
-          Minimum run length.
-        dim : str
-          Dimension along which to calculate consecutive run (default: 'time').
-        ufunc_1dim : Union[str, bool]
-          Use the 1d 'ufunc' version of this function : default (auto) will attempt to select optimal
-          usage based on number of data points.  Using 1D_ufunc=True is typically more efficient
-          for dataarray with a small number of gridpoints.
-        Returns
-        -------
-        xr.DataArray
-          Number of distinct runs of a minimum length (int).
-        """
+    Parameters
+    ----------
+    da: xr.DataArray
+      Input N-dimensional DataArray (boolean)
+    window : int
+      Minimum run length.
+    dim : str
+      Dimension along which to calculate consecutive run (default: 'time').
+    ufunc_1dim : Union[str, bool]
+      Use the 1d 'ufunc' version of this function : default (auto) will attempt to select optimal
+      usage based on number of data points.  Using 1D_ufunc=True is typically more efficient
+      for dataarray with a small number of gridpoints.
+    Returns
+    -------
+    xr.DataArray
+      Number of distinct runs of a minimum length (int).
+    """
     if ufunc_1dim == "auto":
         npts = get_npts(da)
         ufunc_1dim = npts <= npts_opt
@@ -159,24 +159,24 @@ def windowed_run_count(
 ) -> xr.DataArray:
     """Return the number of consecutive true values in array for runs at least as long as given duration.
 
-        Parameters
-        ----------
-        da: xr.DataArray
-          Input N-dimensional DataArray (boolean).
-        window : int
-          Minimum run length.
-        dim : str
-          Dimension along which to calculate consecutive run (default: 'time').
-        ufunc_1dim : Union[str, bool]
-          Use the 1d 'ufunc' version of this function : default (auto) will attempt to select optimal
-          usage based on number of data points. Using 1D_ufunc=True is typically more efficient
-          for dataarray with a small number of gridpoints.
-
-        Returns
-        -------
-        xr.DataArray
-          Total number of true values part of a consecutive runs of at least `window` long.
-        """
+    Parameters
+    ----------
+    da: xr.DataArray
+      Input N-dimensional DataArray (boolean).
+    window : int
+      Minimum run length.
+    dim : str
+      Dimension along which to calculate consecutive run (default: 'time').
+    ufunc_1dim : Union[str, bool]
+      Use the 1d 'ufunc' version of this function : default (auto) will attempt to select optimal
+      usage based on number of data points. Using 1D_ufunc=True is typically more efficient
+      for dataarray with a small number of gridpoints.
+
+    Returns
+    -------
+    xr.DataArray
+      Total number of true values part of a consecutive runs of at least `window` long.
+    """
     if ufunc_1dim == "auto":
         npts = get_npts(da)
         ufunc_1dim = npts <= npts_opt
@@ -276,8 +276,8 @@ def last_run(
       DateTimeAccessor object to use (ex: 'dayofyear').
     ufunc_1dim : Union[str, bool]
       Use the 1d 'ufunc' version of this function : default (auto) will attempt to select optimal
-      usage based on number of data points.  Using 1D_ufunc=True is typically more efficient
-      for dataarray with a small number of gridpoints.
+      usage based on number of data points.  Using `1D_ufunc=True` is typically more efficient
+      for a DataArray with a small number of grid points.
 
     Returns
     -------
@@ -293,6 +293,163 @@ def last_run(
     return out
 
 
+def run_length_with_date(
+    da: xr.DataArray, window: int, date: str = "07-01", dim: str = "time",
+):
+    """Return the length of the longest consecutive run of True values found
+    to be semi-continuous before and after a given date.
+
+    Parameters
+    ----------
+    da : xr.DataArray
+      Input N-dimensional DataArray (boolean)
+    window : int
+      Minimum duration of consecutive run to accumulate values.
+    date:
+      The date that a run must include to be considered valid.
+    dim : str
+      Dimension along which to calculate consecutive run (default: 'time').
+
+    Returns
+    -------
+    out : xr.DataArray
+      Length of longest run of True values along a given dimension inclusive of a given date.
+
+    Notes
+    -----
+    The run can include holes of False or NaN values, so long as they do not exceed the window size.
+    """
+    include_date = datetime.strptime(date, "%m-%d").timetuple().tm_yday
+
+    def _crld(rlda, rlwindow, rldoy, rldim):
+        mid_index = np.where(rlda.time.dt.dayofyear == rldoy)[0]
+        if (
+            mid_index.size == 0
+        ):  # The date is not within the group. Happens at boundaries.
+            all_nans = xr.full_like(rlda.isel(time=0), np.nan)
+            all_nans.attrs = {}
+            return all_nans
+        end = first_run(
+            (~rlda).where(rlda.time >= rlda.time[mid_index][0]),
+            window=rlwindow,
+            dim=rldim,
+        )
+        beg = first_run(rlda, window=rlwindow, dim=rldim)
+        sl = end - beg
+        sl = xr.where(
+            beg.isnull() & end.notnull(), 0, sl
+        )  # If series is never triggered
+        sl = xr.where(
+            beg.notnull() & end.isnull(), rlda.time.size - beg, sl
+        )  # If series is not ended by end of resample time frequency
+        return sl.where(sl >= 0)
+
+    return _crld(rlda=da, rlwindow=window, rldoy=include_date, rldim=dim)
+
+
+def run_end_after_date(
+    da: xr.DataArray,
+    window: int,
+    date: str = "07-01",
+    dim: str = "time",
+    coord: str = "dayofyear",
+):
+    """Return the index of the first item after the end of a run after a given date.
+
+    Parameters
+    ----------
+    da : xr.DataArray
+      Input N-dimensional DataArray (boolean)
+    window : int
+      Minimum duration of consecutive run to accumulate values.
+    date : str
+      The date after which to look for the end of a run.
+    dim : str
+      Dimension along which to calculate consecutive run (default: 'time').
+    coord : Optional[str]
+      If not False, the function returns values along `dim` instead of indexes.
+      If `dim` has a datetime dtype, `coord` can also be a str of the name of the
+      DateTimeAccessor object to use (ex: 'dayofyear').
+
+    Returns
+    -------
+    out : xr.DataArray
+      Index (or coordinate if `coord` is not False) of last item in last valid run. Returns np.nan if there are no valid run.
+    """
+    after_date = datetime.strptime(date, "%m-%d").timetuple().tm_yday
+
+    def _cred(rlda, rlwindow, rldoy, rldim, rlcoord):
+        mid_idx = np.where(rlda.time.dt.dayofyear == rldoy)[0]
+        if (
+            mid_idx.size == 0
+        ):  # The date is not within the group. Happens at boundaries.
+            all_nans = xr.full_like(rlda.isel(time=0), np.nan)
+            all_nans.attrs = {}
+            return all_nans
+        end = first_run(
+            (~rlda).where(rlda.time >= rlda.time[mid_idx][0]),
+            window=rlwindow,
+            dim=rldim,
+            coord=rlcoord,
+        )
+        beg = first_run(
+            rlda.where(rlda.time < rlda.time[mid_idx][0]), window=rlwindow, dim=rldim,
+        )
+        end = xr.where(
+            end.isnull() & beg.notnull(), rlda.time.isel(time=-1).dt.dayofyear, end
+        )
+        return end.where(beg.notnull())
+
+    return _cred(rlda=da, rlwindow=window, rldoy=after_date, rldim=dim, rlcoord=coord)
+
+
+def last_run_before_date(
+    da: xr.DataArray,
+    window: int,
+    date: str = "07-01",
+    dim: str = "time",
+    coord: str = "dayofyear",
+):
+    """Return the index of the last item of the last run before a given date.
+
+    Parameters
+    ----------
+    da : xr.DataArray
+      Input N-dimensional DataArray (boolean)
+    window : int
+      Minimum duration of consecutive run to accumulate values.
+    date : str
+      The date before which to look for the last event.
+    dim : str
+      Dimension along which to calculate consecutive run (default: 'time').
+    coord : Optional[str]
+      If not False, the function returns values along `dim` instead of indexes.
+      If `dim` has a datetime dtype, `coord` can also be a str of the name of the
+      DateTimeAccessor object to use (ex: 'dayofyear').
+
+    Returns
+    -------
+    out : xr.DataArray
+      Index (or coordinate if `coord` is not False) of last item in last valid run. Returns np.nan if there are no valid run.
+    """
+    before_date = datetime.strptime(date, "%m-%d").timetuple().tm_yday
+
+    def _clrbd(rlda, rlwindow, rlbefore_date, rldim, rlcoord):
+        mid_idx = np.where(rlda.time.dt.dayofyear == rlbefore_date)[0]
+        if (
+            mid_idx.size == 0
+        ):  # The date is not within the group. Happens at boundaries.
+            all_nans = xr.full_like(rlda.isel(time=0), np.nan)
+            all_nans.attrs = dict()
+            return all_nans
+        run = rlda.where(rlda.time <= rlda.time[mid_idx][0])
+        return last_run(run, window=rlwindow, dim=rldim, coord=rlcoord)
+
+    return _clrbd(
+        rlda=da, rlwindow=window, rlbefore_date=before_date, rldim=dim, rlcoord=coord
+    )
+
+
 def rle_1d(
     arr: Union[int, float, bool, Sequence[Union[int, float, bool]]]
 ) -> Tuple[np.array, np.array, np.array]:
@@ -317,7 +474,6 @@ def rle_1d(
     >>> a = [1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3]
     >>> rle_1d(a)
     (array([1, 2, 3]), array([2, 4, 6]), array([0, 2, 6]))
-
     """
     ia = np.asarray(arr)
     n = len(ia)
@@ -394,14 +550,13 @@ def windowed_run_count_1d(arr: Sequence[bool], window: int) -> int:
     return np.where(v * rl >= window, rl, 0).sum()
 
 
-def windowed_run_events_1d(arr, window: int):
+def windowed_run_events_1d(arr: Sequence[bool], window: int):
     """Return the number of runs of a minimum length.
 
     Parameters
     ----------
-    arr : bool array
-      Input array
-
+    arr : Sequence[bool]
+      Input array (bool)
     window : int
       Minimum run length.
 
@@ -414,14 +569,14 @@ def windowed_run_events_1d(arr, window: int):
     return (v * rl >= window).sum()
 
 
-def windowed_run_count_ufunc(x, window) -> xr.apply_ufunc:
+def windowed_run_count_ufunc(x: Sequence[bool], window: int) -> xr.apply_ufunc:
     """Dask-parallel version of windowed_run_count_1d, ie the number of consecutive true values in
     array for runs at least as long as given duration.
 
     Parameters
     ----------
-    x : bool array
-      Input array
+    x : Sequence[bool]
+      Input array (bool)
     window : int
       Minimum duration of consecutive run to accumulate values.
 
@@ -502,7 +657,7 @@ def first_run_ufunc(x: xr.DataArray, window: int, dim: str = "time",) -> xr.appl
     x : xr.DataArray
       Input array (bool)
     window : int
-    index: Optional[str]
+    dim: Optional[str]
 
     Returns
     -------
@@ -549,8 +704,5 @@ def lazy_indexing(da: xr.DataArray, index: xr.DataArray):
     index = index.fillna(0).astype(int)
     func = partial(_index_from_1d_array, da)
 
-    if isinstance(index, dsk.Array):
-        out = dsk.map_blocks(func, index, dtype=index.dtype)
-    else:
-        out = func(index)
+    out = index.map_blocks(func)
     return out.where(~invalid)
diff --git a/xclim/locales/__init__.py b/xclim/locales/__init__.py
index 6fb4fc45..29910f07 100644
--- a/xclim/locales/__init__.py
+++ b/xclim/locales/__init__.py
@@ -316,7 +316,7 @@ def generate_local_dict(locale: str, init_english: bool = False):
     if best_locale is not None:
         locname, attrs = get_local_dict(best_locale)
         for ind_name in attrs.copy().keys():
-            if ind_name not in indicators:
+            if ind_name != "attrs_mapping" and ind_name not in indicators:
                 attrs.pop(ind_name)
     else:
         attrs = {}
diff --git a/xclim/locales/fr.json b/xclim/locales/fr.json
index 0beeffc7..b9170273 100644
--- a/xclim/locales/fr.json
+++ b/xclim/locales/fr.json
@@ -1,26 +1,95 @@
 {
   "attrs_mapping": {
-    "modifiers": ["m", "f", "mpl", "fpl"],
-    "YS": ["annuel", "annuelle", "annuels", "annuelles"],
-    "MS": ["mensuel", "mensuelle", "mensuels", "mensuelles"],
-    "QS-DEC": ["saisonnier", "saisonnière", "saisonniers", "saisonnières"],
-    "DJF": ["hivernal", "hivernale", "hivernaux", "hivernales"],
-    "MAM": ["printanier", "printanière", "printaniers", "printanières"],
-    "JJA": ["estival", "estivale", "estivaux", "estivales"],
-    "SON": ["automnal", "automnale", "automnaux", "automnales"],
-    "norm": ["normal", "normale", "normaux", "normales"],
-    "m1": ["janvier"],
-    "m2": ["février"],
-    "m3": ["mars"],
-    "m4": ["avril"],
-    "m5": ["mai"],
-    "m6": ["juin"],
-    "m7": ["juillet"],
-    "m8": ["août"],
-    "m9": ["septembre"],
-    "m10": ["octobre"],
-    "m11": ["novembre"],
-    "m12": ["decembre"]
+    "modifiers": [
+      "m",
+      "f",
+      "mpl",
+      "fpl"
+    ],
+    "YS": [
+      "annuel",
+      "annuelle",
+      "annuels",
+      "annuelles"
+    ],
+    "MS": [
+      "mensuel",
+      "mensuelle",
+      "mensuels",
+      "mensuelles"
+    ],
+    "QS-DEC": [
+      "saisonnier",
+      "saisonnière",
+      "saisonniers",
+      "saisonnières"
+    ],
+    "DJF": [
+      "hivernal",
+      "hivernale",
+      "hivernaux",
+      "hivernales"
+    ],
+    "MAM": [
+      "printanier",
+      "printanière",
+      "printaniers",
+      "printanières"
+    ],
+    "JJA": [
+      "estival",
+      "estivale",
+      "estivaux",
+      "estivales"
+    ],
+    "SON": [
+      "automnal",
+      "automnale",
+      "automnaux",
+      "automnales"
+    ],
+    "norm": [
+      "normal",
+      "normale",
+      "normaux",
+      "normales"
+    ],
+    "m1": [
+      "janvier"
+    ],
+    "m2": [
+      "février"
+    ],
+    "m3": [
+      "mars"
+    ],
+    "m4": [
+      "avril"
+    ],
+    "m5": [
+      "mai"
+    ],
+    "m6": [
+      "juin"
+    ],
+    "m7": [
+      "juillet"
+    ],
+    "m8": [
+      "août"
+    ],
+    "m9": [
+      "septembre"
+    ],
+    "m10": [
+      "octobre"
+    ],
+    "m11": [
+      "novembre"
+    ],
+    "m12": [
+      "decembre"
+    ]
   },
   "atmos.rain_frzgr": {
     "long_name": "Nombre de jours de pluie sur sol gelé",
@@ -406,5 +475,40 @@
     "comment": "",
     "title": "Étendue de la banquise",
     "abstract": "Mesure de l'étendue de toutes les régions où la concentration de glace de mer excède un seuil."
+  },
+  "atmos.dry_days": {
+    "long_name": "Nombre de jours secs (precip < {thresh})",
+    "description": "Nombre {freq} de jours où la précipitation journalière est sous {thresh}.",
+    "comment": "",
+    "title": "Jours secs",
+    "abstract": "Nombre de jours où la précipitation journalière est sous un seuil."
+  },
+  "atmos.hot_spell_frequency": {
+    "long_name": "Nombre de périodes chaudes (Tmax > {thresh_tasmax} durant >= {window} jours)",
+    "description": "Nombre {freq} de périodes chaudes. Une période chaude se produit lorsque la température journalière maxmale excède {thresh_tasmax} durant au moins {window} jours.",
+    "comment": "",
+    "title": "Fréquences des périodes chaudes",
+    "abstract": "Nombre de périodes chaudes dans un période donnée. Une période chaude se produit lorsque la température journalière maxmale excède un seuil spécifique durant un minimum de jours."
+  },
+  "atmos.last_spring_frost": {
+    "long_name": "Jour de l'année du dernier gel printanier",
+    "description": "Jour de l'année du dernier gel printanier, defini comme le dernier jour où la température journalière moyenne reste sous {thresh} durant au moins {window} jours avant une certaine date.",
+    "comment": "",
+    "title": "Jour de l'année du dernier gel printanier",
+    "abstract": "Dernier jour où la température est inférieur à un seuil durant un certain nombre de jours, limité par une date finale."
+  },
+  "atmos.consecutive_frost_free_days": {
+    "long_name": "Nombre maximal de jours consécutifs avec Tmin > {thresh}",
+    "description": "Nombre maximal {freq} de jours consécutifs où la température journalière minimale est sous {thresh}.",
+    "comment": "",
+    "title": "Nombre maximal de jours de gel consécutifs (Tn > 0℃)",
+    "abstract": "La nombre maximal de jours consécutifs où la température journalière minimale est sous un certain seuil."
+  },
+  "atmos.growing_season_end": {
+    "long_name": "Jour de l'année de la fin de la saison de végétation",
+    "description": "Jour de l'année de la fin de la saison de végétation. La saison de végétation est définie comme l'intervalle entre la première série de {window} jours où la température journalière est au-dessus de {thresh} et la première série (après le {mid_date}) de {window} jours où elle est sous {thresh}.",
+    "comment": "",
+    "title": "Fin de la saison de végétation",
+    "abstract": "Jour de l'année de la fin de la saison de végétation. La saison de végétation est définie comme l'intervalle entre la première série de N jours où la température journalière est au-dessus d'un seuil et la première série (après une certaine date) de N jours où elle est sous ce même seuil"
   }
 }


Test patch:
diff --git a/tests/test_indices.py b/tests/test_indices.py
index f3a3ad36..f658ff25 100644
--- a/tests/test_indices.py
+++ b/tests/test_indices.py
@@ -13,7 +13,6 @@
 # For correctness, I think it would be useful to use a small dataset and run the original ICCLIM indicators on it,
 # saving the results in a reference netcdf dataset. We could then compare the hailstorm output to this reference as
 # a first line of defense.
-# import cftime
 import calendar
 import os
 
@@ -142,64 +141,55 @@ class TestColdSpellDays:
 
 
 class TestConsecutiveFrostDays:
-    @staticmethod
-    def time_series(values):
-        coords = pd.date_range(
-            "7/1/2000", periods=len(values), freq=pd.DateOffset(days=1)
-        )
-        return xr.DataArray(
-            values,
-            coords=[coords],
-            dims="time",
-            attrs={
-                "standard_name": "air_temperature",
-                "cell_methods": "time: minimum within days",
-                "units": "K",
-            },
-        )
-
-    def test_one_freeze_day(self):
-        a = self.time_series(np.array([3, 4, 5, -1, 3]) + K2C)
+    def test_one_freeze_day(self, tasmin_series):
+        a = tasmin_series(np.array([3, 4, 5, -1, 3]) + K2C)
         cfd = xci.consecutive_frost_days(a)
         assert cfd == 1
         assert cfd.time.dt.year == 2000
 
-    def test_no_freeze(self):
-        a = self.time_series(np.array([3, 4, 5, 1, 3]) + K2C)
+    def test_no_freeze(self, tasmin_series):
+        a = tasmin_series(np.array([3, 4, 5, 1, 3]) + K2C)
         cfd = xci.consecutive_frost_days(a)
         assert cfd == 0
 
-    def test_all_year_freeze(self):
-        a = self.time_series(np.zeros(365) - 10 + K2C)
+    def test_all_year_freeze(self, tasmin_series):
+        a = tasmin_series(np.zeros(365) - 10 + K2C)
         cfd = xci.consecutive_frost_days(a)
         assert cfd == 365
 
 
-class TestCoolingDegreeDays:
-    @staticmethod
-    def time_series(values):
-        coords = pd.date_range(
-            "7/1/2000", periods=len(values), freq=pd.DateOffset(days=1)
-        )
-        return xr.DataArray(
-            values,
-            coords=[coords],
-            dims="time",
-            attrs={
-                "standard_name": "air_temperature",
-                "cell_methods": "time: mean within days",
-                "units": "K",
-            },
-        )
+class TestMaximumConsecutiveFrostFreeDays:
+    def test_one_freeze_day(self, tasmin_series):
+        a = tasmin_series(np.array([3, 4, 5, -1, 3]) + K2C)
+        ffd = xci.maximum_consecutive_frost_free_days(a)
+        assert ffd == 3
+        assert ffd.time.dt.year == 2000
+
+    def test_two_freeze_days_with_threshold(self, tasmin_series):
+        a = tasmin_series(np.array([3, 4, 5, -0.8, -2, 3]) + K2C)
+        ffd = xci.maximum_consecutive_frost_free_days(a, thresh="-1 degC")
+        assert ffd == 4
+
+    def test_no_freeze(self, tasmin_series):
+        a = tasmin_series(np.array([3, 4, 5, 1, 3]) + K2C)
+        ffd = xci.maximum_consecutive_frost_free_days(a)
+        assert ffd == 5
+
+    def test_all_year_freeze(self, tasmin_series):
+        a = tasmin_series(np.zeros(365) - 10 + K2C)
+        ffd = xci.maximum_consecutive_frost_free_days(a)
+        assert np.all(ffd) == 0
+
 
-    def test_no_cdd(self):
-        a = self.time_series(np.array([10, 15, -5, 18]) + K2C)
+class TestCoolingDegreeDays:
+    def test_no_cdd(self, tas_series):
+        a = tas_series(np.array([10, 15, -5, 18]) + K2C)
         cdd = xci.cooling_degree_days(a)
         assert cdd == 0
         assert cdd.units == "C days"
 
-    def test_cdd(self):
-        a = self.time_series(np.array([20, 25, -15, 19]) + K2C)
+    def test_cdd(self, tas_series):
+        a = tas_series(np.array([20, 25, -15, 19]) + K2C)
         cdd = xci.cooling_degree_days(a)
         assert cdd == 10
 
@@ -259,6 +249,16 @@ class TestDailyPrIntensity:
         np.testing.assert_array_almost_equal(out[0], 2.5)
 
 
+class TestLastSpringFrost:
+    def test_simple(self, tas_series):
+        a = np.zeros(365)
+        a[180:270] = 303.15
+        tas = tas_series(a, start="2000/1/1")
+
+        lsf = xci.last_spring_frost(tas)
+        assert lsf == 180
+
+
 class TestDaysOverPrecipThresh:
     def test_simple(self, pr_series, per_doy):
         a = np.zeros(365)
@@ -334,6 +334,28 @@ class TestGrowingDegreeDays:
         assert xci.growing_degree_days(da)[0] == 1
 
 
+class TestGrowingSeasonEnd:
+    @pytest.mark.parametrize(
+        "d1,d2,mid_date,expected",
+        [
+            ("1950-01-01", "1951-01-01", "07-01", np.nan),  # No growing season
+            ("2000-01-01", "2000-12-31", "07-01", 365),  # All year growing season
+            ("2000-07-10", "2001-01-01", "07-01", np.nan),  # End happens before start
+            ("2000-06-15", "2000-07-15", "07-01", 198),  # Normal case
+            ("2000-06-15", "2000-07-25", "07-15", 208),  # PCC Case
+            ("2000-06-15", "2000-07-15", "10-01", 274),  # Late mid_date
+            ("2000-06-15", "2000-07-15", "01-10", np.nan),  # Early mid_date
+        ],
+    )
+    def test_varying_mid_dates(self, tas_series, d1, d2, mid_date, expected):
+        # generate a year of data
+        tas = tas_series(np.zeros(365), start="2000/1/1")
+        warm_period = tas.sel(time=slice(d1, d2))
+        tas = tas.where(~tas.time.isin(warm_period.time), 280)
+        gs_end = xci.growing_season_end(tas, mid_date=mid_date)
+        np.testing.assert_array_equal(gs_end, expected)
+
+
 class TestGrowingSeasonLength:
     @pytest.mark.parametrize(
         "d1,d2,expected",
@@ -386,87 +408,134 @@ class TestHeatWaveIndex:
 
 
 class TestHeatWaveFrequency:
-    def test_1d(self, tasmax_series, tasmin_series):
+    @pytest.mark.parametrize(
+        "thresh_tasmin,thresh_tasmax,window,expected",
+        [
+            ("22 C", "30 C", 3, 2),  # Some HW
+            ("22 C", "30 C", 4, 1),  # No HW
+            ("10 C", "10 C", 3, 1),  # One long HW
+            ("40 C", "40 C", 3, 0),  # Windowed
+        ],
+    )
+    def test_1d(
+        self,
+        tasmax_series,
+        tasmin_series,
+        thresh_tasmin,
+        thresh_tasmax,
+        window,
+        expected,
+    ):
         tn = tasmin_series(np.asarray([20, 23, 23, 23, 23, 22, 23, 23, 23, 23]) + K2C)
         tx = tasmax_series(np.asarray([29, 31, 31, 31, 29, 31, 31, 31, 31, 31]) + K2C)
 
-        # some hw
-        hwf = xci.heat_wave_frequency(
-            tn, tx, thresh_tasmin="22 C", thresh_tasmax="30 C"
-        )
-        np.testing.assert_allclose(hwf.values, 2)
-        hwf = xci.heat_wave_frequency(
-            tn, tx, thresh_tasmin="22 C", thresh_tasmax="30 C", window=4
-        )
-        np.testing.assert_allclose(hwf.values, 1)
-        # one long hw
-        hwf = xci.heat_wave_frequency(
-            tn, tx, thresh_tasmin="10 C", thresh_tasmax="10 C"
-        )
-        np.testing.assert_allclose(hwf.values, 1)
-        # no hw
         hwf = xci.heat_wave_frequency(
-            tn, tx, thresh_tasmin="40 C", thresh_tasmax="40 C"
+            tn,
+            tx,
+            thresh_tasmin=thresh_tasmin,
+            thresh_tasmax=thresh_tasmax,
+            window=window,
         )
-        np.testing.assert_allclose(hwf.values, 0)
+        np.testing.assert_allclose(hwf.values, expected)
 
 
 class TestHeatWaveMaxLength:
-    def test_1d(self, tasmax_series, tasmin_series):
+    @pytest.mark.parametrize(
+        "thresh_tasmin,thresh_tasmax,window,expected",
+        [
+            ("22 C", "30 C", 3, 4),  # Some HW
+            ("10 C", "10 C", 3, 10),  # One long HW
+            ("40 C", "40 C", 3, 0),  # No HW
+            ("22 C", "30 C", 5, 0),  # Windowed
+        ],
+    )
+    def test_1d(
+        self,
+        tasmax_series,
+        tasmin_series,
+        thresh_tasmin,
+        thresh_tasmax,
+        window,
+        expected,
+    ):
         tn = tasmin_series(np.asarray([20, 23, 23, 23, 23, 22, 23, 23, 23, 23]) + K2C)
         tx = tasmax_series(np.asarray([29, 31, 31, 31, 29, 31, 31, 31, 31, 31]) + K2C)
 
-        # some hw
         hwml = xci.heat_wave_max_length(
-            tn, tx, thresh_tasmin="22 C", thresh_tasmax="30 C"
+            tn,
+            tx,
+            thresh_tasmin=thresh_tasmin,
+            thresh_tasmax=thresh_tasmax,
+            window=window,
         )
-        np.testing.assert_allclose(hwml.values, 4)
-
-        # one long hw
-        hwml = xci.heat_wave_max_length(
-            tn, tx, thresh_tasmin="10 C", thresh_tasmax="10 C"
-        )
-        np.testing.assert_allclose(hwml.values, 10)
-
-        # no hw
-        hwml = xci.heat_wave_max_length(
-            tn, tx, thresh_tasmin="40 C", thresh_tasmax="40 C"
-        )
-        np.testing.assert_allclose(hwml.values, 0)
-
-        hwml = xci.heat_wave_max_length(
-            tn, tx, thresh_tasmin="22 C", thresh_tasmax="30 C", window=5
-        )
-        np.testing.assert_allclose(hwml.values, 0)
+        np.testing.assert_allclose(hwml.values, expected)
 
 
 class TestHeatWaveTotalLength:
-    def test_1d(self, tasmax_series, tasmin_series):
+    @pytest.mark.parametrize(
+        "thresh_tasmin,thresh_tasmax,window,expected",
+        [
+            ("22 C", "30 C", 3, 7),  # Some HW
+            ("10 C", "10 C", 3, 10),  # One long HW
+            ("40 C", "40 C", 3, 0),  # No HW
+            ("22 C", "30 C", 5, 0),  # Windowed
+        ],
+    )
+    def test_1d(
+        self,
+        tasmax_series,
+        tasmin_series,
+        thresh_tasmin,
+        thresh_tasmax,
+        window,
+        expected,
+    ):
         tn = tasmin_series(np.asarray([20, 23, 23, 23, 23, 22, 23, 23, 23, 23]) + K2C)
         tx = tasmax_series(np.asarray([29, 31, 31, 31, 29, 31, 31, 31, 31, 31]) + K2C)
 
         # some hw
         hwml = xci.heat_wave_total_length(
-            tn, tx, thresh_tasmin="22 C", thresh_tasmax="30 C"
+            tn,
+            tx,
+            thresh_tasmin=thresh_tasmin,
+            thresh_tasmax=thresh_tasmax,
+            window=window,
         )
-        np.testing.assert_allclose(hwml.values, 7)
+        np.testing.assert_allclose(hwml.values, expected)
 
-        # one long hw
-        hwml = xci.heat_wave_total_length(
-            tn, tx, thresh_tasmin="10 C", thresh_tasmax="10 C"
-        )
-        np.testing.assert_allclose(hwml.values, 10)
 
-        # no hw
-        hwml = xci.heat_wave_total_length(
-            tn, tx, thresh_tasmin="40 C", thresh_tasmax="40 C"
-        )
-        np.testing.assert_allclose(hwml.values, 0)
+class TestHotSpellFrequency:
+    @pytest.mark.parametrize(
+        "thresh_tasmax,window,expected",
+        [
+            ("30 C", 3, 2),  # Some HS
+            ("30 C", 4, 1),  # One long HS
+            ("10 C", 3, 1),  # No HS
+            ("40 C", 5, 0),  # Windowed
+        ],
+    )
+    def test_1d(self, tasmax_series, thresh_tasmax, window, expected):
+        tx = tasmax_series(np.asarray([29, 31, 31, 31, 29, 31, 31, 31, 31, 31]) + K2C)
 
-        hwml = xci.heat_wave_total_length(
-            tn, tx, thresh_tasmin="22 C", thresh_tasmax="30 C", window=5
-        )
-        np.testing.assert_allclose(hwml.values, 0)
+        hsf = xci.hot_spell_frequency(tx, thresh_tasmax=thresh_tasmax, window=window)
+        np.testing.assert_allclose(hsf.values, expected)
+
+
+class TestHotSpellMaxLength:
+    @pytest.mark.parametrize(
+        "thresh_tasmax,window,expected",
+        [
+            ("30 C", 3, 5),  # Some HS
+            ("10 C", 3, 10),  # One long HS
+            ("40 C", 3, 0),  # No HS
+            ("30 C", 5, 5),  # Windowed
+        ],
+    )
+    def test_1d(self, tasmax_series, thresh_tasmax, window, expected):
+        tx = tasmax_series(np.asarray([29, 31, 31, 31, 29, 31, 31, 31, 31, 31]) + K2C)
+
+        hsml = xci.hot_spell_max_length(tx, thresh_tasmax=thresh_tasmax, window=window)
+        np.testing.assert_allclose(hsml.values, expected)
 
 
 class TestTnDaysBelow:
@@ -752,47 +821,20 @@ class TestTas:
 
 
 class TestTxMin:
-    @staticmethod
-    def time_series(values):
-        coords = pd.date_range(
-            "7/1/2000", periods=len(values), freq=pd.DateOffset(days=1)
-        )
-        return xr.DataArray(
-            values,
-            coords=[coords],
-            dims="time",
-            attrs={
-                "standard_name": "air_temperature",
-                "cell_methods": "time: maximum within days",
-                "units": "K",
-            },
-        )
+    def test_simple(self, tasmax_series):
+        a = tasmax_series(np.array([20, 25, -15, 19]))
+        txm = xci.tx_min(a, freq="YS")
+        assert txm == -15
 
 
 class TestTxMean:
-    @staticmethod
-    def time_series(values):
-        coords = pd.date_range(
-            "7/1/2000", periods=len(values), freq=pd.DateOffset(days=1)
-        )
-        return xr.DataArray(
-            values,
-            coords=[coords],
-            dims="time",
-            attrs={
-                "standard_name": "air_temperature",
-                "cell_methods": "time: maximum within days",
-                "units": "K",
-            },
-        )
-
-    def test_attrs(self):
-        a = self.time_series(np.array([320, 321, 322, 323, 324]))
+    def test_attrs(self, tasmax_series):
+        a = tasmax_series(np.array([320, 321, 322, 323, 324]))
         txm = xci.tx_mean(a, freq="YS")
         assert txm == 322
         assert txm.units == "K"
 
-        a = self.time_series(np.array([20, 21, 22, 23, 24]))
+        a = tasmax_series(np.array([20, 21, 22, 23, 24]))
         a.attrs["units"] = "C"
         txm = xci.tx_mean(a, freq="YS")
 
@@ -801,24 +843,8 @@ class TestTxMean:
 
 
 class TestTxMax:
-    @staticmethod
-    def time_series(values):
-        coords = pd.date_range(
-            "7/1/2000", periods=len(values), freq=pd.DateOffset(days=1)
-        )
-        return xr.DataArray(
-            values,
-            coords=[coords],
-            dims="time",
-            attrs={
-                "standard_name": "air_temperature",
-                "cell_methods": "time: maximum within days",
-                "units": "K",
-            },
-        )
-
-    def test_simple(self):
-        a = self.time_series(np.array([20, 25, -15, 19]))
+    def test_simple(self, tasmax_series):
+        a = tasmax_series(np.array([20, 25, -15, 19]))
         txm = xci.tx_max(a, freq="YS")
         assert txm == 25
 
@@ -986,12 +1012,38 @@ class TestWinterRainRatio:
 
 # I'd like to parametrize some of these tests so we don't have to write individual tests for each indicator.
 class TestTG:
-    def test_cmip3(self, cmip3_day_tas):
+    @staticmethod
+    @pytest.fixture(scope="session")
+    def cmip3_day_tas():
+        # xr.set_options(enable_cftimeindex=False)
+        ds = xr.open_dataset(
+            os.path.join(
+                TESTS_DATA, "cmip3", "tas.sresb1.giss_model_e_r.run1.atm.da.nc"
+            )
+        )
+        yield ds.tas
+        ds.close()
+
+    def test_cmip3_tgmean(self, cmip3_day_tas):
         pytest.importorskip("xarray", "0.11.4")
         xci.tg_mean(cmip3_day_tas)
 
-    def compare_against_icclim(self, cmip3_day_tas):
-        pass
+    def test_cmip3_tgmin(self, cmip3_day_tas):
+        pytest.importorskip("xarray", "0.11.4")
+        xci.tg_min(cmip3_day_tas)
+
+    def test_cmip3_tgmax(self, cmip3_day_tas):
+        pytest.importorskip("xarray", "0.11.4")
+        xci.tg_max(cmip3_day_tas)
+
+    def test_indice_against_icclim(self, cmip3_day_tas):
+        pytest.importorskip("xarray", "0.11.4")
+        from xclim import icclim
+
+        ind = xci.tg_mean(cmip3_day_tas)
+        icclim = icclim.TG(cmip3_day_tas)
+
+        np.testing.assert_array_equal(icclim, ind)
 
 
 @pytest.mark.skipif(
@@ -1040,16 +1092,6 @@ class TestFireWeatherIndex:
         )
 
 
-@pytest.fixture(scope="session")
-def cmip3_day_tas():
-    # xr.set_options(enable_cftimeindex=False)
-    ds = xr.open_dataset(
-        os.path.join(TESTS_DATA, "cmip3", "tas.sresb1.giss_model_e_r.run1.atm.da.nc")
-    )
-    yield ds.tas
-    ds.close()
-
-
 class TestWindConversion:
     da_uas = xr.DataArray(
         np.array([[3.6, -3.6], [-1, 0]]),
diff --git a/tests/test_run_length.py b/tests/test_run_length.py
index c1b16e59..072f5052 100644
--- a/tests/test_run_length.py
+++ b/tests/test_run_length.py
@@ -174,25 +174,6 @@ class TestFirstRun:
         out = rl.first_run(runs, window=1, dim="time", coord=coord, ufunc_1dim=use_1dim)
         np.testing.assert_array_equal(out.load(), expected)
 
-    @pytest.mark.parametrize(
-        "coord,expected",
-        [(False, 39), (True, np.datetime64("2000-02-09")), ("dayofyear", 40)],
-    )
-    @pytest.mark.parametrize(
-        "use_dask,use_1dim", [(True, False), (False, False), (False, True)]
-    )
-    def test_last_run(self, tas_series, coord, expected, use_dask, use_1dim):
-        t = np.zeros(60)
-        t[30:40] = 2
-        tas = tas_series(t, start="2000-01-01")
-        runs = xr.concat((tas, tas), dim="dim0")
-
-        if use_dask:
-            runs = runs.chunk({"time": 10, "dim0": 1})
-
-        out = rl.last_run(runs, window=1, dim="time", coord=coord, ufunc_1dim=use_1dim)
-        np.testing.assert_array_equal(out.load(), expected)
-
 
 class TestWindowedRunEvents:
     nc_pr = os.path.join(TESTS_DATA, "NRCANdaily", "nrcan_canada_daily_pr_1990.nc")
@@ -228,3 +209,99 @@ class TestWindowedRunCount:
             rl.windowed_run_count, window=4, dim="time", ufunc_1dim=False
         )
         np.testing.assert_array_equal(lt_orig, lt_Ndim)
+
+
+class TestLastRun:
+    @pytest.mark.parametrize(
+        "coord,expected",
+        [(False, 39), (True, np.datetime64("2000-02-09")), ("dayofyear", 40)],
+    )
+    @pytest.mark.parametrize(
+        "use_dask,use_1dim", [(True, False), (False, False), (False, True)]
+    )
+    def test_simple(self, tas_series, coord, expected, use_dask, use_1dim):
+        t = np.zeros(60)
+        t[30:40] = 2
+        tas = tas_series(t, start="2000-01-01")
+        runs = xr.concat((tas, tas), dim="dim0")
+
+        if use_dask:
+            runs = runs.chunk({"time": 10, "dim0": 1})
+
+        out = rl.last_run(runs, window=1, dim="time", coord=coord, ufunc_1dim=use_1dim)
+        np.testing.assert_array_equal(out.load(), expected)
+
+
+class TestRunsWithDates:
+    @pytest.mark.parametrize(
+        "date,end,expected",
+        [
+            ("07-01", 210, 70),
+            ("07-01", 190, 50),
+            ("04-01", 150, np.NaN),  # date falls early
+            ("11-01", 150, 164),  # date ends late
+        ],
+    )
+    @pytest.mark.parametrize("use_dask", [True, False])
+    def test_run_length_with_date(self, tas_series, date, end, expected, use_dask):
+        t = np.zeros(360)
+        t[140:end] = 1
+        tas = tas_series(t, start="2000-01-01")
+        runs = xr.concat((tas, tas), dim="dim0")
+        runs = runs == 1
+
+        if use_dask:
+            runs = runs.chunk({"time": 10, "dim0": 1})
+
+        out = rl.run_length_with_date(runs, window=1, dim="time", date=date,)
+        np.testing.assert_array_equal(np.mean(out.load()), expected)
+
+    @pytest.mark.parametrize(
+        "coord,date,end,expected",
+        [
+            ("dayofyear", "07-01", 210, 211),
+            (False, "07-01", 190, 190),
+            ("dayofyear", "04-01", 150, np.NaN),  # date falls early
+            ("dayofyear", "11-01", 150, 305),  # date ends late
+        ],
+    )
+    @pytest.mark.parametrize("use_dask", [True, False])
+    def test_run_end_after_date(self, tas_series, coord, date, end, expected, use_dask):
+        t = np.zeros(360)
+        t[140:end] = 1
+        tas = tas_series(t, start="2000-01-01")
+        runs = xr.concat((tas, tas), dim="dim0")
+        runs = runs == 1
+
+        if use_dask:
+            runs = runs.chunk({"time": 10, "dim0": 1})
+
+        out = rl.run_end_after_date(runs, window=1, date=date, dim="time", coord=coord)
+        np.testing.assert_array_equal(np.mean(out.load()), expected)
+
+    @pytest.mark.parametrize(
+        "coord,date,end,expected",
+        [
+            ("dayofyear", "07-01", 210, 182),
+            (False, "07-01", 190, 181),
+            ("dayofyear", "04-01", 150, np.NaN),  # date falls early
+            ("dayofyear", "11-01", 150, 150),  # date ends late
+        ],
+    )
+    @pytest.mark.parametrize("use_dask", [True, False])
+    def test_last_run_before_date(
+        self, tas_series, coord, date, end, expected, use_dask
+    ):
+        t = np.zeros(360)
+        t[140:end] = 1
+        tas = tas_series(t, start="2000-01-01")
+        runs = xr.concat((tas, tas), dim="dim0")
+        runs = runs == 1
+
+        if use_dask:
+            runs = runs.chunk({"time": 10, "dim0": 1})
+
+        out = rl.last_run_before_date(
+            runs, window=1, date=date, dim="time", coord=coord
+        )
+        np.testing.assert_array_equal(np.mean(out.load()), expected)