    def __post_process_output(
        self, j, final_blended_forecast_single_member, worker_state
    ):
        """
        Apply post-processing steps to refine the final blended forecast. This
        involves masking, filling missing data with the blended NWP forecast,
        and applying probability matching to ensure consistency.

        **Steps:**

        1. **Use Mask and Fill Missing Data:**
           - Areas without reliable radar extrapolation are filled using the
             blended NWP forecast to maintain spatial coherence.

        2. **Lagrangian Blended Probability Matching:**
           - Uses the latest extrapolated radar rainfall field blended with
             the NWP model(s) forecast as a reference.
           - Ensures that the statistical distribution of the final forecast
             remains consistent with the benchmark dataset.

        3. **Blend the Extrapolated Rainfall Field with NWP Forecasts:**
           - The extrapolated rainfall field is used only for post-processing.
           - The forecast is blended using predefined weights at scale level 2.
           - This ensures that both extrapolated and modeled precipitation
             contribute appropriately to the final output.

        4. **Apply Probability Matching:**
           - Adjusts the final precipitation distribution using either empirical
             cumulative distribution functions (CDF) or mean adjustments to
             match the reference dataset.

        The final processed forecast is stored in `final_blended_forecast_single_member`.
        """

        weights_probability_matching = worker_state.weights[
            :-1, 1
        ]  # Weights without noise, level 2
        weights_probability_matching_normalized = weights_probability_matching / np.sum(
            weights_probability_matching
        )
        # And the weights for outside the radar domain
        weights_probability_matching_mod_only = worker_state.weights_model_only[
            :-1, 1
        ]  # Weights without noise, level 2
        weights_probability_matching_normalized_mod_only = (
            weights_probability_matching_mod_only
            / np.sum(weights_probability_matching_mod_only)
        )
        # Stack the fields
        if self.__config.blend_nwp_members:
            precip_forecast_probability_matching_final = np.concatenate(
                (
                    worker_state.precip_extrapolated_probability_matching[
                        None, worker_state.subtimestep_index
                    ],
                    worker_state.precip_models_timestep,
                ),
                axis=0,
            )
        else:
            precip_forecast_probability_matching_final = np.concatenate(
                (
                    worker_state.precip_extrapolated_probability_matching[
                        None, worker_state.subtimestep_index
                    ],
                    worker_state.precip_models_timestep[None, j],
                ),
                axis=0,
            )
        # Blend it
        precip_forecast_probability_matching_blended = np.sum(
            weights_probability_matching_normalized.reshape(
                weights_probability_matching_normalized.shape[0], 1, 1
            )
            * precip_forecast_probability_matching_final,
            axis=0,
        )
        if self.__config.blend_nwp_members:
            precip_forecast_probability_matching_blended_mod_only = np.sum(
                weights_probability_matching_normalized_mod_only.reshape(
                    weights_probability_matching_normalized_mod_only.shape[0],
                    1,
                    1,
                )
                * worker_state.precip_models_timestep,
                axis=0,
            )
        else:
            precip_forecast_probability_matching_blended_mod_only = (
                worker_state.precip_models_timestep[j]
            )

        # The extrapolation components are NaN outside the advected
        # radar domain. This results in NaN values in the blended
        # forecast outside the radar domain. Therefore, fill these
        # areas with the "..._mod_only" blended forecasts, consisting
        # of the NWP and noise components.

        nan_indices = np.isnan(worker_state.final_blended_forecast_recomposed)
        if self.__config.smooth_radar_mask_range != 0:
            # Compute the smooth dilated mask
            new_mask = blending.utils.compute_smooth_dilated_mask(
                nan_indices,
                max_padding_size_in_px=self.__config.smooth_radar_mask_range,
            )

            # Ensure mask values are between 0 and 1
            mask_model = np.clip(new_mask, 0, 1)
            mask_radar = np.clip(1 - new_mask, 0, 1)

            # Handle NaNs in precip_forecast_new and precip_forecast_new_mod_only by setting NaNs to 0 in the blending step
            precip_forecast_recomposed_mod_only_no_nan = np.nan_to_num(
                worker_state.final_blended_forecast_recomposed_mod_only, nan=0
            )
            precip_forecast_recomposed_no_nan = np.nan_to_num(
                worker_state.final_blended_forecast_recomposed, nan=0
            )

            # Perform the blending of radar and model inside the radar domain using a weighted combination
            worker_state.final_blended_forecast_recomposed = np.nansum(
                [
                    mask_model * precip_forecast_recomposed_mod_only_no_nan,
                    mask_radar * precip_forecast_recomposed_no_nan,
                ],
                axis=0,
            )

            precip_forecast_probability_matching_blended = np.nansum(
                [
                    precip_forecast_probability_matching_blended * mask_radar,
                    precip_forecast_probability_matching_blended_mod_only * mask_model,
                ],
                axis=0,
            )
        else:
            worker_state.final_blended_forecast_recomposed[nan_indices] = (
                worker_state.final_blended_forecast_recomposed_mod_only[nan_indices]
            )
            nan_indices = np.isnan(precip_forecast_probability_matching_blended)
            precip_forecast_probability_matching_blended[nan_indices] = (
                precip_forecast_probability_matching_blended_mod_only[nan_indices]
            )

        # Finally, fill the remaining nan values, if present, with
        # the minimum value in the forecast
        nan_indices = np.isnan(worker_state.final_blended_forecast_recomposed)
        worker_state.final_blended_forecast_recomposed[nan_indices] = np.nanmin(
            worker_state.final_blended_forecast_recomposed
        )
        nan_indices = np.isnan(precip_forecast_probability_matching_blended)
        precip_forecast_probability_matching_blended[nan_indices] = np.nanmin(
            precip_forecast_probability_matching_blended
        )

        # Apply the masking and prob. matching
        precip_field_mask_temp = None
        if self.__config.mask_method is not None:
            # apply the precipitation mask to prevent generation of new
            # precipitation into areas where it was not originally
            # observed
            precip_forecast_min_value = (
                worker_state.final_blended_forecast_recomposed.min()
            )
            if self.__config.mask_method == "incremental":
                # The incremental mask is slightly different from
                # the implementation in the non-blended steps.py, as
                # it is not based on the last forecast, but instead
                # on R_pm_blended. Therefore, the buffer does not
                # increase over time.
                # Get the mask for this forecast
                precip_field_mask = (
                    precip_forecast_probability_matching_blended
                    >= self.__config.precip_threshold
                )
                # Buffer the mask
                # Convert the precipitation field mask into an 8-bit unsigned integer mask
                obs_mask_uint8 = precip_field_mask.astype("uint8")

                # Perform an initial binary dilation using the provided structuring element
                dilated_mask = binary_dilation(obs_mask_uint8, self.__params.struct)

                # Create a binary structure element for incremental dilations
                struct_element = generate_binary_structure(2, 1)

                # Initialize a floating-point mask to accumulate dilations for a smooth transition
                accumulated_mask = dilated_mask.astype(float)

                # Iteratively dilate the mask and accumulate the results to create a grayscale rim
                for _ in range(self.__params.mask_rim):
                    dilated_mask = binary_dilation(dilated_mask, struct_element)
                    accumulated_mask += dilated_mask

                # Normalize the accumulated mask values between 0 and 1
                precip_field_mask = accumulated_mask / np.max(accumulated_mask)
                # Get the final mask
                worker_state.final_blended_forecast_recomposed = (
                    precip_forecast_min_value
                    + (
                        worker_state.final_blended_forecast_recomposed
                        - precip_forecast_min_value
                    )
                    * precip_field_mask
                )
                precip_field_mask_temp = (
                    worker_state.final_blended_forecast_recomposed
                    > precip_forecast_min_value
                )
            elif self.__config.mask_method == "obs":
                # The mask equals the most recent benchmark
                # rainfall field
                precip_field_mask_temp = (
                    precip_forecast_probability_matching_blended
                    >= self.__config.precip_threshold
                )

            # Set to min value outside of mask
            worker_state.final_blended_forecast_recomposed[~precip_field_mask_temp] = (
                precip_forecast_min_value
            )

        # If probmatching_method is not None, resample the distribution from
        # both the extrapolation cascade and the model (NWP) cascade and use
        # that for the probability matching.
        if (
            self.__config.probmatching_method is not None
            and self.__config.resample_distribution
        ):
            arr1 = worker_state.precip_extrapolated_probability_matching[
                worker_state.subtimestep_index
            ]
            arr2 = worker_state.precip_models_timestep[j]
            # resample weights based on cascade level 2.
            # Areas where one of the fields is nan are not included.
            precip_forecast_probability_matching_resampled = (
                probmatching.resample_distributions(
                    first_array=arr1,
                    second_array=arr2,
                    probability_first_array=weights_probability_matching_normalized[0],
                    randgen=worker_state.randgen_probmatching[j],
                )
            )
        else:
            precip_forecast_probability_matching_resampled = (
                precip_forecast_probability_matching_blended.copy()
            )

        if self.__config.probmatching_method == "cdf":
            # nan indices in the extrapolation nowcast
            nan_indices = np.isnan(
                worker_state.precip_extrapolated_probability_matching[
                    worker_state.subtimestep_index
                ]
            )
            # Adjust the CDF of the forecast to match the resampled distribution combined from
            # extrapolation and model fields.
            # Rainfall outside the pure extrapolation domain is not taken into account.
            if np.any(np.isfinite(worker_state.final_blended_forecast_recomposed)):
                worker_state.final_blended_forecast_recomposed = (
                    probmatching.nonparam_match_empirical_cdf(
                        worker_state.final_blended_forecast_recomposed,
                        precip_forecast_probability_matching_resampled,
                        nan_indices,
                    )
                )
                precip_forecast_probability_matching_resampled = None
        elif self.__config.probmatching_method == "mean":
            # Use R_pm_blended as benchmark field and
            mean_probabiltity_matching_forecast = np.mean(
                precip_forecast_probability_matching_resampled[
                    precip_forecast_probability_matching_resampled
                    >= self.__config.precip_threshold
                ]
            )
            no_rain_mask = (
                worker_state.final_blended_forecast_recomposed
                >= self.__config.precip_threshold
            )
            mean_precip_forecast = np.mean(
                worker_state.final_blended_forecast_recomposed[no_rain_mask]
            )
            worker_state.final_blended_forecast_recomposed[no_rain_mask] = (
                worker_state.final_blended_forecast_recomposed[no_rain_mask]
                - mean_precip_forecast
                + mean_probabiltity_matching_forecast
            )
            precip_forecast_probability_matching_resampled = None

        final_blended_forecast_single_member.append(
            worker_state.final_blended_forecast_recomposed
        )
        return final_blended_forecast_single_member
Apply post-processing steps to refine the final blended forecast. This
involves masking, filling missing data with the blended NWP forecast,
and applying probability matching to ensure consistency.

**Steps:**

1. **Use Mask and Fill Missing Data:**
- Areas without reliable radar extrapolation are filled using the
blended NWP forecast to maintain spatial coherence.

2. **Lagrangian Blended Probability Matching:**
- Uses the latest extrapolated radar rainfall field blended with
the NWP model(s) forecast as a reference.
- Ensures that the statistical distribution of the final forecast
remains consistent with the benchmark dataset.

3. **Blend the Extrapolated Rainfall Field with NWP Forecasts:**
- The extrapolated rainfall field is used only for post-processing.
- The forecast is blended using predefined weights at scale level 2.
- This ensures that both extrapolated and modeled precipitation
contribute appropriately to the final output.

4. **Apply Probability Matching:**
- Adjusts the final precipitation distribution using either empirical
cumulative distribution functions (CDF) or mean adjustments to
match the reference dataset.

The final processed forecast is stored in `final_blended_forecast_single_member`.
pysteps/blending/steps.py