diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index cb9b63479..10cd6dbc4 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,6 +1,6 @@
 repos:
 -   repo: https://github.com/psf/black
-    rev: 25.1.0
+    rev: 26.1.0
     hooks:
     - id: black
       language_version: python3
diff --git a/examples/anvil_nowcast.py b/examples/anvil_nowcast.py
index f7fdd4033..898f56fb3 100644
--- a/examples/anvil_nowcast.py
+++ b/examples/anvil_nowcast.py
@@ -9,6 +9,7 @@
 
 Load the libraries.
 """
+
 from datetime import datetime, timedelta
 import warnings
 
diff --git a/examples/ens_kalman_filter_blended_forecast.py b/examples/ens_kalman_filter_blended_forecast.py
index b0a9cdf39..e139b1b0c 100644
--- a/examples/ens_kalman_filter_blended_forecast.py
+++ b/examples/ens_kalman_filter_blended_forecast.py
@@ -30,7 +30,6 @@
 from pysteps.visualization import plot_precip_field
 import pysteps_nwp_importers
 
-
 ################################################################################
 # Read the radar images and the NWP forecast
 # ------------------------------------------
diff --git a/examples/optical_flow_methods_convergence.py b/examples/optical_flow_methods_convergence.py
index 569c6709f..c5948ac4c 100644
--- a/examples/optical_flow_methods_convergence.py
+++ b/examples/optical_flow_methods_convergence.py
@@ -17,6 +17,7 @@
 
 Let's first load the libraries that we will use.
 """
+
 from datetime import datetime
 import time
 
diff --git a/examples/plot_custom_precipitation_range.py b/examples/plot_custom_precipitation_range.py
index 4f1a98df1..6e590608d 100644
--- a/examples/plot_custom_precipitation_range.py
+++ b/examples/plot_custom_precipitation_range.py
@@ -18,7 +18,6 @@
 from pysteps.visualization import plot_precip_field
 from pysteps.datasets import download_pysteps_data, create_default_pystepsrc
 
-
 ###############################################################################
 # Download the data if it is not available
 # ----------------------------------------
diff --git a/examples/plot_ensemble_verification.py b/examples/plot_ensemble_verification.py
index 5ebf2aaff..4e7d76d04 100644
--- a/examples/plot_ensemble_verification.py
+++ b/examples/plot_ensemble_verification.py
@@ -18,7 +18,6 @@
 from pysteps.utils import conversion, dimension, transformation
 from pysteps.visualization import plot_precip_field
 
-
 ###############################################################################
 # Read precipitation field
 # ------------------------
diff --git a/examples/plot_linear_blending.py b/examples/plot_linear_blending.py
index 91cc267c0..4fe5cae19 100644
--- a/examples/plot_linear_blending.py
+++ b/examples/plot_linear_blending.py
@@ -19,7 +19,6 @@
 from pysteps.utils import conversion
 from pysteps.visualization import plot_precip_field
 
-
 ################################################################################
 # Read the radar images and the NWP forecast
 # ------------------------------------------
diff --git a/examples/steps_blended_forecast.py b/examples/steps_blended_forecast.py
index cd9562203..484dc4207 100644
--- a/examples/steps_blended_forecast.py
+++ b/examples/steps_blended_forecast.py
@@ -19,7 +19,6 @@
 from pysteps import io, rcparams, blending, nowcasts
 from pysteps.visualization import plot_precip_field
 
-
 ################################################################################
 # Read the radar images and the NWP forecast
 # ------------------------------------------
diff --git a/examples/thunderstorm_detection_and_tracking.py b/examples/thunderstorm_detection_and_tracking.py
index f44d28b90..597c199f1 100644
--- a/examples/thunderstorm_detection_and_tracking.py
+++ b/examples/thunderstorm_detection_and_tracking.py
@@ -23,6 +23,7 @@
 
 @author: feldmann-m
 """
+
 ################################################################################
 # Import all required functions
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/pysteps/blending/ens_kalman_filter_methods.py b/pysteps/blending/ens_kalman_filter_methods.py
index 5124a87e7..161b7242a 100644
--- a/pysteps/blending/ens_kalman_filter_methods.py
+++ b/pysteps/blending/ens_kalman_filter_methods.py
@@ -63,7 +63,6 @@
     forecast. Defaults to True.
 """
 
-
 import numpy as np
 
 from pysteps import utils
diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py
index d93af1db3..8f9bce2a5 100644
--- a/pysteps/blending/steps.py
+++ b/pysteps/blending/steps.py
@@ -43,6 +43,7 @@
     calculate_weights_spn
     blend_means_sigmas
 """
+
 import math
 import time
 from copy import copy, deepcopy
@@ -137,6 +138,14 @@ class StepsBlendingConfig:
       The calculation method of the blending weights. Options are the method
       by :cite:`BPS2006` and the covariance-based method by :cite:`SPN2013`.
       Defaults to bps.
+    timestep_start_full_nwp_weight: int, optional.
+      The timestep, which should be smaller than timesteps, at which a linear
+      transition takes place from the calculated weights to full NWP weight
+      (and zero extrapolation and noise weight) to ensure the blending
+      procedure becomes equal to the NWP forecast(s) at the last timestep
+      of the blending procedure. If not provided, the blending stick to the
+      theoretical weights provided by the chosen weights_method for a given
+      lead time and skill of each blending component.
     conditional: bool, optional
       If set to True, compute the statistics of the precipitation field
       conditionally by excluding pixels where the values are below the threshold
@@ -290,6 +299,7 @@ class StepsBlendingConfig:
     ar_order: int
     velocity_perturbation_method: str | None
     weights_method: str
+    timestep_start_full_nwp_weight: int
     conditional: bool
     probmatching_method: str | None
     mask_method: str | None
@@ -406,6 +416,8 @@ class StepsBlendingState:
     # Final outputs
     final_blended_forecast: np.ndarray | None = None
     final_blended_forecast_non_perturbed: np.ndarray | None = None
+    weights: np.ndarray | None = None
+    weights_model_only: np.ndarray | None = None
 
     # Timing and indexing
     time_prev_timestep: list[float] | None = None
@@ -644,7 +656,7 @@ def __blended_nowcast_main_loop(self):
             def worker(j):
                 worker_state = copy(self.__state)
                 self.__determine_NWP_skill_for_next_timestep(t, j, worker_state)
-                self.__determine_weights_per_component(worker_state)
+                self.__determine_weights_per_component(t, worker_state)
                 self.__regress_extrapolation_and_noise_cascades(j, worker_state, t)
                 self.__perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep(
                     t, j, worker_state
@@ -664,6 +676,7 @@ def worker(j):
                                 worker_state,
                             )
                         )
+
                     final_blended_forecast_all_members_one_timestep[j] = (
                         final_blended_forecast_single_member
                     )
@@ -889,6 +902,18 @@ def __check_inputs(self):
             if self.__config.mask_method == "incremental":
                 raise ValueError("mask_method='incremental' but timestep=None")
 
+        if self.__config.timestep_start_full_nwp_weight is not None:
+            if self.__config.timestep_start_full_nwp_weight < 0:
+                raise ValueError(
+                    "timestep_start_full_nwp_weight cannot be smaller than zero"
+                )
+
+        if self.__config.timestep_start_full_nwp_weight is not None:
+            if self.__config.timestep_start_full_nwp_weight >= self.__timesteps[-1]:
+                raise ValueError(
+                    "timestep_start_full_nwp_weight cannot be the same or larger than the total number of timesteps in this forecast"
+                )
+
     def __print_forecast_info(self):
         """
         Print information about the forecast setup, including inputs, methods, and parameters.
@@ -2051,29 +2076,48 @@ def __determine_NWP_skill_for_next_timestep(self, t, j, worker_state):
                 axis=0,
             )
 
-    def __determine_weights_per_component(self, worker_state):
+    def __determine_weights_per_component(self, t, worker_state):
         """
         Compute blending weights for each component based on the selected method
         ('bps' or 'spn'). Weights are determined for both full blending and
         model-only scenarios, accounting for correlations and covariance.
         """
+        start_smoothing_to_final_weights = False
+        if self.__config.timestep_start_full_nwp_weight is not None:
+            if t > self.__config.timestep_start_full_nwp_weight:
+                start_smoothing_to_final_weights = True
         # Weights following the bps method. These are needed for the velocity
         # weights prior to the advection step. If weights method spn is
         # selected, weights will be overwritten with those weights prior to
         # blending step.
         # weight = [(extr_field, n_model_fields, noise), n_cascade_levels, ...]
-        worker_state.weights = calculate_weights_bps(
-            worker_state.rho_final_blended_forecast
-        )
+        if not start_smoothing_to_final_weights:
+            worker_state.weights = calculate_weights_bps(
+                worker_state.rho_final_blended_forecast
+            )
+        else:
+            worker_state.weights = calculate_end_weights(
+                previous_weights=self.__state.weights,
+                timestep=t,
+                n_timesteps=self.__timesteps[-1],
+                start_full_nwp_weight=self.__config.timestep_start_full_nwp_weight,
+                model_only=False,
+            )
 
         # The model only weights
-        if self.__config.weights_method == "bps":
+        if (
+            self.__config.weights_method == "bps"
+            and not start_smoothing_to_final_weights
+        ):
             # Determine the weights of the components without the extrapolation
             # cascade, in case this is no data or outside the mask.
             worker_state.weights_model_only = calculate_weights_bps(
                 worker_state.rho_final_blended_forecast[1:, :]
             )
-        elif self.__config.weights_method == "spn":
+        elif (
+            self.__config.weights_method == "spn"
+            and not start_smoothing_to_final_weights
+        ):
             # Only the weights of the components without the extrapolation
             # cascade will be determined here. The full set of weights are
             # determined after the extrapolation step in this method.
@@ -2112,13 +2156,23 @@ def __determine_weights_per_component(self, worker_state):
             else:
                 # Same as correlation and noise is 1 - correlation
                 worker_state.weights_model_only = calculate_weights_bps(
-                    worker_state.rho_final_blended_forecast[1:, :]
+                    worker_state.rho_final_blended_forecast[1:, :],
                 )
+        elif start_smoothing_to_final_weights:
+            worker_state.weights_model_only = calculate_end_weights(
+                previous_weights=self.__state.weights_model_only,
+                timestep=t,
+                n_timesteps=self.__timesteps[-1],
+                start_full_nwp_weight=self.__config.timestep_start_full_nwp_weight,
+                model_only=True,
+            )
         else:
             raise ValueError(
                 "Unknown weights method %s: must be 'bps' or 'spn'"
                 % self.__config.weights_method
             )
+        self.__state.weights = worker_state.weights
+        self.__state.weights_model_only = worker_state.weights_model_only
 
     def __regress_extrapolation_and_noise_cascades(self, j, worker_state, t):
         """
@@ -2774,8 +2828,16 @@ def __blend_cascades(self, t_sub, j, worker_state):
             )
 
         # First determine the blending weights if method is spn. The
-        # weights for method bps have already been determined.
-        if self.__config.weights_method == "spn":
+        # weights for method bps have already been determined.'
+        start_smoothing_to_final_weights = False
+        if self.__config.timestep_start_full_nwp_weight is not None:
+            if t_sub >= self.__config.timestep_start_full_nwp_weight:
+                start_smoothing_to_final_weights = True
+
+        if (
+            self.__config.weights_method == "spn"
+            and not start_smoothing_to_final_weights
+        ):
             worker_state.weights = np.zeros(
                 (
                     cascade_stack_all_components.shape[0],
@@ -2800,6 +2862,8 @@ def __blend_cascades(self, t_sub, j, worker_state):
                     covariance=covariance_nwp_models,
                 )
 
+            self.__state.weights = worker_state.weights
+
         # Create weights_with_noise to ensure there is always a 3D weights field, even
         # if self.__config.nowcasting_method is "external_nowcast" and n_ens_members is 1.
         worker_state.weights_with_noise = worker_state.weights.copy()
@@ -3237,6 +3301,7 @@ def forecast(
     ar_order=2,
     vel_pert_method="bps",
     weights_method="bps",
+    timestep_start_full_nwp_weight=None,
     conditional=False,
     probmatching_method="cdf",
     mask_method="incremental",
@@ -3380,6 +3445,14 @@ def forecast(
       The calculation method of the blending weights. Options are the method
       by :cite:`BPS2006` and the covariance-based method by :cite:`SPN2013`.
       Defaults to bps.
+    timestep_start_full_nwp_weight: int, optional.
+      The timestep, which should be smaller than timesteps, at which a linear
+      transition takes place from the calculated weights to full (1.0) NWP weight
+      (and zero extrapolation and noise weight) to ensure the blending
+      procedure becomes equal to the NWP forecast(s) at the last timestep
+      of the blending procedure. If not provided, the blending stick to the
+      theoretical weights provided by the chosen weights_method for a given
+      lead time and skill of each blending component.
     conditional: bool, optional
       If set to True, compute the statistics of the precipitation field
       conditionally by excluding pixels where the values are below the threshold
@@ -3569,6 +3642,7 @@ def forecast(
         ar_order=ar_order,
         velocity_perturbation_method=vel_pert_method,
         weights_method=weights_method,
+        timestep_start_full_nwp_weight=timestep_start_full_nwp_weight,
         conditional=conditional,
         probmatching_method=probmatching_method,
         mask_method=mask_method,
@@ -3678,6 +3752,7 @@ def calculate_weights_bps(correlations):
         # total_ratios: [scale, ...] - the denominator of eq. 11 & 12 in BPS2006
         weights = correlations * np.sqrt(ratios / total_ratios)
         # weights: [component, scale, ...]
+
         # Calculate the weight of the noise component.
         # Original BPS2006 method in the following two lines (eq. 13)
         total_square_weights = np.sum(np.square(weights), axis=0)
@@ -3779,6 +3854,112 @@ def calculate_weights_spn(correlations, covariance):
     return weights
 
 
+# TODO: Where does this piece of code best fit: in utils or inside the class?
+def calculate_end_weights(
+    previous_weights, timestep, n_timesteps, start_full_nwp_weight, model_only=False
+):
+    """Calculate the linear transition from the previous weights to the final weights
+    (1.0 for NWP and 0.0 for the extrapolation and noise components). This method uses
+    the BPS weights determination method to determine the corresponding noise.
+
+    Parameters
+    ----------
+    previous_weights : array-like
+      The weights from the previous timestep. This weight will be used to ensure
+      a linear transition takes place from the last weights at the timestep of
+      start_full_nwp_weight and the final weights (1.0 for NWP and 0.0 for
+      the extrapolation and noise components).
+    timestep : int
+      The timestep or sub timestep for which the weight is calculated. Only
+      used when start_full_nwp_weight is not None.
+    n_timesteps: int
+      The total number of forecast timesteps in the forecast.
+    start_full_nwp_weight : int
+      The timestep, which should be smaller than timesteps, at which a linear
+      transition takes place from the calculated weights to full NWP weight
+      (and zero extrapolation and noise weight) to ensure the blending
+      procedure becomes equal to the NWP forecast(s) at the last timestep
+      of the blending procedure. If not provided, the blending stick to the
+      theoretical weights provided by the chosen weights_method for a given
+      lead time and skill of each blending component.
+    model_only : bool
+      If set to True, the weights will only be determined for the model and
+      noise components.
+
+    Returns
+    -------
+    weights : array-like
+      Array of shape [component+1, scale_level, ...]
+      containing the weights to be used in STEPS blending for
+      each original component plus an addtional noise component, scale level,
+      and optionally along [y, x] dimensions.
+
+    References
+    ----------
+    :cite:`BPS2006`
+
+    Notes
+    -----
+    The weights in the BPS method can sum op to more than 1.0.
+    """
+    weights = previous_weights[:-1, :].copy()
+    if not model_only:
+        if timestep > start_full_nwp_weight and timestep < n_timesteps:
+            weights[0, :] = weights[0, :] - (
+                (timestep - start_full_nwp_weight)
+                / (n_timesteps - start_full_nwp_weight)
+                * weights[0, :]
+            )
+            weights[1:, :] = (
+                1.0
+                / weights[1:, :].shape[0]
+                * (
+                    weights[1:, :]
+                    + (
+                        (timestep - start_full_nwp_weight)
+                        / (n_timesteps - start_full_nwp_weight)
+                        * (1.0 - weights[1:, :])
+                    )
+                )
+            )
+        elif timestep > start_full_nwp_weight and timestep == n_timesteps:
+            weights[0, :] = 0.0
+            # If one model or model member is provided to blend together,
+            # the weight equals 1.0, otherwise the sum of the weights
+            # equals 1.0.
+            weights[1:, :] = 1.0 / weights[1:, :].shape[0]
+
+    else:
+        if timestep > start_full_nwp_weight and timestep < n_timesteps:
+            weights = (
+                1.0
+                / weights.shape[0]
+                * (
+                    weights
+                    + (
+                        (timestep - start_full_nwp_weight)
+                        / (n_timesteps - start_full_nwp_weight)
+                        * (1.0 - weights)
+                    )
+                )
+            )
+        elif timestep > start_full_nwp_weight and timestep == n_timesteps:
+            weights[:] = 1.0 / weights.shape[0]
+
+    if weights.shape[0] > 1:
+        # Calculate the weight of the noise component.
+        # Original BPS2006 method in the following two lines (eq. 13)
+        total_square_weights = np.sum(np.square(weights), axis=0)
+        noise_weight = np.sqrt(1.0 - total_square_weights)
+        # Finally, add the noise_weights to the weights variable.
+        weights = np.concatenate((weights, noise_weight[None, ...]), axis=0)
+    else:
+        noise_weight = 1.0 - weights
+        weights = np.concatenate((weights, noise_weight), axis=0)
+
+    return weights
+
+
 # TODO: Where does this piece of code best fit: in utils or inside the class?
 def blend_means_sigmas(means, sigmas, weights):
     """Calculate the blended means and sigmas, the normalization parameters
diff --git a/pysteps/datasets.py b/pysteps/datasets.py
index 91abc62ed..173d5d5d8 100644
--- a/pysteps/datasets.py
+++ b/pysteps/datasets.py
@@ -14,6 +14,7 @@
     info
     load_dataset
 """
+
 import gzip
 import json
 import os
diff --git a/pysteps/decorators.py b/pysteps/decorators.py
index 44fbaebdb..23c85cc84 100644
--- a/pysteps/decorators.py
+++ b/pysteps/decorators.py
@@ -14,6 +14,7 @@
     prepare_interpolator
     memoize
 """
+
 import inspect
 import uuid
 import warnings
diff --git a/pysteps/downscaling/rainfarm.py b/pysteps/downscaling/rainfarm.py
index b40e359df..98ef1b8cd 100644
--- a/pysteps/downscaling/rainfarm.py
+++ b/pysteps/downscaling/rainfarm.py
@@ -169,7 +169,7 @@ def _compute_kernel_radius(ds_factor):
 def _make_tophat_kernel(ds_factor):
     """Compute 2d uniform (tophat) kernel"""
     radius = _compute_kernel_radius(ds_factor)
-    (mx, my) = np.mgrid[-radius : radius + 0.01, -radius : radius + 0.01]
+    mx, my = np.mgrid[-radius : radius + 0.01, -radius : radius + 0.01]
     tophat = ((mx**2 + my**2) <= radius**2).astype(float)
     return tophat / tophat.sum()
 
diff --git a/pysteps/io/interface.py b/pysteps/io/interface.py
index b1dc74855..917362532 100644
--- a/pysteps/io/interface.py
+++ b/pysteps/io/interface.py
@@ -12,6 +12,7 @@
 
     get_method
 """
+
 from importlib.metadata import entry_points
 
 from pysteps.decorators import postprocess_import
diff --git a/pysteps/motion/interface.py b/pysteps/motion/interface.py
index eac81aa46..52742389a 100644
--- a/pysteps/motion/interface.py
+++ b/pysteps/motion/interface.py
@@ -23,6 +23,7 @@
 
     get_method
 """
+
 import numpy as np
 
 from pysteps.motion.constant import constant
diff --git a/pysteps/postprocessing/interface.py b/pysteps/postprocessing/interface.py
index 130f7cab7..10a7dc2c6 100644
--- a/pysteps/postprocessing/interface.py
+++ b/pysteps/postprocessing/interface.py
@@ -16,6 +16,7 @@
 
     get_method
 """
+
 import importlib
 from importlib.metadata import entry_points
 
diff --git a/pysteps/tests/test_blending_clim.py b/pysteps/tests/test_blending_clim.py
index 0cb867a73..6508a959f 100644
--- a/pysteps/tests/test_blending_clim.py
+++ b/pysteps/tests/test_blending_clim.py
@@ -12,7 +12,6 @@
 
 from pysteps.blending.clim import save_skill, calc_clim_skill
 
-
 random.seed(12356)
 n_cascade_levels = 7
 model_names = ["alaro13", "arome13"]
diff --git a/pysteps/tests/test_blending_skill_scores.py b/pysteps/tests/test_blending_skill_scores.py
index 9403531d2..cec22be9a 100644
--- a/pysteps/tests/test_blending_skill_scores.py
+++ b/pysteps/tests/test_blending_skill_scores.py
@@ -11,7 +11,6 @@
     clim_regr_values,
 )
 
-
 # Set the climatological correlations values
 clim_cor_values_8lev = np.array(
     [0.848, 0.537, 0.237, 0.065, 0.020, 0.0044, 0.0052, 0.0040]
diff --git a/pysteps/tests/test_blending_steps.py b/pysteps/tests/test_blending_steps.py
index f914638f7..fe9e2863b 100644
--- a/pysteps/tests/test_blending_steps.py
+++ b/pysteps/tests/test_blending_steps.py
@@ -10,64 +10,69 @@
 
 # fmt:off
 steps_arg_values = [
-    (1, 3, 4, 8, 'steps', None, None, False, "spn", True, 4, False, False, 0, False, None, None),
-    (1, 3, 4, 8,'steps', "obs", None, False, "spn", True, 4, False, False, 0, False, None, None),
-    (1, 3, 4, 8,'steps', "incremental", None, False, "spn", True, 4, False, False, 0, False, None, None),
-    (1, 3, 4, 8,'steps', None, "mean", False, "spn", True, 4, False, False, 0, False, None, None),
-    (1, 3, 4, 8,'steps', None, "mean", False, "spn", True, 4, False, False, 0, True, None, None),
-    (1, 3, 4, 8,'steps', None, "cdf", False, "spn", True, 4, False, False, 0, False, None, None),
-    (1, [1, 2, 3], 4, 8,'steps', None, "cdf", False, "spn", True, 4, False, False, 0, False, None, None),
-    (1, 3, 4, 8,'steps', "incremental", "cdf", False, "spn", True, 4, False, False, 0, False, None, None),
-    (1, 3, 4, 6,'steps', "incremental", "cdf", False, "bps", True, 4, False, False, 0, False, None, None),
-    (1, 3, 4, 6,'steps', "incremental", "cdf", False, "bps", False, 4, False, False, 0, False, None, None),
-    (1, 3, 4, 6,'steps', "incremental", "cdf", False, "bps", False, 4, False, False, 0, True, None, None),
-    (1, 3, 4, 9,'steps', "incremental", "cdf", False, "spn", True, 4, False, False, 0, False, None, None),
-    (2, 3, 10, 8,'steps', "incremental", "cdf", False, "spn", True, 10, False, False, 0, False, None, None),
-    (5, 3, 5, 8,'steps', "incremental", "cdf", False, "spn", True, 5, False, False, 0, False, None, None),
-    (1, 10, 1, 8,'steps', "incremental", "cdf", False, "spn", True, 1, False, False, 0, False, None, None),
-    (2, 3, 2, 8,'steps', "incremental", "cdf", True, "spn", True, 2, False, False, 0, False, None, None),
-    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, False, False, 0, False, None, None),
-    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, False, False, 0, False, "bps", None),
-    #    Test the case where the radar image contains no rain.
-    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, True, False, 0, False, None, None),
-    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, True, False, 0, False, None, None),
-    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, True, False, 0, True, None, None),
-    #   Test the case where the NWP fields contain no rain.
-    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, False, True, 0, False, None, None),
-    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, False, True, 0, True, None, None),
+    (1, 3, 4, 8, 'steps', None, None, False, "spn", True, 4, False, False, 0, False, None, None, None),
+    (1, 3, 4, 8,'steps', "obs", None, False, "spn", True, 4, False, False, 0, False, None, None, None),
+    (1, 3, 4, 8,'steps', "incremental", None, False, "spn", True, 4, False, False, 0, False, None, None, None),
+    (1, 3, 4, 8,'steps', None, "mean", False, "spn", True, 4, False, False, 0, False, None, None, None),
+    (1, 3, 4, 8,'steps', None, "mean", False, "spn", True, 4, False, False, 0, True, None, None, None),
+    (1, 3, 4, 8,'steps', None, "cdf", False, "spn", True, 4, False, False, 0, False, None, None, None),
+    (1, [1, 2, 3], 4, 8,'steps', None, "cdf", False, "spn", True, 4, False, False, 0, False, None, None, None),
+    (1, 3, 4, 8,'steps', "incremental", "cdf", False, "spn", True, 4, False, False, 0, False, None, None, None),
+    (1, 3, 4, 6,'steps', "incremental", "cdf", False, "bps", True, 4, False, False, 0, False, None, None, None),
+    (1, 3, 4, 6,'steps', "incremental", "cdf", False, "bps", False, 4, False, False, 0, False, None, None, None),
+    (1, 3, 4, 6,'steps', "incremental", "cdf", False, "bps", False, 4, False, False, 0, True, None, None, None),
+    (1, 3, 4, 9,'steps', "incremental", "cdf", False, "spn", True, 4, False, False, 0, False, None, None, None),
+    (2, 3, 10, 8,'steps', "incremental", "cdf", False, "spn", True, 10, False, False, 0, False, None, None, None),
+    (5, 3, 5, 8,'steps', "incremental", "cdf", False, "spn", True, 5, False, False, 0, False, None, None, None),
+    (1, 10, 1, 8,'steps', "incremental", "cdf", False, "spn", True, 1, False, False, 0, False, None, None, None),
+    (2, 3, 2, 8,'steps', "incremental", "cdf", True, "spn", True, 2, False, False, 0, False, None, None, None),
+    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, False, False, 0, False, None, None, None),
+    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, False, False, 0, False, "bps", None, None),
+    # Test the case where the radar image contains no rain.
+    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, True, False, 0, False, None, None, None),
+    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, True, False, 0, False, None, None, None),
+    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, True, False, 0, True, None, None, None),
+    # Test the case where the NWP fields contain no rain.
+    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, False, True, 0, False, None, None, None),
+    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, False, True, 0, True, None, None, None),
     # Test the case where both the radar image and the NWP fields contain no rain.
-    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, True, True, 0, False, None, None),
-    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, True, True, 0, False, None, None),
-    (5, 3, 5, 6,'steps', "obs", "mean", True, "spn", True, 5, True, True, 0, False, None, None),
+    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, True, True, 0, False, None, None, None),
+    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, True, True, 0, False, None, None, None),
+    (5, 3, 5, 6,'steps', "obs", "mean", True, "spn", True, 5, True, True, 0, False, None, None, None),
+    # Test cases where we apply timestep_start_full_nwp_weight
+    (1, 10, 2, 6,'steps', "incremental", "cdf", False, "bps", False, 2, False, False, 0, True, None, None, 5),
+    (1, 10, 2, 6,'steps', "incremental", "cdf", False, "spn", False, 2, False, False, 0, False, None, None, 5),
     # Test for smooth radar mask
-    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, False, False, 80, False, None, None),
-    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, False, False, 80, False, None, None),
-    (5, 3, 5, 6,'steps', "obs", "mean", False, "spn", False, 5, False, False, 80, False, None, None),
-    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, False, True, 80, False, None, None),
-    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, True, False, 80, True, None, None),
-    (5, 3, 5, 6,'steps', "obs", "mean", False, "spn", False, 5, True, True, 80, False, None, None),
-    (5, [1, 2, 3], 5, 6,'steps', "obs", "mean", False, "spn", False, 5, True, True, 80, False, None, None),
-    (5, [1, 3], 5, 6,'steps', "obs", "mean", False, "spn", False, 5, True, True, 80, False, None, None),
+    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, False, False, 80, False, None, None, None),
+    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, False, False, 80, False, None, None, None),
+    (5, 3, 5, 6,'steps', "obs", "mean", False, "spn", False, 5, False, False, 80, False, None, None, None),
+    (1, 3, 6, 8,'steps', None, None, False, "spn", True, 6, False, True, 80, False, None, None, None),
+    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "spn", False, 5, True, False, 80, True, None, None, None),
+    (5, 3, 5, 6,'steps', "obs", "mean", False, "spn", False, 5, True, True, 80, False, None, None, None),
+    (5, [1, 2, 3], 5, 6,'steps', "obs", "mean", False, "spn", False, 5, True, True, 80, False, None, None, None),
+    (5, [1, 3], 5, 6,'steps', "obs", "mean", False, "spn", False, 5, True, True, 80, False, None, None, None),
     # Test the usage of a max_mask_rim in the mask_kwargs
-    (1, 3, 6, 8,'steps', None, None, False, "bps", True, 6, False, False, 80, False, None, 40),
-    (5, 3, 5, 6,'steps', "obs", "mean", False, "bps", False, 5, False, False, 80, False, None, 40),
-    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "bps", False, 5, False, False, 80, False, None, 25),
-    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "bps", False, 5, False, False, 80, False, None, 40),
-    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "bps", False, 5, False, False, 80, False, None, 60),
+    (1, 3, 6, 8,'steps', None, None, False, "bps", True, 6, False, False, 80, False, None, 40, None),
+    (5, 3, 5, 6,'steps', "obs", "mean", False, "bps", False, 5, False, False, 80, False, None, 40, None),
+    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "bps", False, 5, False, False, 80, False, None, 25, None),
+    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "bps", False, 5, False, False, 80, False, None, 40, None),
+    (5, 3, 5, 6,'steps', "incremental", "cdf", False, "bps", False, 5, False, False, 80, False, None, 60, None),
     #Test the externally provided nowcast
-    (1, 10, 1, 8,'external_nowcast_det', None, None, False, "spn", True, 1, False, False, 0, False, None, None),
-    (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "bps", True, 1, False, False, 0, False, None, None),
-    (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "spn", True, 1, False, False, 80, False, None, None),
-    (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "bps", True, 1, True, False, 0, False, None, None),
-    (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "spn", True, 1, False, True, 0, False, None, None),
-    (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "bps", True, 1, True, True, 0, False, None, None),
-    (1, 10, 1, 8,'external_nowcast_det', "incremental", "cdf", False, "spn", True, 1, False, False, 0, True, None, None),
-    (1, 10, 1, 8,'external_nowcast_det', "incremental", "obs", False, "bps", True, 1, False, False, 0, False, None, None),
-    (5, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None),
-    (5, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None),
-    (1, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None),
-    (1, 10, 1, 8,'external_nowcast_ens', "incremental", "cdf", False, "bps", True, 5, False, False, 0, False, None, None),
-    (5, 10, 1, 8,'external_nowcast_ens', "incremental", "obs", False, "spn", True, 5, False, False, 0, False, None, None)
+    (1, 10, 1, 8,'external_nowcast_det', None, None, False, "spn", True, 1, False, False, 0, False, None, None, None),
+    (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "bps", True, 1, False, False, 0, False, None, None, None),
+    (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "spn", True, 1, False, False, 80, False, None, None, None),
+    (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "bps", True, 1, True, False, 0, False, None, None, None),
+    (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "spn", True, 1, False, True, 0, False, None, None, None),
+    (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "bps", True, 1, True, True, 0, False, None, None, None),
+    (1, 10, 1, 8,'external_nowcast_det', "incremental", "cdf", False, "spn", True, 1, False, False, 0, True, None, None, None),
+    (1, 10, 1, 8,'external_nowcast_det', "incremental", "obs", False, "bps", True, 1, False, False, 0, False, None, None, None),
+    (1, 10, 1, 8,'external_nowcast_det', "incremental", None, False, "bps", True, 1, False, False, 0, False, None, None, 5),
+    (5, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None, None),
+    (5, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None, None),
+    (1, 10, 5, 8,'external_nowcast_ens', "incremental", None, False, "spn", True, 5, False, False, 0, False, None, None, None),
+    (1, 10, 1, 8,'external_nowcast_ens', "incremental", "cdf", False, "bps", True, 5, False, False, 0, False, None, None, None),
+    (5, 10, 1, 8,'external_nowcast_ens', "incremental", "obs", False, "spn", True, 5, False, False, 0, False, None, None, None),
+    (1, 10, 5, 8,'external_nowcast_ens', "incremental", "cdf", False, "bps", True, 5, False, False, 0, False, None, None, 5)
 ]
 
 # fmt:on
@@ -90,6 +95,7 @@
     "resample_distribution",
     "vel_pert_method",
     "max_mask_rim",
+    "timestep_start_full_nwp_weight",
 )
 
 
@@ -112,6 +118,7 @@ def test_steps_blending(
     resample_distribution,
     vel_pert_method,
     max_mask_rim,
+    timestep_start_full_nwp_weight,
 ):
     pytest.importorskip("cv2")
 
@@ -350,7 +357,7 @@ def test_steps_blending(
     assert nwp_velocity.ndim == 5, "nwp_velocity must be a five-dimensional array"
 
     ###
-    # The nowcasting
+    # The blending
     ###
     precip_forecast = blending.steps.forecast(
         precip=radar_precip,
@@ -373,6 +380,7 @@ def test_steps_blending(
         ar_order=2,
         vel_pert_method=vel_pert_method,
         weights_method=weights_method,
+        timestep_start_full_nwp_weight=timestep_start_full_nwp_weight,
         conditional=False,
         probmatching_method=probmatching_method,
         mask_method=mask_method,
diff --git a/pysteps/tests/test_datasets.py b/pysteps/tests/test_datasets.py
index 72d6f3660..93eefaaaa 100644
--- a/pysteps/tests/test_datasets.py
+++ b/pysteps/tests/test_datasets.py
@@ -13,7 +13,6 @@
 )
 from pysteps.exceptions import DirectoryNotEmpty
 
-
 _datasets_opt_deps = dict(
     fmi=["pyproj"],
     mch=["PIL"],
diff --git a/pysteps/tests/test_noise_fftgenerators.py b/pysteps/tests/test_noise_fftgenerators.py
index cecaf8ca4..272497649 100644
--- a/pysteps/tests/test_noise_fftgenerators.py
+++ b/pysteps/tests/test_noise_fftgenerators.py
@@ -3,7 +3,6 @@
 from pysteps.noise import fftgenerators
 from pysteps.tests.helpers import get_precipitation_fields
 
-
 PRECIP = get_precipitation_fields(
     num_prev_files=0,
     num_next_files=0,
diff --git a/pysteps/tests/test_nowcasts_linda.py b/pysteps/tests/test_nowcasts_linda.py
index da5369e7e..a9d599972 100644
--- a/pysteps/tests/test_nowcasts_linda.py
+++ b/pysteps/tests/test_nowcasts_linda.py
@@ -7,7 +7,6 @@
 from pysteps.nowcasts.linda import forecast
 from pysteps.tests.helpers import get_precipitation_fields
 
-
 linda_arg_names = (
     "timesteps",
     "add_perturbations",
diff --git a/pysteps/tests/test_plt_animate.py b/pysteps/tests/test_plt_animate.py
index 6b9892ae3..39df985c3 100644
--- a/pysteps/tests/test_plt_animate.py
+++ b/pysteps/tests/test_plt_animate.py
@@ -9,7 +9,6 @@
 from pysteps.tests.helpers import get_precipitation_fields
 from pysteps.visualization.animations import animate
 
-
 PRECIP, METADATA = get_precipitation_fields(
     num_prev_files=2,
     num_next_files=0,
diff --git a/pysteps/tests/test_plt_motionfields.py b/pysteps/tests/test_plt_motionfields.py
index d0c7e6414..c6bb548d9 100644
--- a/pysteps/tests/test_plt_motionfields.py
+++ b/pysteps/tests/test_plt_motionfields.py
@@ -8,7 +8,6 @@
 from pysteps.visualization import plot_precip_field, quiver, streamplot
 from pysteps.tests.helpers import get_precipitation_fields
 
-
 arg_names_quiver = (
     "source",
     "axis",
diff --git a/pysteps/tests/test_plugins_support.py b/pysteps/tests/test_plugins_support.py
index 280bc5c75..3e063d0be 100644
--- a/pysteps/tests/test_plugins_support.py
+++ b/pysteps/tests/test_plugins_support.py
@@ -4,13 +4,13 @@
 
 https://github.com/pySTEPS/cookiecutter-pysteps-plugin
 """
+
 import os
 import pytest
 import subprocess
 import sys
 import tempfile
 
-
 __ = pytest.importorskip("cookiecutter")
 from cookiecutter.main import cookiecutter
 
diff --git a/pysteps/tests/test_utils_interpolate.py b/pysteps/tests/test_utils_interpolate.py
index eee0e6b2c..6b2c55443 100644
--- a/pysteps/tests/test_utils_interpolate.py
+++ b/pysteps/tests/test_utils_interpolate.py
@@ -4,7 +4,6 @@
 
 from pysteps.utils import get_method
 
-
 interp_methods = (
     "idwinterp2d",
     "rbfinterp2d",
diff --git a/pysteps/tests/test_utils_pca.py b/pysteps/tests/test_utils_pca.py
index be8488286..4b1fcb09d 100644
--- a/pysteps/tests/test_utils_pca.py
+++ b/pysteps/tests/test_utils_pca.py
@@ -4,7 +4,6 @@
 import numpy as np
 from pysteps.utils import pca
 
-
 pca_arg_values = (
     (10, 10),
     (20, 20),
diff --git a/pysteps/tests/test_verification_detcatscores.py b/pysteps/tests/test_verification_detcatscores.py
index ba14b23e6..2126df6bb 100644
--- a/pysteps/tests/test_verification_detcatscores.py
+++ b/pysteps/tests/test_verification_detcatscores.py
@@ -6,7 +6,6 @@
 
 from pysteps.verification import det_cat_fct
 
-
 # CREATE A LARGE DATASET TO MATCH
 # EXAMPLES IN
 # http://www.cawcr.gov.au/projects/verification/
diff --git a/pysteps/utils/cleansing.py b/pysteps/utils/cleansing.py
index 701ec78bb..7c4276c57 100644
--- a/pysteps/utils/cleansing.py
+++ b/pysteps/utils/cleansing.py
@@ -11,6 +11,7 @@
     decluster
     detect_outliers
 """
+
 import warnings
 
 import numpy as np
diff --git a/pysteps/utils/images.py b/pysteps/utils/images.py
index ab7bd8277..6f2376e2f 100644
--- a/pysteps/utils/images.py
+++ b/pysteps/utils/images.py
@@ -10,6 +10,7 @@
 
     morph_opening
 """
+
 import numpy as np
 from numpy.ma.core import MaskedArray
 
diff --git a/pysteps/utils/reprojection.py b/pysteps/utils/reprojection.py
index 1d070550d..8acd34698 100644
--- a/pysteps/utils/reprojection.py
+++ b/pysteps/utils/reprojection.py
@@ -11,6 +11,7 @@
 
     reproject_grids
 """
+
 from pysteps.exceptions import MissingOptionalDependency
 from scipy.interpolate import griddata
 
diff --git a/pysteps/verification/salscores.py b/pysteps/verification/salscores.py
index 75abe6124..f2d97b77c 100644
--- a/pysteps/verification/salscores.py
+++ b/pysteps/verification/salscores.py
@@ -13,6 +13,7 @@
     sal_amplitude
     sal_location
 """
+
 from math import sqrt, hypot
 
 import numpy as np
diff --git a/pysteps/visualization/precipfields.py b/pysteps/visualization/precipfields.py
index 8b03ed807..69bf384ba 100644
--- a/pysteps/visualization/precipfields.py
+++ b/pysteps/visualization/precipfields.py
@@ -11,6 +11,7 @@
     plot_precip_field
     get_colormap
 """
+
 import copy
 import warnings