From 454cde5a714b5e0675c6fe8868af4c645bbde44c Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Fri, 8 Sep 2023 14:06:07 -0400
Subject: [PATCH 01/11] first pass at wind_power_potential

---
 tests/test_indices.py        |  16 ++++
 xclim/indices/_conversion.py | 153 ++++++++++++++++++++++++++++++++++-
 2 files changed, 168 insertions(+), 1 deletion(-)

diff --git a/tests/test_indices.py b/tests/test_indices.py
index 3f7176a38..418068c00 100644
--- a/tests/test_indices.py
+++ b/tests/test_indices.py
@@ -3626,3 +3626,19 @@ def test_late_frost_days_lat(self, tasmin_series):
         tasmin = tasmin_series(np.array([-1, 1, 2, -4, 0]) + K2C, start="30/3/2023")
         lfd = xci.frost_days(tasmin, date_bounds=("04-01", "06-30"))
         np.testing.assert_allclose(lfd, 1)
+
+
+class TestWindProfile:
+    def test_simple(self, sfcWind_series):
+        a = np.linspace(0, 100)
+        v = xci.wind_profile(sfcWind_series(a), h="100 m", h_r="10 m")
+        np.testing.assert_allclose(v, a * 10 ** (1 / 7))
+
+
+class TestWindPowerPotential:
+    def test_simple(self, sfcWind_series):
+        a = [2, 6, 20, 30]
+        p = xci.wind_power_potential(
+            sfcWind_series(a, units="m/s"), cut_in="4 m/s", rated="8 m/s"
+        )
+        np.testing.assert_allclose(p, [0, 0.5**3, 1, 0])
diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py
index 2ae948ffe..564aab2fa 100644
--- a/xclim/indices/_conversion.py
+++ b/xclim/indices/_conversion.py
@@ -5,7 +5,7 @@
 
 import numpy as np
 import xarray as xr
-from numba import float32, float64, vectorize  # noqa
+from numba import float32, float64, njit, vectorize  # noqa
 
 from xclim.core.calendar import date_range, datetime_to_decimal_year
 from xclim.core.units import (
@@ -53,6 +53,8 @@
     "uas_vas_2_sfcwind",
     "universal_thermal_climate_index",
     "wind_chill_index",
+    "wind_power_potential",
+    "wind_profile",
 ]
 
 
@@ -2017,3 +2019,152 @@ def mean_radiant_temperature(
         0.25,
     )
     return mrt.assign_attrs({"units": "K"})
+
+
+@declare_units(wind_speed="[speed]", h="[length]", h_r="[length]")
+def wind_profile(
+    wind_speed: xr.DataArray,
+    h: Quantified,
+    h_r: Quantified,
+    method: str = "power_law",
+    **kwds,
+):
+    r"""Wind speed at a given height estimated from the wind speed at a reference height.
+
+    Estimate the wind speed based on a power law profile relating wind speed to height above the surface:
+
+    .. math::
+
+        v = v_h \\left( \frac{h}{h_r} \right)^{\alpha}
+
+    Parameters
+    ----------
+    wind_speed : xarray.DataArray
+        Wind speed at the reference height.
+    h : Quantified
+        Height at which to compute the wind speed.
+    h_r : Quantified
+        Reference height.
+    method : {"power_law"}
+        Method to use. Currently only "power_law" is implemented.
+    kwds : dict
+        Additional keyword arguments to pass to the method. For `power_law`, this is `alpha`, which takes a default
+        value of 1/7.
+    """
+    # Convert units to meters
+    h = convert_units_to(h, "m")
+    h_r = convert_units_to(h_r, "m")
+
+    if method == "power_law":
+        alpha = kwds.pop("alpha", 1 / 7)
+        out = wind_speed * (h / h_r) ** alpha
+        out.attrs["units"] = wind_speed.attrs["units"]
+        return out
+    else:
+        raise NotImplementedError(f"Method {method} not implemented.")
+
+
+@declare_units(
+    wind_speed="[speed]",
+    air_density="[air_density]",
+    cut_in="[speed]",
+    rated="[speed]",
+    cut_out="[speed]",
+)
+def wind_power_potential(
+    wind_speed: xr.DataArray,
+    air_density: xr.DataArray = None,
+    cut_in: Quantified = "3.5 m/s",
+    rated: Quantified = "13 m/s",
+    cut_out: Quantified = "25 m/s",
+) -> xr.DataArray:
+    r"""Wind power potential estimated from a semi-idealized wind power production factor.
+
+    The actual power production of a wind farm can be estimated by multiplying its nominal (nameplate) capacity by the
+    wind power potential, which depends on wind speed at the hub height, the turbine specifications and air density.
+
+    Parameters
+    ----------
+    wind_speed : xarray.DataArray
+        Wind speed at the hub height. Use the `wind_profile` function to estimate from the surface wind speed. Note
+        that the temporal resolution of wind time series has a significant influence on the results. Mean daily wind
+        speeds yield lower values than hourly wind speeds.
+    air_density: xarray.DataArray
+        Air density at the hub height. Defaults to 1.225 kg/m³. This is worth changing if applying in cold or
+        mountainous regions with non-standard air density.
+    cut_in : Quantified
+        Cut-in wind speed. Default is 3.5 m/s.
+    rated : Quantified
+        Rated wind speed. Default is 13 m/s.
+    cut_out : Quantified
+        Cut-out wind speed. Default is 25 m/s.
+
+    Returns
+    -------
+    xr.DataArray
+      The power production factor. Multiply by the nominal capacity to get the actual power production.
+
+    See Also
+    --------
+    wind_profile : Estimate wind speed at the hub height from the surface wind speed.
+
+    Notes
+    -----
+    This estimate of wind power production is based on a semi-idealized power curve with four wind regimes specified
+    by the cut-in wind speed (:math:`u_i`), the rated speed (:math:`u_r`) and the cut-out speed (:math:`u_o`).
+    Power production is zero for wind speeds below the cut-in speed, increases cubically between the cut-in
+    and rated speed, is constant between the rated and cut-out speed, and is zero for wind speeds above the cut-out
+    speed to avoid damage to the turbine (10.1088/1748-9326/aab211):
+
+    .. math::
+
+        \begin{cases}
+        0,  &  v < u_i \\
+        (v - u_i)^3 / (u_r - u_i)^3,  & u_i <= v < u_r \\
+        1, & u_r <= v < u_o \\
+        0, v >= u_o
+        \\end{cases}
+
+    For non-standard air density, the wind speed is scaled using
+    :math:`v_n = v \\left( \frac{\rho}{\rho_0} \right)^{1/3}`.
+
+    Note that the temporal resolution of the wind speed time series has a significant influence on the results,
+    but percent changes in the wind power potential projections are similar across resolutions.
+    https://doi.org/10.1016/j.renene.2020.02.090
+    """
+    # Convert units
+    cut_in = convert_units_to(cut_in, wind_speed)
+    rated = convert_units_to(rated, wind_speed)
+    cut_out = convert_units_to(cut_out, wind_speed)
+
+    # Correct wind speed for air density
+    if air_density is not None:
+        default_air_density = convert_units_to("1.225 kg/m^3", air_density)
+        f = (air_density / default_air_density) ** (1 / 3)
+    else:
+        f = 1
+
+    v = wind_speed * f
+
+    if 1:
+        out = xr.zeros_like(v)
+        out = out + xr.where(
+            (v >= cut_in) * (v < rated), ((v - cut_in) / (rated - cut_in)) ** 3, 0
+        )
+        out = out + xr.where((v >= rated) * (v < cut_out), 1, 0)
+    else:
+        out = xr.apply_ufunc(_wind_power_factor, wind_speed, cut_in, rated, cut_out)
+
+    out.attrs["units"] = ""
+    return out
+
+
+@njit
+def _wind_power_factor(v, cut_in, rated, cut_out):
+    if v < cut_in:
+        return 0
+    if v < rated:
+        return ((v - cut_in) / (rated - cut_in)) ** 3
+    if v < cut_out:
+        return 1
+    return 0

From 0a2d40415462be733c877c312fdfa6f44ddefcd2 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Fri, 8 Sep 2023 14:42:02 -0400
Subject: [PATCH 02/11] Replaced wind power production logic by numba.vectorize
 function

---
 tests/test_indices.py        |  5 +++++
 xclim/indices/_conversion.py | 19 ++++++-------------
 2 files changed, 11 insertions(+), 13 deletions(-)

diff --git a/tests/test_indices.py b/tests/test_indices.py
index 418068c00..abd606001 100644
--- a/tests/test_indices.py
+++ b/tests/test_indices.py
@@ -3642,3 +3642,8 @@ def test_simple(self, sfcWind_series):
             sfcWind_series(a, units="m/s"), cut_in="4 m/s", rated="8 m/s"
         )
         np.testing.assert_allclose(p, [0, 0.5**3, 1, 0])
+
+    # def test_benchmark(self, sfcWind_series, benchmark):
+    #     a = np.random.rand(50000) * 30
+    #     v = sfcWind_series(a, units="m/s")
+    #     benchmark(xci.wind_power_potential, v)
diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py
index 564aab2fa..d5ca712c5 100644
--- a/xclim/indices/_conversion.py
+++ b/xclim/indices/_conversion.py
@@ -2146,25 +2146,18 @@ def wind_power_potential(
 
     v = wind_speed * f
 
-    if 1:
-        out = xr.zeros_like(v)
-        out = out + xr.where(
-            (v >= cut_in) * (v < rated), ((v - cut_in) / (rated - cut_in)) ** 3, 0
-        )
-        out = out + xr.where((v >= rated) * (v < cut_out), 1, 0)
-    else:
-        out = xr.apply_ufunc(_wind_power_factor, wind_speed, cut_in, rated, cut_out)
-
+    out = xr.apply_ufunc(_wind_power_factor, v, cut_in, rated, cut_out)
     out.attrs["units"] = ""
     return out
 
 
-@njit
+@vectorize
 def _wind_power_factor(v, cut_in, rated, cut_out):
+    """Wind power factor function"""
     if v < cut_in:
-        return 0
+        return 0.0
     if v < rated:
         return ((v - cut_in) / (rated - cut_in)) ** 3
     if v < cut_out:
-        return 1
-    return 0
+        return 1.0
+    return 0.0

From 36a542288c451eaf7194ba46a7340b9660725c97 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Fri, 8 Sep 2023 14:47:20 -0400
Subject: [PATCH 03/11] Update CHANGES

---
 CHANGES.rst | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/CHANGES.rst b/CHANGES.rst
index c059d400c..c833b2b0a 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -2,6 +2,16 @@
 Changelog
 =========
 
+v0.46.0 (unreleased)
+--------------------
+Contributors to this version: David Huard (:user:`huard`),
+
+New indicators
+^^^^^^^^^^^^^^
+* Add ``wind_profile`` to estimate the wind speed at different heights based on wind speed at a reference height.
+* Add ``wind_power_potential`` to estimate the potential for wind power production given wind speed at the turbine hub height and turbine specifications. (:issue:`1458`)
+
+
 v0.45.0 (2023-09-05)
 --------------------
 Contributors to this version: David Huard (:user:`huard`), Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user:`juliettelavoie`), Gabriel Rondeau-Genesse (:user:`RondeauG`), Marco Braun (:user:`vindelico`), Éric Dupuis (:user:`coxipi`).

From 2267c24457b39ac9053f209b7a891953c2f12795 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Fri, 8 Sep 2023 15:51:16 -0400
Subject: [PATCH 04/11] update bibtex. Debug math equations

---
 docs/references.bib          | 27 ++++++++++++++++++++++++
 xclim/indices/_conversion.py | 40 +++++++++++++++++++++++-------------
 2 files changed, 53 insertions(+), 14 deletions(-)

diff --git a/docs/references.bib b/docs/references.bib
index e269a50ff..1ccf861a5 100644
--- a/docs/references.bib
+++ b/docs/references.bib
@@ -2046,3 +2046,30 @@ @misc{usda_2012
 	author = {{USDA Agricultural Research Service}},
 	year = {2012},
 }
+
+@article{chen_2020,
+title = {Impacts of climate change on wind resources over North America based on NA-CORDEX},
+journal = {Renewable Energy},
+volume = {153},
+pages = {1428-1438},
+year = {2020},
+issn = {0960-1481},
+doi = {https://doi.org/10.1016/j.renene.2020.02.090},
+url = {https://www.sciencedirect.com/science/article/pii/S0960148120302858},
+author = {Liang Chen},
+keywords = {Climate change, Wind energy, NA-CORDEX, North America}
+}
+
+@article{tobin_2018,
+doi = {10.1088/1748-9326/aab211},
+url = {https://dx.doi.org/10.1088/1748-9326/aab211},
+year = {2018},
+month = {apr},
+publisher = {IOP Publishing},
+volume = {13},
+number = {4},
+pages = {044024},
+author = {I Tobin and W Greuell and S Jerez and F Ludwig and R Vautard and M T H van Vliet and F-M Bréon},
+title = {Vulnerabilities and resilience of European power generation to 1.5°C, 2°C and 3°C warming},
+journal = {Environmental Research Letters}
+}
diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py
index d5ca712c5..bcb76123c 100644
--- a/xclim/indices/_conversion.py
+++ b/xclim/indices/_conversion.py
@@ -2031,11 +2031,7 @@ def wind_profile(
 ):
     r"""Wind speed at a given height estimated from the wind speed at a reference height.
 
-    Estimate the wind speed based on a power law profile relating wind speed to height above the surface:
-
-    .. math::
-
-        v = v_h \\left( \frac{h}{h_r} \right)^{\alpha}
+    Estimate the wind speed based on a power law profile relating wind speed to height above the surface.
 
     Parameters
     ----------
@@ -2050,6 +2046,18 @@ def wind_profile(
     kwds : dict
         Additional keyword arguments to pass to the method. For `power_law`, this is `alpha`, which takes a default
         value of 1/7.
+
+    Notes
+    -----
+    The power law profile is given by
+
+    .. math::
+
+        v = v_r \left( \frac{h}{h_r} \right)^{\alpha},
+
+    where :math:`v_r` is the wind speed at the reference height, :math:`h` is the height at which the wind speed is
+    desired, and :math:`h_r` is the reference height.
+
     """
     # Convert units to meters
     h = convert_units_to(h, "m")
@@ -2114,23 +2122,27 @@ def wind_power_potential(
     by the cut-in wind speed (:math:`u_i`), the rated speed (:math:`u_r`) and the cut-out speed (:math:`u_o`).
     Power production is zero for wind speeds below the cut-in speed, increases cubically between the cut-in
     and rated speed, is constant between the rated and cut-out speed, and is zero for wind speeds above the cut-out
-    speed to avoid damage to the turbine (10.1088/1748-9326/aab211):
+    speed to avoid damage to the turbine :cite:p:`tobin_2018`:
 
     .. math::
 
         \begin{cases}
         0,  &  v < u_i \\
-        (v - u_i)^3 / (u_r - u_i)^3,  & u_i <= v < u_r \\
-        1, & u_r <= v < u_o \\
-        0, v >= u_o
-        \\end{cases}
+        (v - u_i)^3 / (u_r - u_i)^3,  & u_i ≤ v < u_r \\
+        1, & u_r ≤ v < u_o \\
+        0, & v ≥ u_o
+        \end{cases}
 
-    For non-standard air density, the wind speed is scaled using
-    :math:`v_n = v \\left( \frac{\rho}{\rho_0} \right)^{1/3}`.
+    For non-standard air density (:math:`\rho`), the wind speed is scaled using
+    :math:`v_n = v \left( \frac{\rho}{\rho_0} \right)^{1/3}`.
 
     Note that the temporal resolution of the wind speed time series has a significant influence on the results,
-    but percent changes in the wind power potential projections are similar across resolutions.
-    https://doi.org/10.1016/j.renene.2020.02.090
+    but percent changes in the wind power potential projections are similar across resolutions :cite:p:`chen_2020`.
+
+    References
+    ----------
+    :cite:cts:`chen_2020,tobin_2018`.
+
     """
     # Convert units
     cut_in = convert_units_to(cut_in, wind_speed)

From fcfe1a828585c30e8d55f6d2ac085a7bb5cb9339 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Fri, 8 Sep 2023 16:03:47 -0400
Subject: [PATCH 05/11] added indicators and their translation

---
 xclim/data/fr.json                    |  8 ++++++++
 xclim/indicators/atmos/_conversion.py | 26 ++++++++++++++++++++++++++
 2 files changed, 34 insertions(+)

diff --git a/xclim/data/fr.json b/xclim/data/fr.json
index 646a1f6a9..963a71593 100644
--- a/xclim/data/fr.json
+++ b/xclim/data/fr.json
@@ -792,6 +792,14 @@
     "long_name": "Indice forêt météo",
     "description": "Évaluation numérique de l'intensité du feu."
   },
+  "WIND_POWER_POTENTIAL": {
+    "title": "Potentiel de production éolienne",
+    "abstract": "Estimation à partir d'une loi de puissance semi-idéalisée de la production d'énergie éolienne selon la vitesse du vent."
+  },
+  "WIND_PROFILE": {
+    "title": "Profil du vent",
+    "abstract": "Estimation de la vitesse du vent à partir du vent à une hauteur de référence."
+  },
   "WIND_SPEED_FROM_VECTOR": {
     "title": "Vitesse et direction du vent à partir du vecteur vent",
     "abstract": "Calcul de la magnitude et la direction de la vitesse du vent à partir des deux composantes ouest-est et sud-nord."
diff --git a/xclim/indicators/atmos/_conversion.py b/xclim/indicators/atmos/_conversion.py
index d1771e88f..9988ea6d4 100644
--- a/xclim/indicators/atmos/_conversion.py
+++ b/xclim/indicators/atmos/_conversion.py
@@ -28,6 +28,8 @@
     "water_budget",
     "water_budget_from_tas",
     "wind_chill_index",
+    "wind_power_potential",
+    "wind_profile",
     "wind_speed_from_vector",
     "wind_vector_from_speed",
 ]
@@ -123,6 +125,30 @@ def cfcheck(self, **das):
     compute=indices.sfcwind_2_uas_vas,
 )
 
+wind_power_potential = Converter(
+    title="Wind power potential",
+    identifier="wind_power_potential",
+    units="W m-2",
+    long_name="Wind power potential",
+    description="Wind power potential using a semi-idealized turbine power curve using a cut_in speed of {cut_in}, "
+    "a rated speed of {rated}, and a cut_out speed of {cut_out}.",
+    cell_methods="",
+    abstract="Calculation of the wind power potential using a semi-idealized turbine power curve.",
+    compute=indices.wind_power_potential,
+)
+
+wind_profile = Converter(
+    title="Wind profile",
+    identifier="wind_profile",
+    var_name="wind_speed",
+    units="m s-1",
+    standard_name="wind_speed",
+    long_name="Wind speed at {h}",
+    description="Wind speed at {h} computed from the wind speed at {h_r} using a power law profile.",
+    cell_methods="",
+    abstract="Calculation of the wind speed at a given height from the wind speed at a reference height.",
+    compute=indices.wind_profile,
+)
 
 saturation_vapor_pressure = Converter(
     title="Saturation vapour pressure (e_sat)",

From 5f834a26c2e709af698e1b300075cb01946cfd14 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Fri, 8 Sep 2023 16:10:50 -0400
Subject: [PATCH 06/11] added smoke tests for the indicators

---
 tests/test_atmos.py                   | 8 ++++++++
 xclim/indicators/atmos/_conversion.py | 2 +-
 2 files changed, 9 insertions(+), 1 deletion(-)

diff --git a/tests/test_atmos.py b/tests/test_atmos.py
index 2828f92fc..deb36c214 100644
--- a/tests/test_atmos.py
+++ b/tests/test_atmos.py
@@ -253,6 +253,14 @@ def test_wind_chill_index(atmosds):
     )
 
 
+def test_wind_profile(atmosds):
+    atmos.wind_profile(wind_speed=atmosds.sfcWind, h_r="10 m", h="100 m", alpha=1 / 7)
+
+
+def test_wind_power_potential(atmosds):
+    atmos.wind_power_potential(wind_speed=atmosds.sfcWind)
+
+
 class TestDrynessIndex:
     def test_simple(self, atmosds):
         ds = atmosds.isel(location=3)
diff --git a/xclim/indicators/atmos/_conversion.py b/xclim/indicators/atmos/_conversion.py
index 9988ea6d4..bcb206de8 100644
--- a/xclim/indicators/atmos/_conversion.py
+++ b/xclim/indicators/atmos/_conversion.py
@@ -128,7 +128,7 @@ def cfcheck(self, **das):
 wind_power_potential = Converter(
     title="Wind power potential",
     identifier="wind_power_potential",
-    units="W m-2",
+    units="",
     long_name="Wind power potential",
     description="Wind power potential using a semi-idealized turbine power curve using a cut_in speed of {cut_in}, "
     "a rated speed of {rated}, and a cut_out speed of {cut_out}.",

From 3645862d67cea48af068f4804633612d25cc0706 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Mon, 11 Sep 2023 17:23:46 -0400
Subject: [PATCH 07/11] completed fr translations and added the new variables
 to the yaml.

---
 xclim/data/fr.json                    |  4 ++++
 xclim/data/variables.yml              | 12 ++++++++++++
 xclim/indicators/atmos/_conversion.py |  4 ++--
 3 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/xclim/data/fr.json b/xclim/data/fr.json
index 963a71593..96356ae89 100644
--- a/xclim/data/fr.json
+++ b/xclim/data/fr.json
@@ -793,10 +793,14 @@
     "description": "Évaluation numérique de l'intensité du feu."
   },
   "WIND_POWER_POTENTIAL": {
+    "long_name": "Potentiel de production éolienne",
+    "description": "Potentiel de production éolienne estimé par une courbe de puissance semi-idéalisée d'une turbine, paramétrée par une vitesse minimale de démarrage {cut_in}, la vitesse nominale {rated}, et la vitesse d'arrêt {cut_out}.",
     "title": "Potentiel de production éolienne",
     "abstract": "Estimation à partir d'une loi de puissance semi-idéalisée de la production d'énergie éolienne selon la vitesse du vent."
   },
   "WIND_PROFILE": {
+    "long_name": "Vitesse du vent à la hauteur {h}",
+    "description": "Vitesse du vent à une hauteur de {h} calculée de la vitesse à {h_r} par un profil de loi de puissance.",
     "title": "Profil du vent",
     "abstract": "Estimation de la vitesse du vent à partir du vent à une hauteur de référence."
   },
diff --git a/xclim/data/variables.yml b/xclim/data/variables.yml
index 4f14d5e3b..0a7553b0c 100644
--- a/xclim/data/variables.yml
+++ b/xclim/data/variables.yml
@@ -1,4 +1,10 @@
 variables:
+  air_density:
+    canonical_units: kg m-3
+    cell_methods: "time: mean"
+    description: Air density.
+    dimensions: "[density]"
+    standard_name: air_density
   areacello:
     canonical_units: m2
     cell_methods: "area: sum"
@@ -316,6 +322,12 @@ variables:
     description: Northward component of the wind velocity (at the surface).
     dimensions: "[speed]"
     standard_name: northward_wind
+  wind_speed:
+    canonical_units: m s-1
+    cell_methods: "time: mean"
+    description: Wind speed.
+    dimensions: "[speed]"
+    standard_name: wind_speed
   wsgsmax:
     cmip6: False
     canonical_units: m s-1
diff --git a/xclim/indicators/atmos/_conversion.py b/xclim/indicators/atmos/_conversion.py
index bcb206de8..6876755f2 100644
--- a/xclim/indicators/atmos/_conversion.py
+++ b/xclim/indicators/atmos/_conversion.py
@@ -143,8 +143,8 @@ def cfcheck(self, **das):
     var_name="wind_speed",
     units="m s-1",
     standard_name="wind_speed",
-    long_name="Wind speed at {h}",
-    description="Wind speed at {h} computed from the wind speed at {h_r} using a power law profile.",
+    long_name="Wind speed at height {h}",
+    description="Wind speed at a height of {h} computed from the wind speed at {h_r} using a power law profile.",
     cell_methods="",
     abstract="Calculation of the wind speed at a given height from the wind speed at a reference height.",
     compute=indices.wind_profile,

From dff6d3e14203cce3fafd06d5efa8c1789e53c169 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Mon, 11 Sep 2023 17:25:36 -0400
Subject: [PATCH 08/11] typo

---
 CHANGES.rst | 1 -
 1 file changed, 1 deletion(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 63a6b83bd..7bdc79cbc 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -2,7 +2,6 @@
 Changelog
 =========
 
-=======
 v0.46.0
 --------------------
 Contributors to this version: Éric Dupuis (:user:`coxipi`), David Huard (:user:`huard`)

From ce400e15c097cf18b0e1812af95fec546fa93074 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Mon, 18 Sep 2023 13:14:38 -0400
Subject: [PATCH 09/11] Changes based on review comments

---
 tests/test_indices.py        | 16 +++++++++++++---
 xclim/indices/_conversion.py | 16 +++++++++-------
 2 files changed, 22 insertions(+), 10 deletions(-)

diff --git a/tests/test_indices.py b/tests/test_indices.py
index abd606001..4e8356006 100644
--- a/tests/test_indices.py
+++ b/tests/test_indices.py
@@ -3637,11 +3637,21 @@ def test_simple(self, sfcWind_series):
 
 class TestWindPowerPotential:
     def test_simple(self, sfcWind_series):
-        a = [2, 6, 20, 30]
+        v = [2, 6, 20, 30]
         p = xci.wind_power_potential(
-            sfcWind_series(a, units="m/s"), cut_in="4 m/s", rated="8 m/s"
+            sfcWind_series(v, units="m/s"), cut_in="4 m/s", rated="8 m/s"
         )
-        np.testing.assert_allclose(p, [0, 0.5**3, 1, 0])
+        np.testing.assert_allclose(p, [0, (6**3 - 4**3) / (8**3 - 4**3), 1, 0])
+
+        # Test discontinuities at the default thresholds
+        v = np.array([3.5, 15])
+        a = sfcWind_series(v - 1e-7, units="m/s")
+        b = sfcWind_series(v + 1e-7, units="m/s")
+
+        pa = xci.wind_power_potential(a)
+        pb = xci.wind_power_potential(b)
+
+        np.testing.assert_array_almost_equal(pa, pb, decimal=6)
 
     # def test_benchmark(self, sfcWind_series, benchmark):
     #     a = np.random.rand(50000) * 30
diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py
index ddc72a7a3..125151324 100644
--- a/xclim/indices/_conversion.py
+++ b/xclim/indices/_conversion.py
@@ -2044,8 +2044,8 @@ def wind_profile(
     method : {"power_law"}
         Method to use. Currently only "power_law" is implemented.
     kwds : dict
-        Additional keyword arguments to pass to the method. For `power_law`, this is `alpha`, which takes a default
-        value of 1/7.
+        Additional keyword arguments to pass to the method. For power_law, this is alpha, which takes a default value
+        of 1/7, but is highly variable based on topography, surface cover and atmospheric stability.
 
     Notes
     -----
@@ -2086,7 +2086,7 @@ def wind_power_potential(
     rated: Quantified = "13 m/s",
     cut_out: Quantified = "25 m/s",
 ) -> xr.DataArray:
-    r"""Wind power potential estimated from a semi-idealized wind power production factor.
+    r"""Wind power potential estimated from an idealized wind power production factor.
 
     The actual power production of a wind farm can be estimated by multiplying its nominal (nameplate) capacity by the
     wind power potential, which depends on wind speed at the hub height, the turbine specifications and air density.
@@ -2118,7 +2118,7 @@ def wind_power_potential(
 
     Notes
     -----
-    This estimate of wind power production is based on a semi-idealized power curve with four wind regimes specified
+    This estimate of wind power production is based on an idealized power curve with four wind regimes specified
     by the cut-in wind speed (:math:`u_i`), the rated speed (:math:`u_r`) and the cut-out speed (:math:`u_o`).
     Power production is zero for wind speeds below the cut-in speed, increases cubically between the cut-in
     and rated speed, is constant between the rated and cut-out speed, and is zero for wind speeds above the cut-out
@@ -2128,7 +2128,7 @@ def wind_power_potential(
 
         \begin{cases}
         0,  &  v < u_i \\
-        (v - u_i)^3 / (u_r - u_i)^3,  & u_i ≤ v < u_r \\
+        (v^3 - u_i^3) / (u_r^3 - u_i^3),  & u_i ≤ v < u_r \\
         1, & u_r ≤ v < u_o \\
         0, & v ≥ u_o
         \end{cases}
@@ -2137,7 +2137,9 @@ def wind_power_potential(
     :math:`v_n = v \left( \frac{\rho}{\rho_0} \right)^{1/3}`.
 
     Note that the temporal resolution of the wind speed time series has a significant influence on the results,
-    but percent changes in the wind power potential projections are similar across resolutions :cite:p:`chen_2020`.
+    even when aggregated at longer time scales. Total annual power production will differ substantially if estimated
+    from hourly or daily wind speeds, due to the non-linearity of the production factor. Note however that percent
+    changes in the wind power potential projections are similar across resolutions :cite:p:`chen_2020`.
 
     References
     ----------
@@ -2169,7 +2171,7 @@ def _wind_power_factor(v, cut_in, rated, cut_out):
     if v < cut_in:
         return 0.0
     if v < rated:
-        return ((v - cut_in) / (rated - cut_in)) ** 3
+        return (v**3 - cut_in**3) / (rated**3 - cut_in**3)
     if v < cut_out:
         return 1.0
     return 0.0

From 5414af660e9e71386977d672613d881049568aa0 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Mon, 18 Sep 2023 13:18:02 -0400
Subject: [PATCH 10/11] remove benmark code

---
 tests/test_indices.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/tests/test_indices.py b/tests/test_indices.py
index 4e8356006..613669582 100644
--- a/tests/test_indices.py
+++ b/tests/test_indices.py
@@ -3652,8 +3652,3 @@ def test_simple(self, sfcWind_series):
         pb = xci.wind_power_potential(b)
 
         np.testing.assert_array_almost_equal(pa, pb, decimal=6)
-
-    # def test_benchmark(self, sfcWind_series, benchmark):
-    #     a = np.random.rand(50000) * 30
-    #     v = sfcWind_series(a, units="m/s")
-    #     benchmark(xci.wind_power_potential, v)

From 2dccc9394b69e952c5c32a19174025a28533cbe5 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Thu, 28 Sep 2023 11:22:14 -0400
Subject: [PATCH 11/11] Apply suggestions from review. Support sum in
 to_agg_units.

---
 tests/test_atmos.py          | 31 +++++++++++++++++++++++++++++--
 xclim/core/units.py          |  4 ++--
 xclim/indices/_conversion.py | 18 ++++++++++--------
 3 files changed, 41 insertions(+), 12 deletions(-)

diff --git a/tests/test_atmos.py b/tests/test_atmos.py
index deb36c214..5c86b8c29 100644
--- a/tests/test_atmos.py
+++ b/tests/test_atmos.py
@@ -254,11 +254,38 @@ def test_wind_chill_index(atmosds):
 
 
 def test_wind_profile(atmosds):
-    atmos.wind_profile(wind_speed=atmosds.sfcWind, h_r="10 m", h="100 m", alpha=1 / 7)
+    out = atmos.wind_profile(
+        wind_speed=atmosds.sfcWind, h_r="10 m", h="100 m", alpha=1 / 7
+    )
+    assert out.attrs["units"] == "m s-1"
+    assert (out > atmosds.sfcWind).all()
 
 
 def test_wind_power_potential(atmosds):
-    atmos.wind_power_potential(wind_speed=atmosds.sfcWind)
+    out = atmos.wind_power_potential(wind_speed=atmosds.sfcWind)
+    assert out.attrs["units"] == ""
+    assert (out >= 0).all()
+    assert (out <= 1).all()
+
+
+def test_wind_power_potential_from_3h_series():
+    """Test a typical computation workflow from 3-hourly time series to daily power production in MWh."""
+    from xclim.core.units import convert_units_to
+    from xclim.indices.generic import select_resample_op
+    from xclim.testing.helpers import test_timeseries
+
+    w = test_timeseries(
+        np.ones(96) * 15, variable="sfcWind", start="7/1/2000", units="m s-1", freq="3H"
+    )
+    out = atmos.wind_power_potential(wind_speed=w)
+
+    # Multiply with nominal capacity
+    power = out * 100
+    power.attrs["units"] = "MW"
+    annual_power = convert_units_to(
+        select_resample_op(power, op="sum", freq="D"), "MWh"
+    )
+    assert (annual_power == 100 * 24).all()
 
 
 class TestDrynessIndex:
diff --git a/xclim/core/units.py b/xclim/core/units.py
index b79fddb08..474394afb 100644
--- a/xclim/core/units.py
+++ b/xclim/core/units.py
@@ -579,7 +579,7 @@ def to_agg_units(
             units="", is_dayofyear=np.int32(1), calendar=get_calendar(orig)
         )
 
-    elif op in ["count", "integral"]:
+    elif op in ["count", "integral", "sum"]:
         m, freq_u_raw = infer_sampling_units(orig[dim])
         orig_u = str2pint(orig.units)
         freq_u = str2pint(freq_u_raw)
@@ -587,7 +587,7 @@ def to_agg_units(
 
         if op == "count":
             out.attrs["units"] = freq_u_raw
-        elif op == "integral":
+        elif op in ["integral", "sum"]:
             out.attrs["units"] = pint2cfunits(orig_u * freq_u)
     else:
         raise ValueError(
diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py
index 125151324..6e571e082 100644
--- a/xclim/indices/_conversion.py
+++ b/xclim/indices/_conversion.py
@@ -5,7 +5,7 @@
 
 import numpy as np
 import xarray as xr
-from numba import float32, float64, njit, vectorize  # noqa
+from numba import float32, float64, vectorize  # noqa
 
 from xclim.core.calendar import date_range, datetime_to_decimal_year
 from xclim.core.units import (
@@ -2094,9 +2094,7 @@ def wind_power_potential(
     Parameters
     ----------
     wind_speed : xarray.DataArray
-        Wind speed at the hub height. Use the `wind_profile` function to estimate from the surface wind speed. Note
-        that the temporal resolution of wind time series has a significant influence on the results. Mean daily wind
-        speeds yield lower values than hourly wind speeds.
+        Wind speed at the hub height. Use the `wind_profile` function to estimate from the surface wind speed.
     air_density: xarray.DataArray
         Air density at the hub height. Defaults to 1.225 kg/m³. This is worth changing if applying in cold or
         mountainous regions with non-standard air density.
@@ -2136,10 +2134,14 @@ def wind_power_potential(
     For non-standard air density (:math:`\rho`), the wind speed is scaled using
     :math:`v_n = v \left( \frac{\rho}{\rho_0} \right)^{1/3}`.
 
-    Note that the temporal resolution of the wind speed time series has a significant influence on the results,
-    even when aggregated at longer time scales. Total annual power production will differ substantially if estimated
-    from hourly or daily wind speeds, due to the non-linearity of the production factor. Note however that percent
-    changes in the wind power potential projections are similar across resolutions :cite:p:`chen_2020`.
+    The temporal resolution of wind time series has a significant influence on the results: mean daily wind
+    speeds yield lower values than hourly wind speeds. Note however that percent changes in the wind power potential
+    climate projections are similar across resolutions :cite:p:`chen_2020`.
+
+    To compute the power production, multiply the power production factor by the nominal
+    turbine capacity (e.g. 100), set the units attribute (e.g. "MW"), resample and sum with
+    `xclim.indices.generic.select_resample_op(power, op="sum", freq="D")`, then convert to
+    the desired units (e.g. "MWh") using `xclim.core.units.convert_units_to`.
 
     References
     ----------