diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index a10a8c8851f..00251ba883b 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -57,6 +57,8 @@ Bug fixes
   Haacker <https://github.com/j-haacker>`_.
 - Fix ``isel`` for multi-coordinate Xarray indexes (:issue:`10063`, :pull:`10066`).
   By `Benoit Bovy <https://github.com/benbovy>`_.
+- Improve handling of dtype and NaT when encoding/decoding masked and packaged datetimes and timedeltas (:issue:`8957`, :pull:`10050`).
+  By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_.
 
 
 Documentation
diff --git a/xarray/coding/variables.py b/xarray/coding/variables.py
index 77e2f0602ce..a2534df0ec7 100644
--- a/xarray/coding/variables.py
+++ b/xarray/coding/variables.py
@@ -502,10 +502,8 @@ def decode(self, variable: Variable, name: T_Name = None):
                     (is_time_like == "datetime" and self.decode_times)
                     or (is_time_like == "timedelta" and self.decode_timedelta)
                 ) and data.dtype.kind in "iu":
-                    dtype, decoded_fill_value = (
-                        np.int64,
-                        np.iinfo(np.int64).min,
-                    )  # np.dtype(f"{is_time_like}64[s]")
+                    dtype = np.int64
+                    decoded_fill_value = np.iinfo(np.int64).min
                 else:
                     dtype, decoded_fill_value = dtypes.maybe_promote(data.dtype)
 
@@ -624,17 +622,17 @@ def decode(self, variable: Variable, name: T_Name = None) -> Variable:
                 scale_factor = np.asarray(scale_factor).item()
             if np.ndim(add_offset) > 0:
                 add_offset = np.asarray(add_offset).item()
-            # if we have a _FillValue/masked_value we already have the wanted
+            # if we have a _FillValue/masked_value in encoding we already have the wanted
             # floating point dtype here (via CFMaskCoder), so no check is necessary
             # only check in other cases and for time-like
             dtype = data.dtype
             is_time_like = _is_time_like(attrs.get("units"))
-            if (is_time_like == "datetime" and self.decode_times) or (
-                is_time_like == "timedelta" and self.decode_timedelta
+            if (
+                ("_FillValue" not in encoding and "missing_value" not in encoding)
+                or (is_time_like == "datetime" and self.decode_times)
+                or (is_time_like == "timedelta" and self.decode_timedelta)
             ):
                 dtype = _choose_float_dtype(dtype, encoding)
-            if "_FillValue" not in encoding and "missing_value" not in encoding:
-                dtype = _choose_float_dtype(dtype, encoding)
 
             transform = partial(
                 _scale_offset_decoding,