def add_doy(ts, tdim="time"):
    """Add coordinate 'doy' day of the year for a 366 days year.

    Parameters
    ----------
    ts: xarray DataArray
        Timeseries array
    tdim: str, optional
        Name of time dimension (default='time')

    Returns
    -------
    ts: xarray DataArray
        Timeseries array with extra 'doy' coordinate
    """

    # get original dayofyear
    # create filters: from 1st of March onwards and non leap years
    # add extra day if not leap year and march or later
    t = ts[tdim]
    doy_original = t.dt.dayofyear
    march_or_later = t.dt.month >= 3
    not_leap_year = ~t.dt.is_leap_year
    doy = doy_original + (not_leap_year & march_or_later)
    # rechunk and return new doy as coordinate of the "ts" input variable
    ts.coords['doy'] = doy.chunk({tdim: -1})
    return ts