    def __perturb_blend_and_advect_extrapolation_and_noise_to_current_timestep(
        self, t, j, worker_state
    ):
        """
        Apply perturbations, blend motion fields, and advect extrapolated and noise
        cascades to the current time step (or sub-timesteps). This step ensures
        realistic motion updates in nowcasting.
        """
        # Settings and initialize the output
        extrap_kwargs_ = worker_state.extrapolation_kwargs.copy()
        extrap_kwargs_noise = worker_state.extrapolation_kwargs.copy()
        extrap_kwargs_pb = worker_state.extrapolation_kwargs.copy()
        velocity_perturbations_extrapolation = self.__velocity
        # The following should be accessible after this function
        worker_state.precip_extrapolated_decomp = []
        worker_state.noise_extrapolated_decomp = []
        worker_state.precip_extrapolated_probability_matching = []

        # Extrapolate per sub time step
        for t_sub in worker_state.subtimesteps:
            if t_sub > 0:
                t_diff_prev_subtimestep_int = t_sub - int(t_sub)
                if t_diff_prev_subtimestep_int > 0.0:
                    precip_forecast_cascade_subtimestep = [
                        (1.0 - t_diff_prev_subtimestep_int)
                        * worker_state.precip_cascades_prev_subtimestep[j][i][-1, :]
                        + t_diff_prev_subtimestep_int
                        * worker_state.precip_cascades[j][i][-1, :]
                        for i in range(self.__config.n_cascade_levels)
                    ]
                    noise_cascade_subtimestep = [
                        (1.0 - t_diff_prev_subtimestep_int)
                        * worker_state.cascade_noise_prev_subtimestep[j][i][-1, :]
                        + t_diff_prev_subtimestep_int
                        * worker_state.precip_noise_cascades[j][i][-1, :]
                        for i in range(self.__config.n_cascade_levels)
                    ]

                else:
                    precip_forecast_cascade_subtimestep = [
                        worker_state.precip_cascades_prev_subtimestep[j][i][-1, :]
                        for i in range(self.__config.n_cascade_levels)
                    ]
                    noise_cascade_subtimestep = [
                        worker_state.cascade_noise_prev_subtimestep[j][i][-1, :]
                        for i in range(self.__config.n_cascade_levels)
                    ]

                precip_forecast_cascade_subtimestep = np.stack(
                    precip_forecast_cascade_subtimestep
                )
                noise_cascade_subtimestep = np.stack(noise_cascade_subtimestep)

                t_diff_prev_subtimestep = t_sub - worker_state.time_prev_timestep[j]
                worker_state.leadtime_since_start_forecast[j] += t_diff_prev_subtimestep

                # compute the perturbed motion field - include the NWP
                # velocities and the weights. Note that we only perturb
                # the extrapolation velocity field, as the NWP velocity
                # field is present per time step
                if self.__config.velocity_perturbation_method is not None:
                    velocity_perturbations_extrapolation = (
                        self.__velocity
                        + self.__params.generate_velocity_noise(
                            self.__params.velocity_perturbations[j],
                            worker_state.leadtime_since_start_forecast[j]
                            * self.__config.timestep,
                        )
                    )

                # Stack the perturbed extrapolation and the NWP velocities
                if self.__config.blend_nwp_members:
                    velocity_stack_all = np.concatenate(
                        (
                            velocity_perturbations_extrapolation[None, :, :, :],
                            worker_state.velocity_models_timestep,
                        ),
                        axis=0,
                    )
                else:
                    velocity_models = worker_state.velocity_models_timestep[j]
                    velocity_stack_all = np.concatenate(
                        (
                            velocity_perturbations_extrapolation[None, :, :, :],
                            velocity_models[None, :, :, :],
                        ),
                        axis=0,
                    )
                    velocity_models = None

                # Obtain a blended optical flow, using the weights of the
                # second cascade following eq. 24 in BPS2006
                velocity_blended = blending.utils.blend_optical_flows(
                    flows=velocity_stack_all,
                    weights=worker_state.weights[
                        :-1, 1
                    ],  # [(extr_field, n_model_fields), cascade_level=2]
                )

                # Extrapolate both cascades to the next time step
                # First recompose the cascade, advect it and decompose it again
                # This is needed to remove the interpolation artifacts.
                # In addition, the number of extrapolations is greatly reduced
                # A. Radar Rain
                precip_forecast_recomp_subtimestep = blending.utils.recompose_cascade(
                    combined_cascade=precip_forecast_cascade_subtimestep,
                    combined_mean=worker_state.mean_extrapolation,
                    combined_sigma=worker_state.std_extrapolation,
                )
                # Make sure we have values outside the mask
                if self.__params.zero_precip_radar:
                    precip_forecast_recomp_subtimestep = np.nan_to_num(
                        precip_forecast_recomp_subtimestep,
                        copy=True,
                        nan=self.__params.precip_zerovalue,
                        posinf=self.__params.precip_zerovalue,
                        neginf=self.__params.precip_zerovalue,
                    )
                # Put back the mask
                precip_forecast_recomp_subtimestep[self.__params.domain_mask] = np.nan
                worker_state.extrapolation_kwargs["displacement_prev"] = (
                    worker_state.previous_displacement[j]
                )
                (
                    precip_forecast_extrapolated_recomp_subtimestep_temp,
                    worker_state.previous_displacement[j],
                ) = self.__params.extrapolation_method(
                    precip_forecast_recomp_subtimestep,
                    velocity_blended,
                    [t_diff_prev_subtimestep],
                    allow_nonfinite_values=True,
                    **worker_state.extrapolation_kwargs,
                )
                precip_extrapolated_recomp_subtimestep = (
                    precip_forecast_extrapolated_recomp_subtimestep_temp[0].copy()
                )
                temp_mask = ~np.isfinite(precip_extrapolated_recomp_subtimestep)
                # TODO: WHERE DO CAN I FIND THIS -15.0
                precip_extrapolated_recomp_subtimestep[
                    ~np.isfinite(precip_extrapolated_recomp_subtimestep)
                ] = self.__params.precip_zerovalue
                precip_extrapolated_decomp = self.__params.decomposition_method(
                    precip_extrapolated_recomp_subtimestep,
                    self.__params.bandpass_filter,
                    mask=self.__params.mask_threshold,
                    fft_method=self.__params.fft,
                    output_domain=self.__config.domain,
                    normalize=True,
                    compute_stats=True,
                    compact_output=True,
                )["cascade_levels"]
                # Make sure we have values outside the mask
                if self.__params.zero_precip_radar:
                    precip_extrapolated_decomp = np.nan_to_num(
                        precip_extrapolated_decomp,
                        copy=True,
                        nan=np.nanmin(precip_forecast_cascade_subtimestep),
                        posinf=np.nanmin(precip_forecast_cascade_subtimestep),
                        neginf=np.nanmin(precip_forecast_cascade_subtimestep),
                    )
                for i in range(self.__config.n_cascade_levels):
                    precip_extrapolated_decomp[i][temp_mask] = np.nan
                # B. Noise
                noise_cascade_subtimestep_recomp = blending.utils.recompose_cascade(
                    combined_cascade=noise_cascade_subtimestep,
                    combined_mean=worker_state.precip_mean_noise[j],
                    combined_sigma=worker_state.precip_std_noise[j],
                )
                extrap_kwargs_noise["displacement_prev"] = (
                    worker_state.previous_displacement_noise_cascade[j]
                )
                extrap_kwargs_noise["map_coordinates_mode"] = "wrap"
                (
                    noise_extrapolated_recomp_temp,
                    worker_state.previous_displacement_noise_cascade[j],
                ) = self.__params.extrapolation_method(
                    noise_cascade_subtimestep_recomp,
                    velocity_blended,
                    [t_diff_prev_subtimestep],
                    allow_nonfinite_values=True,
                    **extrap_kwargs_noise,
                )
                noise_extrapolated_recomp = noise_extrapolated_recomp_temp[0].copy()
                noise_extrapolated_decomp = self.__params.decomposition_method(
                    noise_extrapolated_recomp,
                    self.__params.bandpass_filter,
                    mask=self.__params.mask_threshold,
                    fft_method=self.__params.fft,
                    output_domain=self.__config.domain,
                    normalize=True,
                    compute_stats=True,
                    compact_output=True,
                )["cascade_levels"]
                for i in range(self.__config.n_cascade_levels):
                    noise_extrapolated_decomp[i] *= self.__params.noise_std_coeffs[i]

                # Append the results to the output lists
                worker_state.precip_extrapolated_decomp.append(
                    precip_extrapolated_decomp.copy()
                )
                worker_state.noise_extrapolated_decomp.append(
                    noise_extrapolated_decomp.copy()
                )
                precip_forecast_cascade_subtimestep = None
                precip_forecast_recomp_subtimestep = None
                precip_forecast_extrapolated_recomp_subtimestep_temp = None
                precip_extrapolated_recomp_subtimestep = None
                precip_extrapolated_decomp = None
                noise_cascade_subtimestep = None
                noise_cascade_subtimestep_recomp = None
                noise_extrapolated_recomp_temp = None
                noise_extrapolated_recomp = None
                noise_extrapolated_decomp = None

                # Finally, also extrapolate the initial radar rainfall
                # field. This will be blended with the rainfall field(s)
                # of the (NWP) model(s) for Lagrangian blended prob. matching
                # min_R = np.min(precip)
                extrap_kwargs_pb["displacement_prev"] = (
                    worker_state.previous_displacement_prob_matching[j]
                )
                # Apply the domain mask to the extrapolation component
                precip_forecast_temp_for_probability_matching = self.__precip.copy()
                precip_forecast_temp_for_probability_matching[
                    self.__params.domain_mask
                ] = np.nan

                (
                    precip_forecast_extrapolated_probability_matching_temp,
                    worker_state.previous_displacement_prob_matching[j],
                ) = self.__params.extrapolation_method(
                    precip_forecast_temp_for_probability_matching,
                    velocity_blended,
                    [t_diff_prev_subtimestep],
                    allow_nonfinite_values=True,
                    **extrap_kwargs_pb,
                )

                worker_state.precip_extrapolated_probability_matching.append(
                    precip_forecast_extrapolated_probability_matching_temp[0]
                )

                worker_state.time_prev_timestep[j] = t_sub

        if len(worker_state.precip_extrapolated_decomp) > 0:
            worker_state.precip_extrapolated_decomp = np.stack(
                worker_state.precip_extrapolated_decomp
            )
            worker_state.noise_extrapolated_decomp = np.stack(
                worker_state.noise_extrapolated_decomp
            )
            worker_state.precip_extrapolated_probability_matching = np.stack(
                worker_state.precip_extrapolated_probability_matching
            )

        # advect the forecast field by one time step if no subtimesteps in the
        # current interval were found
        if not worker_state.subtimesteps:
            t_diff_prev_subtimestep = t + 1 - worker_state.time_prev_timestep[j]
            worker_state.leadtime_since_start_forecast[j] += t_diff_prev_subtimestep

            # compute the perturbed motion field - include the NWP
            # velocities and the weights
            if self.__config.velocity_perturbation_method is not None:
                velocity_perturbations_extrapolation = (
                    self.__velocity
                    + self.__params.generate_velocity_noise(
                        self.__params.velocity_perturbations[j],
                        worker_state.leadtime_since_start_forecast[j]
                        * self.__config.timestep,
                    )
                )

            # Stack the perturbed extrapolation and the NWP velocities
            if self.__config.blend_nwp_members:
                velocity_stack_all = np.concatenate(
                    (
                        velocity_perturbations_extrapolation[None, :, :, :],
                        worker_state.velocity_models_timestep,
                    ),
                    axis=0,
                )
            else:
                velocity_models = worker_state.velocity_models_timestep[j]
                velocity_stack_all = np.concatenate(
                    (
                        velocity_perturbations_extrapolation[None, :, :, :],
                        velocity_models[None, :, :, :],
                    ),
                    axis=0,
                )
                velocity_models = None

            # Obtain a blended optical flow, using the weights of the
            # second cascade following eq. 24 in BPS2006
            velocity_blended = blending.utils.blend_optical_flows(
                flows=velocity_stack_all,
                weights=worker_state.weights[
                    :-1, 1
                ],  # [(extr_field, n_model_fields), cascade_level=2]
            )

            # Extrapolate the extrapolation and noise cascade

            extrap_kwargs_["displacement_prev"] = worker_state.previous_displacement[j]
            extrap_kwargs_noise["displacement_prev"] = (
                worker_state.previous_displacement_noise_cascade[j]
            )
            extrap_kwargs_noise["map_coordinates_mode"] = "wrap"

            (
                _,
                worker_state.previous_displacement[j],
            ) = self.__params.extrapolation_method(
                None,
                velocity_blended,
                [t_diff_prev_subtimestep],
                allow_nonfinite_values=True,
                **extrap_kwargs_,
            )

            (
                _,
                worker_state.previous_displacement_noise_cascade[j],
            ) = self.__params.extrapolation_method(
                None,
                velocity_blended,
                [t_diff_prev_subtimestep],
                allow_nonfinite_values=True,
                **extrap_kwargs_noise,
            )

            # Also extrapolate the radar observation, used for the probability
            # matching and post-processing steps
            extrap_kwargs_pb["displacement_prev"] = (
                worker_state.previous_displacement_prob_matching[j]
            )
            (
                _,
                worker_state.previous_displacement_prob_matching[j],
            ) = self.__params.extrapolation_method(
                None,
                velocity_blended,
                [t_diff_prev_subtimestep],
                allow_nonfinite_values=True,
                **extrap_kwargs_pb,
            )

            worker_state.time_prev_timestep[j] = t + 1

        worker_state.precip_cascades_prev_subtimestep[j] = worker_state.precip_cascades[
            j
        ]
        worker_state.cascade_noise_prev_subtimestep[j] = (
            worker_state.precip_noise_cascades[j]
        )
Apply perturbations, blend motion fields, and advect extrapolated and noise
cascades to the current time step (or sub-timesteps). This step ensures
realistic motion updates in nowcasting.
pysteps/blending/steps.py