diff --git a/CHANGES.rst b/CHANGES.rst
index 7a42043bb..bf01c1f80 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -13,6 +13,16 @@ New features and enhancements
 Breaking changes
 ^^^^^^^^^^^^^^^^
 * `pint` has been pinned below v0.24 until `xclim` can be updated to support the latest version. (:issue:`1771`, :pull:`1772`).
+* Calendar utilities that have an equivalent in `xarray` have been deprecated and will be removed in `xclim` v0.51.0. (:issue:`1010`, :pull:`1761`). This concerns the following members of ``xclim.core.calendar``:
+    - ``convert_calendar`` : Use ``Dataset.convert_calendar``, ``DataArray.convert_calendar`` or ``xr.coding.calendar_ops.convert_calendar``  instead.
+        + If your code passes ``target`` as an array, first convert the source to the target's calendar and then reindex the result to ``target``.
+        + If you were using the ``doy=True`` option, replace it with ``xc.core.calendar.convert_doy(source, target_cal).convert_calendar(target_cal)``.
+        + ``"default"`` is no longer a valid calendar name for any xclim functions and will not be returned by ``get_calendar``. Xarray has a ``use_cftime`` argument, xclim exposes it when the distinction is needed.
+    - ``date_range`` : Use ``xarray.date_range`` instead.
+    - ``date_range_like``: Use ``xarray.date_range_like`` instead.
+    - ``interp_calendar`` : Use ``Dataset.interp_calendar`` or ``xarray.coding.calendar_ops.interp_calendar`` instead.
+    - ``days_in_year`` : Use ``xarray.coding.calendar_ops._days_in_year`` instead.
+    - ``datetime_to_decimal_year`` : Use ``xarray.coding.calendar_ops._datetime_to_decimal_year`` instead.
 
 Internal changes
 ^^^^^^^^^^^^^^^^
diff --git a/tests/test_calendar.py b/tests/test_calendar.py
index 9e9af5630..31a3849a0 100644
--- a/tests/test_calendar.py
+++ b/tests/test_calendar.py
@@ -16,16 +16,11 @@
     common_calendar,
     compare_offsets,
     construct_offset,
-    convert_calendar,
     convert_doy,
-    date_range,
-    datetime_to_decimal_year,
-    days_in_year,
     days_since_to_doy,
     doy_to_days_since,
     ensure_cftime_array,
     get_calendar,
-    interp_calendar,
     max_doy,
     parse_offset,
     percentile_doy,
@@ -236,7 +231,7 @@ def test_adjust_doy_366_to_360():
             "360_day",
             360,
         ),
-        (("NRCANdaily", "nrcan_canada_daily_pr_1990.nc"), "default", 366),
+        (("NRCANdaily", "nrcan_canada_daily_pr_1990.nc"), "proleptic_gregorian", 366),
     ],
 )
 def test_get_calendar(file, cal, maxdoy, open_dataset):
@@ -249,8 +244,8 @@ def test_get_calendar(file, cal, maxdoy, open_dataset):
 @pytest.mark.parametrize(
     "obj,cal",
     [
-        ([pd.Timestamp.now()], "default"),
-        (pd.Timestamp.now(), "default"),
+        ([pd.Timestamp.now()], "standard"),
+        (pd.Timestamp.now(), "standard"),
         (cftime.DatetimeAllLeap(2000, 1, 1), "all_leap"),
         (np.array([cftime.DatetimeNoLeap(2000, 1, 1)]), "noleap"),
         (xr.cftime_range("2000-01-01", periods=4, freq="D"), "standard"),
@@ -266,158 +261,6 @@ def test_get_calendar_errors(obj):
         get_calendar(obj)
 
 
-@pytest.mark.parametrize(
-    "source,target,target_as_str,freq",
-    [
-        ("standard", "noleap", True, "D"),
-        ("noleap", "default", True, "D"),
-        ("noleap", "all_leap", False, "D"),
-        ("proleptic_gregorian", "noleap", False, "4h"),
-        ("default", "noleap", True, "4h"),
-    ],
-)
-def test_convert_calendar(source, target, target_as_str, freq):
-    src = xr.DataArray(
-        date_range("2004-01-01", "2004-12-31", freq=freq, calendar=source),
-        dims=("time",),
-        name="time",
-    )
-    da_src = xr.DataArray(
-        np.linspace(0, 1, src.size), dims=("time",), coords={"time": src}
-    )
-    tgt = xr.DataArray(
-        date_range("2004-01-01", "2004-12-31", freq=freq, calendar=target),
-        dims=("time",),
-        name="time",
-    )
-
-    conv = convert_calendar(da_src, target if target_as_str else tgt)
-
-    assert get_calendar(conv) == target
-
-    if target_as_str and max_doy[source] < max_doy[target]:
-        assert conv.size == src.size
-    elif not target_as_str:
-        assert conv.size == tgt.size
-
-        assert conv.isnull().sum() == max(max_doy[target] - max_doy[source], 0)
-
-
-@pytest.mark.parametrize(
-    "source,target,freq",
-    [
-        ("standard", "360_day", "D"),
-        ("360_day", "default", "D"),
-        ("proleptic_gregorian", "360_day", "4h"),
-    ],
-)
-@pytest.mark.parametrize("align_on", ["date", "year"])
-def test_convert_calendar_360_days(source, target, freq, align_on):
-    src = xr.DataArray(
-        date_range("2004-01-01", "2004-12-30", freq=freq, calendar=source),
-        dims=("time",),
-        name="time",
-    )
-    da_src = xr.DataArray(
-        np.linspace(0, 1, src.size), dims=("time",), coords={"time": src}
-    )
-
-    conv = convert_calendar(da_src, target, align_on=align_on)
-
-    assert get_calendar(conv) == target
-
-    if align_on == "date":
-        np.testing.assert_array_equal(
-            conv.time.resample(time="ME").last().dt.day,
-            [30, 29, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30],
-        )
-    elif target == "360_day":
-        np.testing.assert_array_equal(
-            conv.time.resample(time="ME").last().dt.day,
-            [30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 29],
-        )
-    else:
-        np.testing.assert_array_equal(
-            conv.time.resample(time="ME").last().dt.day,
-            [30, 29, 30, 30, 31, 30, 30, 31, 30, 31, 29, 31],
-        )
-    if source == "360_day" and align_on == "year":
-        assert conv.size == 360 if freq == "D" else 360 * 4
-    else:
-        assert conv.size == 359 if freq == "D" else 359 * 4
-
-
-def test_convert_calendar_360_days_random():
-    da_std = xr.DataArray(
-        np.linspace(0, 1, 366 * 2),
-        dims=("time",),
-        coords={
-            "time": date_range(
-                "2004-01-01", "2004-12-31T23:59:59", freq="12h", calendar="default"
-            )
-        },
-    )
-    da_360 = xr.DataArray(
-        np.linspace(0, 1, 360 * 2),
-        dims=("time",),
-        coords={
-            "time": date_range(
-                "2004-01-01", "2004-12-30T23:59:59", freq="12h", calendar="360_day"
-            )
-        },
-    )
-
-    conv = convert_calendar(da_std, "360_day", align_on="random")
-    assert get_calendar(conv) == "360_day"
-    assert conv.size == 720
-    conv2 = convert_calendar(da_std, "360_day", align_on="random")
-    assert (conv != conv2).any()
-
-    conv = convert_calendar(da_360, "default", align_on="random")
-    assert get_calendar(conv) == "default"
-    assert conv.size == 720
-    assert np.datetime64("2004-02-29") not in conv.time
-    conv2 = convert_calendar(da_360, "default", align_on="random")
-    assert (conv2 != conv).any()
-
-    conv = convert_calendar(da_360, "noleap", align_on="random", missing=np.NaN)
-    conv = conv.where(conv.isnull(), drop=True)
-    nandoys = conv.time.dt.dayofyear[::2]
-    assert all(nandoys < np.array([74, 147, 220, 293, 366]))
-    assert all(nandoys > np.array([0, 73, 146, 219, 292]))
-
-
-@pytest.mark.parametrize(
-    "source,target,freq",
-    [
-        ("standard", "noleap", "D"),
-        ("noleap", "default", "4h"),
-        ("noleap", "all_leap", "ME"),
-        ("360_day", "noleap", "D"),
-        ("noleap", "360_day", "D"),
-    ],
-)
-def test_convert_calendar_missing(source, target, freq):
-    src = xr.DataArray(
-        date_range(
-            "2004-01-01",
-            "2004-12-31" if source != "360_day" else "2004-12-30",
-            freq=freq,
-            calendar=source,
-        ),
-        dims=("time",),
-        name="time",
-    )
-    da_src = xr.DataArray(
-        np.linspace(0, 1, src.size), dims=("time",), coords={"time": src}
-    )
-    out = convert_calendar(da_src, target, missing=0, align_on="date")
-    assert xr.infer_freq(out.time) == freq
-    if source == "360_day":
-        assert out.time[-1].dt.day == 31
-        assert out[-1] == 0
-
-
 def test_convert_calendar_and_doy():
     doy = xr.DataArray(
         [31, 32, 336, 364.23, 365],
@@ -427,68 +270,42 @@ def test_convert_calendar_and_doy():
         },
         attrs={"is_dayofyear": 1, "calendar": "noleap"},
     )
-    out = convert_calendar(doy, "360_day", align_on="date", doy=True)
+    out = convert_doy(doy, target_cal="360_day").convert_calendar(
+        "360_day", align_on="date"
+    )
+    # out = convert_calendar(doy, "360_day", align_on="date", doy=True)
     np.testing.assert_allclose(
         out, [30.575342, 31.561644, 331.39726, 359.240548, 360.0]
     )
     assert out.time.dt.calendar == "360_day"
-    out = convert_calendar(doy, "360_day", align_on="date", doy="date")
+    out = convert_doy(doy, target_cal="360_day", align_on="date").convert_calendar(
+        "360_day", align_on="date"
+    )
     np.testing.assert_array_equal(out, [np.NaN, 31, 332, 360.23, np.NaN])
     assert out.time.dt.calendar == "360_day"
 
 
-@pytest.mark.parametrize(
-    "source,target",
-    [
-        ("standard", "noleap"),
-        ("noleap", "default"),
-        ("standard", "360_day"),
-        ("360_day", "standard"),
-        ("noleap", "all_leap"),
-        ("360_day", "noleap"),
-    ],
-)
-def test_interp_calendar(source, target):
-    src = xr.DataArray(
-        date_range("2004-01-01", "2004-07-30", freq="D", calendar=source),
-        dims=("time",),
-        name="time",
-    )
-    tgt = xr.DataArray(
-        date_range("2004-01-01", "2004-07-30", freq="D", calendar=target),
-        dims=("time",),
-        name="time",
-    )
-    da_src = xr.DataArray(
-        np.linspace(0, 1, src.size), dims=("time",), coords={"time": src}
-    )
-    conv = interp_calendar(da_src, tgt)
-
-    assert conv.size == tgt.size
-    assert get_calendar(conv) == target
-
-    np.testing.assert_almost_equal(conv.max(), 1, 2)
-    assert conv.min() == 0
-
-
 @pytest.mark.parametrize(
     "inp,calout",
     [
         (
             xr.DataArray(
-                date_range("2004-01-01", "2004-01-10", freq="D"),
+                xr.date_range("2004-01-01", "2004-01-10", freq="D"),
                 dims=("time",),
                 name="time",
             ),
             "standard",
         ),
-        (date_range("2004-01-01", "2004-01-10", freq="D"), "standard"),
+        (xr.date_range("2004-01-01", "2004-01-10", freq="D"), "standard"),
         (
-            xr.DataArray(date_range("2004-01-01", "2004-01-10", freq="D")).values,
+            xr.DataArray(xr.date_range("2004-01-01", "2004-01-10", freq="D")).values,
             "standard",
         ),
-        (date_range("2004-01-01", "2004-01-10", freq="D").values, "standard"),
-        (date_range("2004-01-01", "2004-01-10", freq="D", calendar="julian"), "julian"),
+        (xr.date_range("2004-01-01", "2004-01-10", freq="D").values, "standard"),
+        (
+            xr.date_range("2004-01-01", "2004-01-10", freq="D", calendar="julian"),
+            "julian",
+        ),
     ],
 )
 def test_ensure_cftime_array(inp, calout):
@@ -496,45 +313,6 @@ def test_ensure_cftime_array(inp, calout):
     assert get_calendar(out) == calout
 
 
-@pytest.mark.parametrize(
-    "year,calendar,exp",
-    [
-        (2004, "standard", 366),
-        (2004, "noleap", 365),
-        (2004, "all_leap", 366),
-        (1500, "default", 365),
-        (1500, "standard", 366),
-        (1500, "proleptic_gregorian", 365),
-        (2030, "360_day", 360),
-    ],
-)
-def test_days_in_year(year, calendar, exp):
-    assert days_in_year(year, calendar) == exp
-
-
-@pytest.mark.parametrize(
-    "source_cal, exp180",
-    [
-        ("standard", 0.49180328),
-        ("default", 0.49180328),
-        ("noleap", 0.49315068),
-        ("all_leap", 0.49180328),
-        ("360_day", 0.5),
-        (None, 0.49180328),
-    ],
-)
-def test_datetime_to_decimal_year(source_cal, exp180):
-    times = xr.DataArray(
-        date_range(
-            "2004-01-01", "2004-12-30", freq="D", calendar=source_cal or "default"
-        ),
-        dims=("time",),
-        name="time",
-    )
-    decy = datetime_to_decimal_year(times, calendar=source_cal)
-    np.testing.assert_almost_equal(decy[180] - 2004, exp180)
-
-
 def test_clim_mean_doy(tas_series):
     arr = tas_series(np.ones(365 * 10))
     mean, stddev = climatological_mean_doy(arr, window=1)
@@ -552,12 +330,12 @@ def test_clim_mean_doy(tas_series):
 
 def test_doy_to_days_since():
     # simple test
-    time = date_range("2020-07-01", "2022-07-01", freq="YS-JUL")
+    time = xr.date_range("2020-07-01", "2022-07-01", freq="YS-JUL")
     da = xr.DataArray(
         [190, 360, 3],
         dims=("time",),
         coords={"time": time},
-        attrs={"is_dayofyear": 1, "calendar": "default"},
+        attrs={"is_dayofyear": 1, "calendar": "standard"},
     )
 
     out = doy_to_days_since(da)
@@ -583,13 +361,13 @@ def test_doy_to_days_since():
     xr.testing.assert_identical(da, da2)
 
     # with start
-    time = date_range("2020-12-31", "2022-12-31", freq="YE")
+    time = xr.date_range("2020-12-31", "2022-12-31", freq="YE")
     da = xr.DataArray(
         [190, 360, 3],
         dims=("time",),
         coords={"time": time},
         name="da",
-        attrs={"is_dayofyear": 1, "calendar": "default"},
+        attrs={"is_dayofyear": 1, "calendar": "proleptic_gregorian"},
     )
 
     out = doy_to_days_since(da, start="01-02")
@@ -600,13 +378,13 @@ def test_doy_to_days_since():
     xr.testing.assert_identical(da, da2)
 
     # finer freq
-    time = date_range("2020-01-01", "2020-03-01", freq="MS")
+    time = xr.date_range("2020-01-01", "2020-03-01", freq="MS")
     da = xr.DataArray(
         [15, 33, 66],
         dims=("time",),
         coords={"time": time},
         name="da",
-        attrs={"is_dayofyear": 1, "calendar": "default"},
+        attrs={"is_dayofyear": 1, "calendar": "proleptic_gregorian"},
     )
 
     out = doy_to_days_since(da)
diff --git a/tests/test_generic.py b/tests/test_generic.py
index 7ec0771cd..c1eab6bf2 100644
--- a/tests/test_generic.py
+++ b/tests/test_generic.py
@@ -6,7 +6,7 @@
 import pytest
 import xarray as xr
 
-from xclim.core.calendar import date_range, doy_to_days_since, select_time
+from xclim.core.calendar import doy_to_days_since, select_time
 from xclim.indices import generic
 
 K2C = 273.15
@@ -108,8 +108,12 @@ def test_doyminmax(self, q_series):
 class TestAggregateBetweenDates:
     def test_calendars(self):
         # generate test DataArray
-        time_std = date_range("1991-07-01", "1993-06-30", freq="D", calendar="standard")
-        time_365 = date_range("1991-07-01", "1993-06-30", freq="D", calendar="noleap")
+        time_std = xr.date_range(
+            "1991-07-01", "1993-06-30", freq="D", calendar="standard"
+        )
+        time_365 = xr.date_range(
+            "1991-07-01", "1993-06-30", freq="D", calendar="noleap"
+        )
         data_std = xr.DataArray(
             np.ones((time_std.size, 4)),
             dims=("time", "lon"),
@@ -159,13 +163,15 @@ def test_calendars(self):
 
     def test_time_length(self):
         # generate test DataArray
-        time_data = date_range(
+        time_data = xr.date_range(
             "1991-01-01", "1993-12-31", freq="D", calendar="standard"
         )
-        time_start = date_range(
+        time_start = xr.date_range(
             "1990-01-01", "1992-12-31", freq="D", calendar="standard"
         )
-        time_end = date_range("1991-01-01", "1993-12-31", freq="D", calendar="standard")
+        time_end = xr.date_range(
+            "1991-01-01", "1993-12-31", freq="D", calendar="standard"
+        )
         data = xr.DataArray(
             np.ones((time_data.size, 4)),
             dims=("time", "lon"),
@@ -206,7 +212,7 @@ def test_time_length(self):
 
     def test_frequency(self):
         # generate test DataArray
-        time_data = date_range(
+        time_data = xr.date_range(
             "1991-01-01", "1992-05-31", freq="D", calendar="standard"
         )
         data = xr.DataArray(
@@ -280,7 +286,7 @@ def test_frequency(self):
 
     def test_day_of_year_strings(self):
         # generate test DataArray
-        time_data = date_range(
+        time_data = xr.date_range(
             "1990-08-01", "1995-06-01", freq="D", calendar="standard"
         )
         data = xr.DataArray(
@@ -519,7 +525,12 @@ def test_last_occurrence(self, tas_series, op, constrain, expected, should_fail)
 class TestTimeSelection:
     @staticmethod
     def series(start, end, calendar):
-        time = date_range(start, end, calendar=calendar)
+        time = xr.date_range(
+            start,
+            end,
+            calendar=calendar.replace("default", "proleptic_gregorian"),
+            use_cftime=(calendar != "default"),
+        )
         return xr.DataArray([1] * time.size, dims=("time",), coords={"time": time})
 
     def test_select_time_month(self):
diff --git a/tests/test_helpers.py b/tests/test_helpers.py
index 2cb66ee39..699cb3b6b 100644
--- a/tests/test_helpers.py
+++ b/tests/test_helpers.py
@@ -5,7 +5,6 @@
 import pytest
 import xarray as xr
 
-from xclim.core.calendar import date_range
 from xclim.core.units import convert_units_to
 from xclim.indices import helpers
 
@@ -54,7 +53,7 @@ def test_extraterrestrial_radiation(method):
 
 @pytest.mark.parametrize("method", ["spencer", "simple"])
 def test_day_lengths(method):
-    time_data = date_range("1992-12-01", "1994-01-01", freq="D", calendar="standard")
+    time_data = xr.date_range("1992-12-01", "1994-01-01", freq="D", calendar="standard")
     data = xr.DataArray(
         np.ones((time_data.size, 7)),
         dims=("time", "lat"),
@@ -66,12 +65,12 @@ def test_day_lengths(method):
 
     events = dict(
         solstice=[
-            ["1992-12-21", [[18.49, 15.43, 13.93, 12.0, 10.07, 8.57, 5.51]]],
-            ["1993-06-21", [[5.51, 8.57, 10.07, 12.0, 13.93, 15.43, 18.49]]],
-            ["1993-12-21", [[18.49, 15.43, 13.93, 12.0, 10.07, 8.57, 5.51]]],
+            ["1992-12-21", [18.49, 15.43, 13.93, 12.0, 10.07, 8.57, 5.51]],
+            ["1993-06-21", [5.51, 8.57, 10.07, 12.0, 13.93, 15.43, 18.49]],
+            ["1993-12-21", [18.49, 15.43, 13.93, 12.0, 10.07, 8.57, 5.51]],
         ],
         equinox=[
-            ["1993-03-20", [[12] * 7]]
+            ["1993-03-20", [12] * 7]
         ],  # True equinox on 1993-03-20 at 14:41 GMT. Some relative tolerance is needed.
     )
 
diff --git a/tests/test_indices.py b/tests/test_indices.py
index 2babd4419..9629d22d5 100644
--- a/tests/test_indices.py
+++ b/tests/test_indices.py
@@ -21,7 +21,7 @@
 import xarray as xr
 
 from xclim import indices as xci
-from xclim.core.calendar import convert_calendar, date_range, percentile_doy
+from xclim.core.calendar import percentile_doy
 from xclim.core.options import set_options
 from xclim.core.units import ValidationError, convert_units_to, units
 
@@ -258,7 +258,7 @@ def test_corn_heat_units(self, tasmin_series, tasmax_series):
         ],
     )
     def test_bedd(self, method, end_date, deg_days, max_deg_days):
-        time_data = date_range(
+        time_data = xr.date_range(
             "1992-01-01", "1995-06-01", freq="D", calendar="standard"
         )
         tn = xr.DataArray(
@@ -596,9 +596,8 @@ def test_standardized_precipitation_index(
     ):
         ds = open_dataset("sdba/CanESM2_1950-2100.nc").isel(location=1)
         if freq == "D":
-            ds = convert_calendar(
-                ds, "366_day", missing=np.NaN
-            )  # to compare with ``climate_indices``
+            ds = ds.convert_calendar("366_day", missing=np.NaN)
+            # to compare with ``climate_indices``
         pr = ds.pr.sel(time=slice("1998", "2000"))
         pr_cal = ds.pr.sel(time=slice("1950", "1980"))
         fitkwargs = {}
diff --git a/tests/test_missing.py b/tests/test_missing.py
index 197a8446d..b6523de08 100644
--- a/tests/test_missing.py
+++ b/tests/test_missing.py
@@ -6,7 +6,6 @@
 import xarray as xr
 
 from xclim.core import missing
-from xclim.core.calendar import convert_calendar
 
 K2C = 273.15
 
@@ -116,10 +115,10 @@ def test_month(self, tasmin_series):
         miss = missing.missing_any(ts, freq="YS", month=[7, 8])
         np.testing.assert_equal(miss, [False])
 
-    @pytest.mark.parametrize("calendar", ("default", "noleap", "360_day"))
+    @pytest.mark.parametrize("calendar", ("proleptic_gregorian", "noleap", "360_day"))
     def test_season(self, tasmin_series, calendar):
         ts = tasmin_series(np.zeros(360))
-        ts = convert_calendar(ts, calendar, missing=0, align_on="date")
+        ts = ts.convert_calendar(calendar, missing=0, align_on="date")
 
         miss = missing.missing_any(ts, freq="YS", season="MAM")
         np.testing.assert_equal(miss, [False])
diff --git a/tests/test_sdba/test_processing.py b/tests/test_sdba/test_processing.py
index 0dca69583..b0888e3ff 100644
--- a/tests/test_sdba/test_processing.py
+++ b/tests/test_sdba/test_processing.py
@@ -5,7 +5,6 @@
 import pytest
 import xarray as xr
 
-from xclim.core.calendar import date_range
 from xclim.core.units import units
 from xclim.sdba.adjustment import EmpiricalQuantileMapping
 from xclim.sdba.base import Grouper
@@ -188,8 +187,8 @@ def test_reordering():
 
 def test_reordering_with_window():
     time = list(
-        date_range("2000-01-01", "2000-01-04", freq="D", calendar="noleap")
-    ) + list(date_range("2001-01-01", "2001-01-04", freq="D", calendar="noleap"))
+        xr.date_range("2000-01-01", "2000-01-04", freq="D", calendar="noleap")
+    ) + list(xr.date_range("2001-01-01", "2001-01-04", freq="D", calendar="noleap"))
 
     x = xr.DataArray(
         np.arange(1, 9, 1),
diff --git a/xclim/core/bootstrapping.py b/xclim/core/bootstrapping.py
index d38764f2e..280ac316d 100644
--- a/xclim/core/bootstrapping.py
+++ b/xclim/core/bootstrapping.py
@@ -14,7 +14,7 @@
 
 import xclim.core.utils
 
-from .calendar import convert_calendar, parse_offset, percentile_doy
+from .calendar import parse_offset, percentile_doy
 
 BOOTSTRAP_DIM = "_bootstrap"
 
@@ -261,10 +261,10 @@ def build_bootstrap_year_da(
         elif len(source[dim]) == len(bloc):
             out_view.loc[{dim: bloc}] = source.data
         elif len(bloc) == 365:
-            out_view.loc[{dim: bloc}] = convert_calendar(source, "365_day").data
+            out_view.loc[{dim: bloc}] = source.convert_calendar("noleap").data
         elif len(bloc) == 366:
-            out_view.loc[{dim: bloc}] = convert_calendar(
-                source, "366_day", missing=np.NAN
+            out_view.loc[{dim: bloc}] = source.convert_calendar(
+                "366_day", missing=np.NAN
             ).data
         elif len(bloc) < 365:
             # 360 days calendar case or anchored years for both source[dim] and bloc case
diff --git a/xclim/core/calendar.py b/xclim/core/calendar.py
index 3927b1075..4bc1136a4 100644
--- a/xclim/core/calendar.py
+++ b/xclim/core/calendar.py
@@ -10,6 +10,7 @@
 import datetime as pydt
 from collections.abc import Sequence
 from typing import Any, TypeVar
+from warnings import warn
 
 import cftime
 import numpy as np
@@ -55,13 +56,10 @@
     "uniform_calendars",
     "unstack_periods",
     "within_bnds_doy",
-    "yearly_interpolated_doy",
-    "yearly_random_doy",
 ]
 
 # Maximum day of year in each calendar.
 max_doy = {
-    "default": 366,
     "standard": 366,
     "gregorian": 366,
     "proleptic_gregorian": 366,
@@ -74,7 +72,7 @@
 }
 
 # Some xclim.core.utils functions made accessible here for backwards compatibility reasons.
-datetime_classes = {"default": pydt.datetime, **cftime._cftime.DATE_TYPES}  # noqa
+datetime_classes = cftime._cftime.DATE_TYPES
 
 # Names of calendars that have the same number of days for all years
 uniform_calendars = ("noleap", "all_leap", "365_day", "366_day", "360_day")
@@ -83,13 +81,30 @@
 DataType = TypeVar("DataType", xr.DataArray, xr.Dataset)
 
 
-def days_in_year(year: int, calendar: str = "default") -> int:
-    """Return the number of days in the input year according to the input calendar."""
-    return (
-        (datetime_classes[calendar](year + 1, 1, 1) - pydt.timedelta(days=1))
-        .timetuple()
-        .tm_yday
+def _get_usecf_and_warn(calendar: str, xcfunc: str, xrfunc: str):
+    if calendar == "default":
+        calendar = "standard"
+        use_cftime = False
+        msg = " and use use_cftime=False instead of calendar='default' to get numpy objects."
+    else:
+        use_cftime = None
+        msg = ""
+    warn(
+        f"`xclim` function {xcfunc} is deprecated in favour of {xrfunc} and will be removed in v0.51.0. Please adjust your script{msg}.",
+        FutureWarning,
+    )
+    return calendar, use_cftime
+
+
+def days_in_year(year: int, calendar: str = "proleptic_gregorian") -> int:
+    """Deprecated : use :py:func:`xarray.coding.calendar_ops._days_in_year` instead. Passing use_cftime=False instead of calendar='default'.
+
+    Return the number of days in the input year according to the input calendar.
+    """
+    calendar, usecf = _get_usecf_and_warn(
+        calendar, "days_in_year", "xarray.coding.calendar_ops._days_in_year"
     )
+    return xr.coding.calendar_ops._days_in_year(year, calendar, use_cftime=usecf)
 
 
 def doy_from_string(doy: DayOfYearStr, year: int, calendar: str) -> int:
@@ -98,53 +113,17 @@ def doy_from_string(doy: DayOfYearStr, year: int, calendar: str) -> int:
     return datetime_classes[calendar](year, int(MM), int(DD)).timetuple().tm_yday
 
 
-def date_range(
-    *args, calendar: str = "default", **kwargs
-) -> pd.DatetimeIndex | CFTimeIndex:
-    """Wrap a Pandas date_range object.
-
-    Uses pd.date_range (if calendar == 'default') or xr.cftime_range (otherwise).
-    """
-    if calendar == "default":
-        return pd.date_range(*args, **kwargs)
-    return xr.cftime_range(*args, calendar=calendar, **kwargs)
+def date_range(*args, **kwargs) -> pd.DatetimeIndex | CFTimeIndex:
+    """Deprecated : use :py:func:`xarray.date_range` instead. Passing use_cftime=False instead of calendar='default'.
 
+    Wrap a Pandas date_range object.
 
-def yearly_interpolated_doy(
-    time: pd.DatetimeIndex | CFTimeIndex, source_calendar: str, target_calendar: str
-):
-    """Return the nearest day in the target calendar of the corresponding "decimal year" in the source calendar."""
-    yr = int(time.dt.year[0])
-    return np.round(
-        days_in_year(yr, target_calendar)
-        * time.dt.dayofyear
-        / days_in_year(yr, source_calendar)
-    ).astype(int)
-
-
-def yearly_random_doy(
-    time: pd.DatetimeIndex | CFTimeIndex,
-    rng: np.random.Generator,
-    source_calendar: str,
-    target_calendar: str,
-):
-    """Return a day of year in the new calendar.
-
-    Removes Feb 29th and five other days chosen randomly within five sections of 72 days.
+    Uses pd.date_range (if calendar == 'default') or xr.cftime_range (otherwise).
     """
-    yr = int(time.dt.year[0])
-    new_doy = np.arange(360) + 1
-    rm_idx = rng.integers(0, 72, 5) + (np.arange(5) * 72)
-    if source_calendar == "360_day":
-        for idx in rm_idx:
-            new_doy[idx + 1 :] = new_doy[idx + 1 :] + 1
-        if days_in_year(yr, target_calendar) == 366:
-            new_doy[new_doy >= 60] = new_doy[new_doy >= 60] + 1
-    elif target_calendar == "360_day":
-        new_doy = np.insert(new_doy, rm_idx - np.arange(5), -1)
-        if days_in_year(yr, source_calendar) == 366:
-            new_doy = np.insert(new_doy, 60, -1)
-    return new_doy[time.dt.dayofyear - 1]
+    calendar, usecf = _get_usecf_and_warn(
+        kwargs.pop("calendar", "default"), "date_range", "xarray.date_range"
+    )
+    return xr.date_range(*args, calendar=calendar, use_cftime=usecf, **kwargs)
 
 
 def get_calendar(obj: Any, dim: str = "time") -> str:
@@ -169,21 +148,18 @@ def get_calendar(obj: Any, dim: str = "time") -> str:
     Returns
     -------
     str
-      The cftime calendar name or "default" when the data is using numpy's or python's datetime types.
+      The Climate and Forecasting (CF) calendar name.
       Will always return "standard" instead of "gregorian", following CF conventions 1.9.
     """
     if isinstance(obj, (xr.DataArray, xr.Dataset)):
-        if obj[dim].dtype == "O":
-            obj = obj[dim].where(obj[dim].notnull(), drop=True)[0].item()
-        elif "datetime64" in obj[dim].dtype.name:
-            return "default"
+        return obj[dim].dt.calendar
     elif isinstance(obj, xr.CFTimeIndex):
         obj = obj.values[0]
     else:
         obj = np.take(obj, 0)
         # Take zeroth element, overcome cases when arrays or lists are passed.
     if isinstance(obj, pydt.datetime):  # Also covers pandas Timestamp
-        return "default"
+        return "standard"
     if isinstance(obj, cftime.datetime):
         if obj.calendar == "gregorian":
             return "standard"
@@ -255,7 +231,7 @@ def _convert_doy_date(doy: int, year: int, src, tgt):
 
 
 def convert_doy(
-    source: xr.DataArray,
+    source: xr.DataArray | xr.Dataset,
     target_cal: str,
     source_cal: str | None = None,
     align_on: str = "year",
@@ -266,8 +242,9 @@ def convert_doy(
 
     Parameters
     ----------
-    source : xr.DataArray
+    source : xr.DataArray or xr.Dataset
       Day of year data (range [1, 366], max depending on the calendar).
+      If a Dataset, the function is mapped to each variables with attribute `is_day_of_year == 1`.
     target_cal : str
       Name of the calendar to convert to.
     source_cal : str, optional
@@ -282,6 +259,22 @@ def convert_doy(
     dim : str
       Name of the temporal dimension.
     """
+    if isinstance(source, xr.Dataset):
+        return source.map(
+            lambda da: (
+                da
+                if da.attrs.get("is_dayofyear") != 1
+                else convert_doy(
+                    da,
+                    target_cal,
+                    source_cal=source_cal,
+                    align_on=align_on,
+                    missing=missing,
+                    dim=dim,
+                )
+            )
+        )
+
     source_cal = source_cal or source.attrs.get("calendar", get_calendar(source[dim]))
     is_calyear = xr.infer_freq(source[dim]) in ("YS-JAN", "Y-DEC", "YE-DEC")
 
@@ -295,7 +288,7 @@ def convert_doy(
             max_doy_src = max_doy[source_cal]
         else:
             max_doy_src = xr.apply_ufunc(
-                days_in_year,
+                xr.coding.calendar_ops._days_in_year,
                 year_of_the_doy,
                 vectorize=True,
                 dask="parallelized",
@@ -305,7 +298,7 @@ def convert_doy(
             max_doy_tgt = max_doy[target_cal]
         else:
             max_doy_tgt = xr.apply_ufunc(
-                days_in_year,
+                xr.coding.calendar_ops._days_in_year,
                 year_of_the_doy,
                 vectorize=True,
                 dask="parallelized",
@@ -337,196 +330,29 @@ def convert_calendar(
     doy: bool | str = False,
     dim: str = "time",
 ) -> DataType:
-    """Convert a DataArray/Dataset to another calendar using the specified method.
-
-    By default, only converts the individual timestamps, does not modify any data except in dropping invalid/surplus dates or inserting missing dates.
-
-    If the source and target calendars are either no_leap, all_leap or a standard type, only the type of the time array is modified.
-    When converting to a leap year from a non-leap year, the 29th of February is removed from the array.
-    In the other direction and if `target` is a string, the 29th of February will be missing in the output,
-    unless `missing` is specified, in which case that value is inserted.
-
-    For conversions involving `360_day` calendars, see Notes.
-
-    This method is safe to use with sub-daily data as it doesn't touch the time part of the timestamps.
-
-    Parameters
-    ----------
-    source : xr.DataArray or xr.Dataset
-        Input array/dataset with a time coordinate of a valid dtype (datetime64 or a cftime.datetime).
-    target : xr.DataArray or str
-        Either a calendar name or the 1D time coordinate to convert to.
-        If an array is provided, the output will be reindexed using it and in that case, days in `target`
-        that are missing in the converted `source` are filled by `missing` (which defaults to NaN).
-    align_on : {None, 'date', 'year', 'random'}
-        Must be specified when either source or target is a `360_day` calendar, ignored otherwise. See Notes.
-    missing : Any, optional
-        A value to use for filling in dates in the target that were missing in the source.
-        If `target` is a string, default (None) is not to fill values. If it is an array, default is to fill with NaN.
-    doy: bool or {'year', 'date'}
-        If not False, variables flagged as "dayofyear" (with a `is_dayofyear==1` attribute) are converted to the new calendar too.
-        Can be a string, which will be passed as the `align_on` argument of :py:func:`convert_doy`.
-        If True, `year` is passed.
-    dim : str
-        Name of the time coordinate.
-
-    Returns
-    -------
-    xr.DataArray or xr.Dataset
-        Copy of source with the time coordinate converted to the target calendar.
-        If `target` is given as an array, the output is reindexed to it, with fill value `missing`.
-        If `target` was a string and `missing` was None (default), invalid dates in the new calendar are dropped,
-        but missing dates are not inserted.
-        If `target` was a string and `missing` was given, then start, end and frequency of the new time axis are
-        inferred and the output is reindexed to that a new array.
+    """Deprecated : use :py:meth:`xarray.Dataset.convert_calendar` or :py:meth:`xarray.DataArray.convert_calendar`
+    or :py:func:`xarray.coding.calendar_ops.convert_calendar` instead. Passing use_cftime=False instead of calendar='default'.
 
-    Notes
-    -----
-    If one of the source or target calendars is `360_day`, `align_on` must be specified and two options are offered.
-
-    "year"
-      The dates are translated according to their rank in the year (dayofyear), ignoring their original month and day information,
-      meaning that the missing/surplus days are added/removed at regular intervals.
-
-      From a `360_day` to a standard calendar, the output will be missing the following dates (day of year in parentheses):
-        To a leap year:
-          January 31st (31), March 31st (91), June 1st (153), July 31st (213), September 31st (275) and November 30th (335).
-        To a non-leap year:
-          February 6th (36), April 19th (109), July 2nd (183), September 12th (255), November 25th (329).
-
-      From standard calendar to a '360_day', the following dates in the source array will be dropped:
-        From a leap year:
-          January 31st (31), April 1st (92), June 1st (153), August 1st (214), September 31st (275), December 1st (336)
-        From a non-leap year:
-          February 6th (37), April 20th (110), July 2nd (183), September 13th (256), November 25th (329)
-
-      This option is best used on daily and subdaily data.
-
-    "date"
-      The month/day information is conserved and invalid dates are dropped from the output. This means that when
-      converting from a `360_day` to a standard calendar, all 31st (Jan, March, May, July, August, October and December)
-      will be missing as there is no equivalent dates in the `360_day` and the 29th (on non-leap years) and 30th of
-      February will be dropped as there are no equivalent dates in a standard calendar.
-
-      This option is best used with data on a frequency coarser than daily.
-
-    "random"
-      Similar to "year", each day of year of the source is mapped to another day of year of the target. However, instead
-      of having always the same missing days according the source and target years, here 5 days are chosen randomly, one
-      for each fifth of the year. However, February 29th is always missing when converting to a leap year, or its value
-      is dropped when converting from a leap year. This is similar to method used in the
-      :cite:t:`pierce_statistical_2014` dataset.
-
-      This option is best used on daily data.
-
-    References
-    ----------
-    :cite:cts:`pierce_statistical_2014`
-
-    Examples
-    --------
-    This method does not try to fill the missing dates other than with a constant value, passed with `missing`.
-    In order to fill the missing dates with interpolation, one can simply use xarray's method:
-
-    >>> tas_nl = convert_calendar(tas, "noleap")  # For the example
-    >>> with_missing = convert_calendar(tas_nl, "standard", missing=np.NaN)
-    >>> out = with_missing.interpolate_na("time", method="linear")
-
-    Here, if Nans existed in the source data, they will be interpolated too. If that is,
-    for some reason, not wanted, the workaround is to do:
-
-    >>> mask = convert_calendar(tas_nl, "standard").notnull()
-    >>> out2 = out.where(mask)
+    Convert a DataArray/Dataset to another calendar using the specified method.
     """
-    cal_src = get_calendar(source, dim=dim)
-
-    if isinstance(target, str):
-        cal_tgt = target
-    else:
-        cal_tgt = get_calendar(target, dim=dim)
-
-    if cal_src == cal_tgt:
-        return source
-
-    if (cal_src == "360_day" or cal_tgt == "360_day") and align_on not in [
-        "year",
-        "date",
-        "random",
-    ]:
-        raise ValueError(
-            "Argument `align_on` must be specified with either 'date', 'year' or "
-            "'random' when converting to or from a '360_day' calendar."
-        )
-    if cal_src != "360_day" and cal_tgt != "360_day":
-        align_on = None
-
-    if doy:
-        doy_align_on = "year" if doy is True else doy
-        if isinstance(source, xr.DataArray) and source.attrs.get("is_dayofyear") == 1:
-            out = convert_doy(source, cal_tgt, align_on=doy_align_on)
-        else:
-            out = source.map(
-                lambda da: (
-                    da
-                    if da.attrs.get("is_dayofyear") != 1
-                    else convert_doy(da, cal_tgt, align_on=doy_align_on)
-                )
-            )
-    else:
-        out = source.copy()
-
-    # TODO Maybe the 5-6 days to remove could be given by the user?
-    if align_on in ["year", "random"]:
-        if align_on == "year":
-            new_doy = source.time.groupby(f"{dim}.year").map(
-                yearly_interpolated_doy,
-                source_calendar=cal_src,
-                target_calendar=cal_tgt,
-            )
-        else:  # align_on == "random"
-            new_doy = source.time.groupby(f"{dim}.year").map(
-                yearly_random_doy,
-                rng=np.random.default_rng(),
-                source_calendar=cal_src,
-                target_calendar=cal_tgt,
-            )
-
-        # Convert the source datetimes, but override the doy with our new doys
-        out[dim] = xr.DataArray(
-            [
-                _convert_datetime(datetime, new_doy=doy, calendar=cal_tgt)
-                for datetime, doy in zip(source[dim].indexes[dim], new_doy)
-            ],
-            dims=(dim,),
-            name=dim,
-        )
-        # Remove NaN that where put on invalid dates in target calendar
-        out = out.where(out[dim].notnull(), drop=True)
-        # Remove duplicate timestamps, happens when reducing the number of days
-        out = out.isel({dim: np.unique(out[dim], return_index=True)[1]})
-    else:
-        time_idx = source[dim].indexes[dim]
-        out[dim] = xr.DataArray(
-            [_convert_datetime(time, calendar=cal_tgt) for time in time_idx],
-            dims=(dim,),
-            name=dim,
-        )
-        # Remove NaN that where put on invalid dates in target calendar
-        out = out.where(out[dim].notnull(), drop=True)
-
-    if isinstance(target, str) and missing is not None:
-        target = date_range_like(source[dim], cal_tgt)
-
     if isinstance(target, xr.DataArray):
-        out = out.reindex(
-            {dim: target}, fill_value=missing if missing is not None else np.nan
+        raise NotImplementedError(
+            "In `xclim` v0.50.0, `convert_calendar` is a direct copy of `xarray.coding.calendar_ops.convert_calendar`. "
+            "To retrieve the previous behaviour with target as a DataArray, convert the source first then reindex to the target."
         )
-
-    # Copy attrs but change remove `calendar` is still present.
-    out[dim].attrs.update(source[dim].attrs)
-    out[dim].attrs.pop("calendar", None)
-
-    return out
+    if doy is not False:
+        raise NotImplementedError(
+            "In `xclim` v0.50.0, `convert_calendar` is a direct copy of `xarray.coding.calendar_ops.convert_calendar`. "
+            "To retrieve the previous behaviour of doy=True, do convert_doy(obj, target_cal).convert_cal(target_cal)."
+        )
+    target, usecf = _get_usecf_and_warn(
+        target,
+        "convert_calendar",
+        "xarray.coding.calendar_ops.convert_calendar or obj.convert_calendar",
+    )
+    return xr.coding.calendar_ops.convert_calendar(
+        source, target, dim=dim, align_on=align_on, missing=missing
+    )
 
 
 def interp_calendar(
@@ -534,38 +360,14 @@ def interp_calendar(
     target: xr.DataArray,
     dim: str = "time",
 ) -> xr.DataArray | xr.Dataset:
-    """Interpolates a DataArray/Dataset to another calendar based on decimal year measure.
+    """Deprecated : use :py:func:`xarray.coding.calendar_ops.interp_calendar` instead.
 
-    Each timestamp in source and target are first converted to their decimal year equivalent
-    then source is interpolated on the target coordinate. The decimal year is the number of
-    years since 0001-01-01 AD.
-    Ex: '2000-03-01 12:00' is 2000.1653 in a standard calendar or 2000.16301 in a 'noleap' calendar.
-
-    This method should be used with daily data or coarser. Sub-daily result will have a modified day cycle.
-
-    Parameters
-    ----------
-    source : xr.DataArray or xr.Dataset
-        The source data to interpolate, must have a time coordinate of a valid dtype (np.datetime64 or cftime objects)
-    target : xr.DataArray
-        The target time coordinate of a valid dtype (np.datetime64 or cftime objects)
-    dim : str
-        The time coordinate name.
-
-    Return
-    ------
-    xr.DataArray or xr.Dataset
-        The source interpolated on the decimal years of target,
+    Interpolates a DataArray/Dataset to another calendar based on decimal year measure.
     """
-    cal_src = get_calendar(source, dim=dim)
-    cal_tgt = get_calendar(target, dim=dim)
-
-    out = source.copy()
-    out[dim] = datetime_to_decimal_year(source[dim], calendar=cal_src).drop_vars(dim)
-    target_idx = datetime_to_decimal_year(target, calendar=cal_tgt)
-    out = out.interp(time=target_idx)
-    out[dim] = target
-    return out
+    _, _ = _get_usecf_and_warn(
+        "standard", "interp_calendar", "xarray.coding.calendar_ops.interp_calendar"
+    )
+    return xr.coding.calendar_ops.interp_calendar(source, target, dim=dim)
 
 
 def ensure_cftime_array(time: Sequence) -> np.ndarray | Sequence[cftime.datetime]:
@@ -602,37 +404,18 @@ def ensure_cftime_array(time: Sequence) -> np.ndarray | Sequence[cftime.datetime
 
 
 def datetime_to_decimal_year(times: xr.DataArray, calendar: str = "") -> xr.DataArray:
-    """Convert a datetime xr.DataArray to decimal years according to its calendar or the given one.
-
-    Decimal years are the number of years since 0001-01-01 00:00:00 AD.
-    Ex: '2000-03-01 12:00' is 2000.1653 in a standard calendar, 2000.16301 in a "noleap" or 2000.16806 in a "360_day".
-
-    Parameters
-    ----------
-    times : xr.DataArray
-    calendar : str
+    """Deprecated : use :py:func:`xarray.coding.calendar_ops_datetime_to_decimal_year` instead.
 
-    Returns
-    -------
-    xr.DataArray
+    Convert a datetime xr.DataArray to decimal years according to its calendar or the given one.
     """
-    calendar = calendar or get_calendar(times)
-    if calendar == "default":
-        calendar = "standard"
-
-    def _make_index(time) -> xr.DataArray:
-        year = int(time.dt.year[0])
-        doys = cftime.date2num(
-            ensure_cftime_array(time), f"days since {year:04d}-01-01", calendar=calendar
-        )
-        return xr.DataArray(
-            year + doys / days_in_year(year, calendar),
-            dims=time.dims,
-            coords=time.coords,
-            name="time",
-        )
-
-    return times.groupby("time.year").map(_make_index)
+    _, _ = _get_usecf_and_warn(
+        "standard",
+        "datetime_to_decimal_year",
+        "xarray.coding.calendar_ops._datetime_to_decimal_year",
+    )
+    return xr.coding.calendar_ops._datetime_to_decimal_year(
+        times, dim="time", calendar=calendar
+    )
 
 
 @update_xclim_history
@@ -1231,7 +1014,10 @@ def _doy_days_since_doys(
     base_doy = base.dt.dayofyear
 
     doy_max = xr.apply_ufunc(
-        days_in_year, base.dt.year, vectorize=True, kwargs={"calendar": calendar}
+        xr.coding.calendar_ops._days_in_year,
+        base.dt.year,
+        vectorize=True,
+        kwargs={"calendar": calendar},
     )
 
     if start is not None:
@@ -1297,7 +1083,7 @@ def doy_to_days_since(
     """
     base_calendar = get_calendar(da)
     calendar = calendar or da.attrs.get("calendar", base_calendar)
-    dac = convert_calendar(da, calendar)
+    dac = da.convert_calendar(calendar)
 
     base_doy, start_doy, doy_max = _doy_days_since_doys(dac.time, start)
 
@@ -1317,7 +1103,7 @@ def doy_to_days_since(
 
     out.attrs.pop("is_dayofyear", None)
     out.attrs.update(calendar=calendar)
-    return convert_calendar(out, base_calendar).rename(da.name)
+    return out.convert_calendar(base_calendar).rename(da.name)
 
 
 def days_since_to_doy(
@@ -1366,7 +1152,7 @@ def days_since_to_doy(
     base_calendar = get_calendar(da)
     calendar = calendar or da.attrs.get("calendar", base_calendar)
 
-    dac = convert_calendar(da, calendar)
+    dac = da.convert_calendar(calendar)
 
     _, start_doy, doy_max = _doy_days_since_doys(dac.time, start)
 
@@ -1380,116 +1166,20 @@ def days_since_to_doy(
         {k: v for k, v in da.attrs.items() if k not in ["units", "calendar"]}
     )
     out.attrs.update(calendar=calendar, is_dayofyear=1)
-    return convert_calendar(out, base_calendar).rename(da.name)
+    return out.convert_calendar(base_calendar).rename(da.name)
 
 
 def date_range_like(source: xr.DataArray, calendar: str) -> xr.DataArray:
-    """Generate a datetime array with the same frequency, start and end as another one, but in a different calendar.
-
-    Parameters
-    ----------
-    source : xr.DataArray
-        1D datetime coordinate DataArray
-    calendar : str
-        New calendar name.
-
-    Raises
-    ------
-    ValueError
-        If the source's frequency was not found.
+    """Deprecated : use :py:func:`xarray.date_range_like` instead. Passing use_cftime=False instead of calendar='default'.
 
-    Returns
-    -------
-    xr.DataArray
-        1D datetime coordinate with the same start, end and frequency as the source, but in the new calendar.
-            The start date is assumed to exist in the target calendar.
-            If the end date doesn't exist, the code tries 1 and 2 calendar days before.
-            Exception when the source is in 360_day and the end of the range is the 30th of a 31-days month,
-            then the 31st is appended to the range.
+    Generate a datetime array with the same frequency, start and end as another one, but in a different calendar.
     """
-    freq = xr.infer_freq(source)
-    if freq is None:
-        raise ValueError(
-            "`date_range_like` was unable to generate a range as the source frequency was not inferrable."
-        )
-
-    src_cal = get_calendar(source)
-    if src_cal == calendar:
-        return source
-
-    index = source.indexes[source.dims[0]]
-    end_src = index[-1]
-    end = _convert_datetime(end_src, calendar=calendar)
-    if end is np.nan:  # Day is invalid, happens at the end of months.
-        end = _convert_datetime(end_src.replace(day=end_src.day - 1), calendar=calendar)
-        if end is np.nan:  # Still invalid : 360_day to non-leap february.
-            end = _convert_datetime(
-                end_src.replace(day=end_src.day - 2), calendar=calendar
-            )
-    if src_cal == "360_day" and end_src.day == 30 and end.daysinmonth == 31:
-        # For the specific case of daily data from 360_day source, the last day is expected to be "missing"
-        end = end.replace(day=31)
-
-    return xr.DataArray(
-        date_range(
-            _convert_datetime(index[0], calendar=calendar),
-            end,
-            freq=freq,
-            calendar=calendar,
-        ),
-        dims=source.dims,
-        name=source.dims[0],
+    calendar, usecf = _get_usecf_and_warn(
+        calendar, "date_range_like", "xarray.date_range_like"
+    )
+    return xr.coding.calendar_ops.date_range_like(
+        source=source, calendar=calendar, use_cftime=usecf
     )
-
-
-def _convert_datetime(
-    datetime: pydt.datetime | cftime.datetime,
-    new_doy: float | int | None = None,
-    calendar: str = "default",
-) -> cftime.datetime | pydt.datetime | float:
-    """Convert a datetime object to another calendar.
-
-    Nanosecond information are lost as cftime.datetime doesn't support them.
-
-    Parameters
-    ----------
-    datetime : datetime.datetime or cftime.datetime
-        A datetime object to convert.
-    new_doy : float or int, optional
-        Allows for redefining the day of year (thus ignoring month and day information from the source datetime).
-        -1 is understood as a nan.
-    calendar : str
-        The target calendar
-
-    Returns
-    -------
-    Union[cftime.datetime, datetime.datetime, np.nan]
-        A datetime object of the target calendar with the same year, month, day and time as the source
-        (month and day according to `new_doy` if given).
-        If the month and day doesn't exist in the target calendar, returns np.nan. (Ex. 02-29 in "noleap")
-    """
-    if new_doy in [np.nan, -1]:
-        return np.nan
-    if new_doy is not None:
-        new_date = cftime.num2date(
-            new_doy - 1,
-            f"days since {datetime.year}-01-01",
-            calendar=calendar if calendar != "default" else "standard",
-        )
-    else:
-        new_date = datetime
-    try:
-        return datetime_classes[calendar](
-            datetime.year,
-            new_date.month,
-            new_date.day,
-            datetime.hour,
-            datetime.minute,
-            datetime.second,
-            datetime.microsecond,
-        )
-    except ValueError:
-        return np.nan
 
 
 def select_time(
@@ -1597,7 +1287,7 @@ def _get_doys(_start, _end, _inclusive):
         if calendar not in uniform_calendars:
             # For non-uniform calendars, we can't simply convert dates to doys
             # conversion to all_leap is safe for all non-uniform calendar as it doesn't remove any date.
-            time = convert_calendar(time, "all_leap")
+            time = time.convert_calendar("all_leap")
             # values of time are the _old_ calendar
             # and the new calendar is in the coordinate
             calendar = "all_leap"
diff --git a/xclim/core/missing.py b/xclim/core/missing.py
index a7973f914..9a73bc1a9 100644
--- a/xclim/core/missing.py
+++ b/xclim/core/missing.py
@@ -28,13 +28,7 @@
 import numpy as np
 import xarray as xr
 
-from .calendar import (
-    date_range,
-    get_calendar,
-    is_offset_divisor,
-    parse_offset,
-    select_time,
-)
+from .calendar import get_calendar, is_offset_divisor, parse_offset, select_time
 from .options import (
     CHECK_MISSING,
     MISSING_METHODS,
@@ -157,11 +151,12 @@ def prepare(self, da, freq, src_timestep, **indexer):
         offset = parse_offset(src_timestep)
         if indexer or offset[1] in "YAQM":
             # Create a full synthetic time series and compare the number of days with the original series.
-            t = date_range(
+            t = xr.date_range(
                 start_time[0],
                 end_time[-1],
                 freq=src_timestep,
                 calendar=get_calendar(da),
+                use_cftime=(start_time.dtype == "O"),
             )
 
             sda = xr.DataArray(data=np.ones(len(t)), coords={"time": t}, dims=("time",))
diff --git a/xclim/core/units.py b/xclim/core/units.py
index 7f432eab7..f06e00040 100644
--- a/xclim/core/units.py
+++ b/xclim/core/units.py
@@ -22,7 +22,7 @@
 from boltons.funcutils import wraps
 from yaml import safe_load
 
-from .calendar import date_range, get_calendar, parse_offset
+from .calendar import get_calendar, parse_offset
 from .options import datacheck
 from .utils import InputKind, Quantified, ValidationError, infer_kind_from_parameter
 
@@ -598,8 +598,12 @@ def _rate_and_amount_converter(
                 label = "upper"
             # We generate "time" with an extra element, so we do not need to repeat the last element below.
             time = xr.DataArray(
-                date_range(
-                    start, periods=len(time) + 1, freq=freq, calendar=get_calendar(time)
+                xr.date_range(
+                    start,
+                    periods=len(time) + 1,
+                    freq=freq,
+                    calendar=get_calendar(time),
+                    use_cftime=(time.dtype == "O"),
                 ),
                 dims=(dim,),
                 name=dim,
diff --git a/xclim/ensembles/_base.py b/xclim/ensembles/_base.py
index 237368678..31988d559 100644
--- a/xclim/ensembles/_base.py
+++ b/xclim/ensembles/_base.py
@@ -13,7 +13,7 @@
 import numpy as np
 import xarray as xr
 
-from xclim.core.calendar import common_calendar, convert_calendar, get_calendar
+from xclim.core.calendar import common_calendar, get_calendar
 from xclim.core.formatting import update_history
 from xclim.core.utils import calc_perc
 
@@ -464,4 +464,4 @@ def _ens_align_datasets(
     if calendar is None:
         calendar = common_calendar(calendars, join="outer")
     cal_kwargs.setdefault("align_on", "date")
-    return [convert_calendar(ds, calendar, **cal_kwargs) for ds in ds_all]
+    return [ds.convert_calendar(calendar, **cal_kwargs) for ds in ds_all]
diff --git a/xclim/indices/generic.py b/xclim/indices/generic.py
index 76d987024..a0e919f62 100644
--- a/xclim/indices/generic.py
+++ b/xclim/indices/generic.py
@@ -17,12 +17,7 @@
 import xarray as xr
 from xarray.coding.cftime_offsets import _MONTH_ABBREVIATIONS  # noqa
 
-from xclim.core.calendar import (
-    convert_calendar,
-    doy_to_days_since,
-    get_calendar,
-    select_time,
-)
+from xclim.core.calendar import doy_to_days_since, get_calendar, select_time
 from xclim.core.units import (
     convert_units_to,
     declare_relative_units,
@@ -825,11 +820,11 @@ def _get_days(_bound, _group, _base_time):
     cal = get_calendar(data, dim="time")
 
     if not isinstance(start, str):
-        start = convert_calendar(start, cal)
+        start = start.convert_calendar(cal)
         start.attrs["calendar"] = cal
         start = doy_to_days_since(start)
     if not isinstance(end, str):
-        end = convert_calendar(end, cal)
+        end = end.convert_calendar(cal)
         end.attrs["calendar"] = cal
         end = doy_to_days_since(end)
 
diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py
index 454819052..057c39305 100644
--- a/xclim/indices/helpers.py
+++ b/xclim/indices/helpers.py
@@ -16,12 +16,11 @@
 import numba as nb
 import numpy as np
 import xarray as xr
-
-from xclim.core.calendar import (
-    datetime_to_decimal_year,
-    ensure_cftime_array,
-    get_calendar,
+from xarray.coding.calendar_ops import (
+    _datetime_to_decimal_year as datetime_to_decimal_year,
 )
+
+from xclim.core.calendar import ensure_cftime_array, get_calendar
 from xclim.core.units import convert_units_to
 from xclim.core.utils import Quantified, _chunk_like
 
diff --git a/xclim/sdba/base.py b/xclim/sdba/base.py
index af7652f70..6d2b6b593 100644
--- a/xclim/sdba/base.py
+++ b/xclim/sdba/base.py
@@ -15,7 +15,7 @@
 import xarray as xr
 from boltons.funcutils import wraps
 
-from xclim.core.calendar import days_in_year, get_calendar
+from xclim.core.calendar import get_calendar
 from xclim.core.options import OPTIONS, SDBA_ENCODE_CF
 from xclim.core.utils import uses_dask
 
@@ -197,7 +197,8 @@ def get_coordinate(self, ds: xr.Dataset | None = None) -> xr.DataArray:
             if ds is not None:
                 cal = get_calendar(ds, dim=self.dim)
                 mdoy = max(
-                    days_in_year(yr, cal) for yr in np.unique(ds[self.dim].dt.year)
+                    xr.coding.calendar_ops._days_in_year(yr, cal)
+                    for yr in np.unique(ds[self.dim].dt.year)
                 )
             else:
                 mdoy = 365
diff --git a/xclim/testing/helpers.py b/xclim/testing/helpers.py
index 94d87fa36..935750101 100644
--- a/xclim/testing/helpers.py
+++ b/xclim/testing/helpers.py
@@ -11,7 +11,7 @@
 import xarray as xr
 from dask.diagnostics import Callback
 
-from xclim.core import calendar
+from xclim.core.calendar import percentile_doy
 from xclim.core.utils import VARIABLES
 from xclim.indices import (
     longwave_upwelling_radiation_from_net_downwelling,
@@ -82,10 +82,10 @@ def generate_atmos(cache_dir: Path):
         branch=TESTDATA_BRANCH,
         engine="h5netcdf",
     ) as ds:
-        tn10 = calendar.percentile_doy(ds.tasmin, per=10)
-        t10 = calendar.percentile_doy(ds.tas, per=10)
-        t90 = calendar.percentile_doy(ds.tas, per=90)
-        tx90 = calendar.percentile_doy(ds.tasmax, per=90)
+        tn10 = percentile_doy(ds.tasmin, per=10)
+        t10 = percentile_doy(ds.tas, per=10)
+        t90 = percentile_doy(ds.tas, per=90)
+        tx90 = percentile_doy(ds.tasmax, per=90)
 
         rsus = shortwave_upwelling_radiation_from_net_downwelling(ds.rss, ds.rsds)
         rlus = longwave_upwelling_radiation_from_net_downwelling(ds.rls, ds.rlds)