diff --git a/profiles_eval/configs/dlwindstat.json b/profiles_eval/configs/dlwindstat.json
index d46eb54..670dcce 100644
--- a/profiles_eval/configs/dlwindstat.json
+++ b/profiles_eval/configs/dlwindstat.json
@@ -15,6 +15,10 @@
       "name_in_file": "height",
       "units": "m"
     },
+    "w_variance": {
+      "name_in_file": "w_variance",
+      "units": "m^2/s^2"
+    },
     "cloud_base": {
       "name_in_file": "dl_cbh",
       "units": "m"
diff --git a/profiles_eval/real_lidar_corrections.py b/profiles_eval/real_lidar_corrections.py
index bd4ae55..813894a 100644
--- a/profiles_eval/real_lidar_corrections.py
+++ b/profiles_eval/real_lidar_corrections.py
@@ -21,11 +21,18 @@ def apply_deadtime_correction(
     raw_signal: np.ndarray,
     correction_counts: np.ndarray,
     correction_factors: np.ndarray,
+    poly_degree: int = 1,
+    n_extrap_samples: int = 3,
 ) -> np.ndarray:
     """
     Apply deadtime correction to raw photon-counting lidar signals using
     a precomputed lookup table (LUT).
-    Out-of-range bins are extrapolated using the nearest boundary factor.
+
+    Values within the LUT range are linearly interpolated.  Values above the
+    maximum LUT count rate are extrapolated via a polynomial fit in log-log
+    space fitted to the last ``n_extrap_samples`` LUT points, or clamped
+    if ``poly_degree=0``.  Values below the minimum LUT count rate are
+    clamped to the first LUT factor.
 
     Parameters
     ----------
@@ -41,6 +48,13 @@ def apply_deadtime_correction(
         1D array of shape (n_lut,) containing the multiplicative correction
         factor at each count rate. A value of 1.0 means no correction;
         values > 1.0 indicate the true signal exceeds the measured signal.
+    poly_degree : int, optional
+        Polynomial degree for log-log extrapolation above the LUT maximum.
+        Default is 1 (power-law). Set to 0 to clamp to the last LUT factor.
+
+    n_extrap_samples : int, optional
+        Number of trailing LUT points used to fit the extrapolation polynomial.
+        Ignored when ``poly_degree=0``. Default is 3.
 
     Returns
     -------
@@ -48,13 +62,24 @@ def apply_deadtime_correction(
         2D array of shape (n_profiles, n_range_bins) with deadtime correction
         applied, in counts/us.
     """
-    interp_factors = np.interp(
-        raw_signal.ravel(),
-        correction_counts,
-        correction_factors,
-    ).reshape(raw_signal.shape)
+    if not isinstance(poly_degree, int) or poly_degree < 0:
+        raise ValueError(f"poly_degree must be a non-negative integer, got {poly_degree}")
+
+    flat = raw_signal.ravel()
+    interp_factors = np.interp(flat, correction_counts, correction_factors)
+
+    if poly_degree > 0:
+        upper_mask = flat > correction_counts[-1]
+        if upper_mask.any():
+            n = min(n_extrap_samples, len(correction_counts))
+            log_x_fit = np.log(correction_counts[-n:])
+            log_y_fit = np.log(correction_factors[-n:])
+            coeffs = np.polyfit(log_x_fit, log_y_fit, poly_degree)
+            interp_factors[upper_mask] = np.exp(
+                np.polyval(coeffs, np.log(flat[upper_mask]))
+            )
 
-    return raw_signal * interp_factors
+    return raw_signal * interp_factors.reshape(raw_signal.shape)
 
 
 def apply_afterpulse_correction(
@@ -212,6 +237,8 @@ def compute_nrb(
     overlap_range: np.ndarray,
     overlap_factors: np.ndarray,
     calibration_constant: float = 1.0,
+    deadtime_poly_degree: int = 1,
+    deadtime_n_extrap_samples: int = 3,
 ) -> np.ndarray:
     """
     Compute the Normalized Relative Backscatter (NRB) from raw lidar data.
@@ -263,6 +290,13 @@ def compute_nrb(
     calibration_constant : float, optional
         Instrument calibration constant C. Use 1.0 for relative NRB.
         Default is 1.0.
+    deadtime_poly_degree : int, optional
+        Polynomial degree for deadtime correction log-log extrapolation above
+        the LUT maximum. Default is 1 (power-law). Set to 0 to clamp to the
+        last LUT factor (nearest-neighbour behaviour).
+    deadtime_n_extrap_samples : int, optional
+        Number of trailing LUT points used to fit the deadtime extrapolation
+        polynomial. Ignored when ``deadtime_poly_degree=0``. Default is 3.
 
     Returns
     -------
@@ -275,6 +309,8 @@ def compute_nrb(
         raw_signal,
         deadtime_correction_counts,
         deadtime_correction_factors,
+        poly_degree=deadtime_poly_degree,
+        n_extrap_samples=deadtime_n_extrap_samples,
     )
 
     # Step 2 — afterpulse correction
@@ -313,6 +349,8 @@ def compute_nrb_dataset(
     calibration_constant: float = 1.0,
     cross_pol: bool = True,
     config_dir: str = "./configs",
+    deadtime_poly_degree: int = 1,
+    deadtime_n_extrap_samples: int = 3,
 ) -> "xr.Dataset":
     """
     Compute co-pol NRB and, optionally, cross-pol NRB and linear depolarization
@@ -353,6 +391,13 @@ def compute_nrb_dataset(
     config_dir : str, optional
         Directory containing the JSON configuration files. Only used when
         ``instrument_type`` is provided. Default is ``"./configs"``.
+    deadtime_poly_degree : int, optional
+        Polynomial degree for deadtime correction log-log extrapolation above
+        the LUT maximum. Default is 1 (power-law). Set to 0 to clamp to the
+        last LUT factor (nearest-neighbour behaviour).
+    deadtime_n_extrap_samples : int, optional
+        Number of trailing LUT points used to fit the deadtime extrapolation
+        polynomial. Ignored when ``deadtime_poly_degree=0``. Default is 3.
 
     Returns
     -------
@@ -387,7 +432,18 @@ def _units(entry):
 
     # Extract all arrays from the dataset upfront
     range_km         = ds[_name(v["range"])].values
-    raw_co           = ds[_name(v["raw_signal_co_pol"])].values
+    
+    # Extract raw co-pol signal and ensure (time, range_or_spatial_dim) layout
+    raw_co_da = ds[_name(v["raw_signal_co_pol"])]
+    # Find the spatial dimension name (could be 'range', 'range_bins', 'height', etc.)
+    dims_list = list(raw_co_da.dims)
+    time_dim = next((d for d in dims_list if d == "time"), None)
+    spatial_dim = next((d for d in dims_list if d != "time"), None)
+    
+    if time_dim and spatial_dim and [time_dim, spatial_dim] != dims_list:
+        raw_co_da = raw_co_da.transpose(time_dim, spatial_dim)
+    raw_co = raw_co_da.values
+    
     background_co    = ds[_name(v["background_co_pol"])].values
     dt_counts        = ds[c["deadtime_counts"]].values
     dt_factors       = ds[c["deadtime_factors"]].values
@@ -411,9 +467,11 @@ def _units(entry):
         overlap_range               = overlap_range,
         overlap_factors             = overlap_factors,
         calibration_constant        = calibration_constant,
+        deadtime_poly_degree        = deadtime_poly_degree,
+        deadtime_n_extrap_samples   = deadtime_n_extrap_samples,
     )
 
-    dims = ("time", "range_bins")
+    dims = ("time", "range")
     data_vars = {
         _name(v["attenuated_backscatter"]): (
             dims, nrb_co,
@@ -422,7 +480,16 @@ def _units(entry):
     }
 
     if cross_pol:
-        raw_cross        = ds[_name(v["raw_signal_cross_pol"])].values
+        # Extract raw cross-pol signal and ensure (time, range_or_spatial_dim) layout
+        raw_cross_da = ds[_name(v["raw_signal_cross_pol"])]
+        dims_list = list(raw_cross_da.dims)
+        time_dim = next((d for d in dims_list if d == "time"), None)
+        spatial_dim = next((d for d in dims_list if d != "time"), None)
+        
+        if time_dim and spatial_dim and [time_dim, spatial_dim] != dims_list:
+            raw_cross_da = raw_cross_da.transpose(time_dim, spatial_dim)
+        raw_cross = raw_cross_da.values
+        
         background_cross = ds[_name(v["background_cross_pol"])].values
         ap_profile_cross = ds[c["afterpulse_profile_cross_pol"]].values
         dc_profile_cross = ds[c["darkcounts_profile_cross_pol"]].values
@@ -440,6 +507,8 @@ def _units(entry):
             overlap_range               = overlap_range,
             overlap_factors             = overlap_factors,
             calibration_constant        = calibration_constant,
+            deadtime_poly_degree        = deadtime_poly_degree,
+            deadtime_n_extrap_samples   = deadtime_n_extrap_samples,
         )
         ldr = nrb_cross / (nrb_co + nrb_cross)
         data_vars[_name(v["attenuated_backscatter_cross_pol"])] = (
@@ -455,8 +524,8 @@ def _units(entry):
     return xr.Dataset(
         data_vars,
         coords={
-            "time":       ds[_name(v["time"])],
-            "range_bins": range_km,
+            "time":  ds[_name(v["time"])],
+            "range": range_km,
         },
     )
 
@@ -472,23 +541,27 @@ def _units(entry):
     # Load data and compute NRB dataset
     # ------------------------------------------------------------------
     FILE = "/data/archive/sgp/sgpminimplC1.b1/sgpminimplC1.b1.20260214.000009.nc"
+    FILE = "/data/archive/sgp/sgpminimplC1.b1/sgpminimplC1.b1.20260612.000004.nc"  # test case for Donna
+    FILE = "/data/archive/kcg/kcgmplpolfsM1.b1/kcgmplpolfsM1.b1.20240601.000009.nc"
+    FILE = "/data/archive/kcg/kcgminimplS1.b1/kcgminimplS1.b1.20240601.000000.nc"  # test case from Damao
     ds     = xr.open_dataset(FILE)
-    result = compute_nrb_dataset(ds, instrument_type="minimpl", config_dir="./configs")
+    result = compute_nrb_dataset(ds, instrument_type="minimpl", config_dir="./configs",
+                                 deadtime_poly_degree=1, deadtime_n_extrap_samples=3)
 
     has_cross    = "ldr" in result
-    max_range_km = 10.0                         # maximum y-axis range (km)
+    max_range_km = 5.0                         # maximum y-axis range (km)
 
     # Build plot dataset — keep only up to max_range_km
-    range_sel = result.range_bins <= max_range_km
-    plot_ds = result.sel(range_bins=range_sel)
+    range_sel = result.range <= max_range_km
+    plot_ds = result.sel(range=range_sel)
 
     # Add log10-transformed fields (raw signal and NRB co); LDR stays linear
     range_mask = ds["range"].values <= max_range_km
     raw_vals   = ds["signal_return_co_pol"].values[:, range_mask]
     plot_ds["raw"] = xr.DataArray(
         np.log10(np.where(raw_vals > 0, raw_vals, np.nan)),
-        dims=["time", "range_bins"],
-        coords={"time": plot_ds.time, "range_bins": plot_ds.range_bins},
+        dims=["time", "range"],
+        coords={"time": plot_ds.time, "range": plot_ds.range},
     )
     plot_ds["log_nrb_co"] = np.log10(plot_ds["nrb_co"].where(plot_ds["nrb_co"] > 0))
 
@@ -513,12 +586,16 @@ def _units(entry):
         vmax = fixed_vmax if fixed_vmax is not None else float(da.quantile(0.99))
 
         da.plot.pcolormesh(
-            ax=ax_c, x="time", y="range_bins",
+            ax=ax_c, x="time", y="range",
             cmap=cmap, vmin=vmin, vmax=vmax,
+            #cmap=cmap, vmin=0.0, vmax=0.08,
             cbar_kwargs={"label": cb_label},
         )
         ax_c.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
-        ax_c.set_ylim(float(plot_ds.range_bins[0]), max_range_km)
+        t_base = plot_ds.time.values[0].astype('datetime64[D]')
+        #ax_c.set_xlim(t_base + np.timedelta64(14, 'h'), t_base + np.timedelta64(17, 'h'))
+        ax_c.set_ylim(float(plot_ds.range[0]), max_range_km)
+        #ax_c.set_ylim((0.5, 1.5))
         ax_c.set_title(f"{title} Curtain")
         ax_c.set_xlabel("Time (UTC)")
         ax_c.set_ylabel("Range (km)")
diff --git a/profiles_eval/real_prof_io.py b/profiles_eval/real_prof_io.py
index 21fa5a5..d91628b 100644
--- a/profiles_eval/real_prof_io.py
+++ b/profiles_eval/real_prof_io.py
@@ -425,6 +425,8 @@ def load_and_process_arm_data(
     range_km: "np.ndarray | None" = None,
     time_step: np.timedelta64 = np.timedelta64(15, "s"),
     data_path_template: str | None = None,
+    deadtime_poly_degree: int = 1,
+    deadtime_n_extrap_samples: int = 3,
 ) -> xr.Dataset:
     """
     Load, optionally correct, and optionally interpolate ARM data.
@@ -483,6 +485,13 @@ def load_and_process_arm_data(
         :func:`find_arm_files`).  When provided, ``data_path`` may be
         ``None`` and the concrete path is derived per-instrument.  This is
         especially convenient when ``instrument_type`` is a list.
+    deadtime_poly_degree : int, optional
+        Polynomial degree for deadtime correction log-log extrapolation above
+        the LUT maximum. Default is 1 (power-law). Set to 0 to clamp to the
+        last LUT factor (nearest-neighbour behaviour).
+    deadtime_n_extrap_samples : int, optional
+        Number of trailing LUT points used to fit the deadtime extrapolation
+        polynomial. Ignored when ``deadtime_poly_degree=0``. Default is 3.
 
     Returns
     -------
@@ -504,6 +513,8 @@ def load_and_process_arm_data(
                     time_min, time_max, safety_delta, config_dir,
                     interpolate, range_km, time_step,
                     data_path_template=data_path_template,
+                    deadtime_poly_degree=deadtime_poly_degree,
+                    deadtime_n_extrap_samples=deadtime_n_extrap_samples,
                 )
             except FileNotFoundError as exc:
                 warnings.warn(f"[{instr}] skipped — {exc}", stacklevel=2)
@@ -652,7 +663,9 @@ def load_and_process_arm_data(
                        else v.get("raw_signal_cross_pol", ""))
 
         ds_corrected = compute_nrb_dataset(
-            ds, instrument_type=instrument_type, config_dir=config_dir
+            ds, instrument_type=instrument_type, config_dir=config_dir,
+            deadtime_poly_degree=deadtime_poly_degree,
+            deadtime_n_extrap_samples=deadtime_n_extrap_samples,
         )
 
         # Propagate source attributes; read from any variable in ds since
diff --git a/profiles_eval/real_prof_main.py b/profiles_eval/real_prof_main.py
index 744260a..4d27462 100644
--- a/profiles_eval/real_prof_main.py
+++ b/profiles_eval/real_prof_main.py
@@ -47,6 +47,7 @@
     "dlwindstat",
     "dlwind",
     "interpolatedsonde",
+    "pblhtbeml",
 ]
 
 # ---------------------------------------------------------------------------
@@ -226,6 +227,12 @@ def _instr_group_key(var):
         default=None,
         help=f"Data path template (default: {DATA_PATH_TEMPLATE})",
     )
+    parser.add_argument(
+        "--output-path",
+        type=str,
+        default="./",
+        help="Output directory for exported files (default: ./)",
+    )
     args = parser.parse_args()
     
     # Parse period_start and period_end from CLI or use defaults
@@ -260,6 +267,11 @@ def _instr_group_key(var):
             # Use resolved site and facility values
             resolved_site = args.site if args.site else SITE
             resolved_facility = args.facility if args.facility else FACILITY
-            out = export_dataset(result, site=resolved_site, facility=resolved_facility)
+            out = export_dataset(
+                result,
+                site=resolved_site,
+                facility=resolved_facility,
+                output_path=args.output_path,
+            )
             print(f"\nExported: {out}")
         t = t_next
diff --git a/profiles_eval/real_prof_plot.py b/profiles_eval/real_prof_plot.py
index ec50e58..f11bcd6 100644
--- a/profiles_eval/real_prof_plot.py
+++ b/profiles_eval/real_prof_plot.py
@@ -8,6 +8,8 @@
 from __future__ import annotations
 
 import warnings
+from collections import OrderedDict
+from pathlib import Path
 
 import matplotlib.colors as mcolors
 import matplotlib.dates as mdates
@@ -25,10 +27,47 @@
 # ---------------------------------------------------------------------------
 
 def _is_log_var(varname: str) -> bool:
-    """Return True if the variable should be rendered with a log colour scale."""
+    """Return True if the variable should be rendered with a log color scale."""
     return "backscatter" in varname or "extinction" in varname
 
 
+def _extract_instr_field(varname: str) -> tuple[str, str]:
+    """Extract (instrument, field_name) from a DataArray variable name.
+    
+    Removes instrument prefix and optional ``_supp_`` marker.
+    
+    Parameters
+    ----------
+    varname : str
+        Full variable name, e.g., ``"mpl_supp_particulate_backscatter"`` or
+        ``"hsrl_molecular_signal_to_noise"``.
+    
+    Returns
+    -------
+    instrument : str
+        Instrument code (e.g., ``"mpl"``, ``"hsrl"``).
+    field : str
+        Field name with instrument prefix and ``_supp_`` removed
+        (e.g., ``"particulate_backscatter"``).
+    
+    Examples
+    --------
+    >>> _extract_instr_field("mpl_supp_particulate_backscatter")
+    ('mpl', 'particulate_backscatter')
+    >>> _extract_instr_field("hsrl_molecular_signal_to_noise")
+    ('hsrl', 'molecular_signal_to_noise')
+    """
+    # Temporarily remove _supp_ for parsing
+    temp_name = varname.replace("_supp_", "_", 1)
+    
+    # Split on first underscore
+    parts = temp_name.split("_", 1)
+    if len(parts) == 2:
+        return parts[0], parts[1]
+    # Fallback if no underscore found
+    return varname, varname
+
+
 def _safe_norm(
     data: np.ndarray,
     log: bool,
@@ -477,6 +516,273 @@ def _finite_ylim(da):
     return fig, axes
 
 
+# ---------------------------------------------------------------------------
+# Multi-instrument comparison plots
+# ---------------------------------------------------------------------------
+
+def plot_profile_curtains(
+    variables: "OrderedDict[str, xr.DataArray]",
+    cmap: str = "viridis",
+    shared_norm: bool = True,
+    ylim: tuple[float, float] | None = None,
+    fig_width: float = 10.0,
+    panel_height: float = 3.0,
+    output_path: "str | Path | None" = None,
+    **kwargs,
+) -> tuple[plt.Figure, np.ndarray]:
+    """Plot curtain panels for multiple instruments from a comparison dict.
+
+    One panel is created per entry in *variables* (HSRL first by convention).
+    Backscatter and extinction variables are rendered with a shared
+    logarithmic color scale; other products use a linear scale.
+
+    Parameters
+    ----------
+    variables : OrderedDict[str, xr.DataArray]
+        ``{legend_label: DataArray}`` as returned by
+        :func:`real_prof_analysis_utils.get_comparison_variables`.
+    cmap : str, optional
+        Colormap name.  Default ``"viridis"``.
+    shared_norm : bool, optional
+        When ``True`` (default), all panels share the same color scale,
+        making inter-instrument differences immediately visible.
+    ylim : tuple of float, optional
+        ``(ymin, ymax)`` range axis limits in km.
+    fig_width : float, optional
+        Figure width in inches.  Default ``10``.
+    panel_height : float, optional
+        Height per panel in inches.  Default ``3``.
+    output_path : str or Path, optional
+        If provided, the figure is saved to this path.
+    **kwargs
+        Forwarded to :func:`matplotlib.pyplot.subplots`.
+
+    Returns
+    -------
+    fig : matplotlib.figure.Figure
+    axes : np.ndarray of matplotlib.axes.Axes
+    """
+    from collections import OrderedDict as _OD
+
+    labels = list(variables.keys())
+    arrays = list(variables.values())
+    n = len(arrays)
+
+    # Determine whether log scale applies from the first variable name
+    first_name = arrays[0].name or ""
+    log = _is_log_var(first_name)
+
+    # Build a shared normalisation across all panels if requested
+    norm_shared: mcolors.Normalize | None = None
+    if shared_norm:
+        all_data = np.concatenate([da.values.ravel() for da in arrays])
+        norm_shared = _safe_norm(all_data, log, varname=first_name)
+
+    fig, _axes = plt.subplots(
+        n, 1, figsize=(fig_width, panel_height * n), squeeze=False, **kwargs
+    )
+    axes: np.ndarray = _axes.ravel()
+
+    for ax, label, da in zip(axes, labels, arrays):
+        panel_norm = norm_shared if shared_norm else _safe_norm(
+            da.values.ravel(), log, varname=da.name or ""
+        )
+        # Build panel title: INSTRUMENT - field_name (with underscores as spaces)
+        instr, field = _extract_instr_field(da.name or "")
+        title = f"{instr.upper()} - {field.replace('_', ' ')}"
+        mesh = _plot_curtain(ax, da, norm=panel_norm, cmap=cmap, title=title)
+        plt.colorbar(mesh, ax=ax, pad=0.02, label=da.attrs.get("units", ""))
+        if ylim is not None:
+            ax.set_ylim(ylim)
+
+    fig.tight_layout()
+
+    if output_path is not None:
+        fig.savefig(output_path, dpi=150, bbox_inches="tight")
+
+    return fig, axes
+
+
+def plot_time_mean_profiles(
+    variables: "OrderedDict[str, xr.DataArray]",
+    time_slice: "slice | None" = None,
+    mask: "xr.DataArray | None" = None,
+    ax: "plt.Axes | None" = None,
+    ylim: "tuple[float, float] | None" = None,
+    figsize: "tuple[float, float]" = (5.0, 8.0),
+    output_path: "str | Path | None" = None,
+) -> tuple[plt.Figure, plt.Axes]:
+    """Plot time-averaged profiles for multiple instruments on a single axes.
+
+    Each instrument is drawn as a separate line.  Backscatter variables use
+    a logarithmic x-axis; SNR and LDR use a linear x-axis.
+
+    Parameters
+    ----------
+    variables : OrderedDict[str, xr.DataArray]
+        ``{legend_label: DataArray}`` as returned by
+        :func:`real_prof_analysis_utils.get_comparison_variables`.
+    time_slice : slice, optional
+        Restrict the time average to a subset, e.g.
+        ``slice("2026-03-10T20:00", "2026-03-10T21:00")``.
+    mask : xr.DataArray, optional
+        Boolean mask on ``(time, range)``.  Masked-out values are excluded
+        from the mean.  Typically the output of
+        :func:`real_prof_analysis_utils.build_lidar_data_mask`.
+    ax : matplotlib.axes.Axes, optional
+        Axes to draw on.  A new figure is created when ``None`` (default).
+    ylim : tuple of float, optional
+        ``(ymin, ymax)`` range axis limits in km.
+    figsize : tuple of float, optional
+        ``(width, height)`` in inches when creating a new figure.
+        Default ``(5, 8)``.
+    output_path : str or Path, optional
+        If provided, the figure is saved to this path.
+
+    Returns
+    -------
+    fig : matplotlib.figure.Figure
+    ax : matplotlib.axes.Axes
+    """
+    if ax is None:
+        fig, ax = plt.subplots(figsize=figsize)
+    else:
+        fig = ax.get_figure()
+
+    first_name = next(iter(variables.values())).name or ""
+    log_x = _is_log_var(first_name)
+
+    for label, da in variables.items():
+        # Optionally restrict to a time window
+        if time_slice is not None:
+            da = da.sel(time=time_slice)
+
+        # Ensure (time, range) layout before averaging
+        if list(da.dims) != ["time", "range"]:
+            da = da.transpose("time", "range")
+
+        # Apply mask: replace masked values with NaN before averaging
+        if mask is not None:
+            m = mask.sel(time=da.time) if "time" in mask.dims else mask
+            da = da.where(m)
+
+        # Time-mean: nanmean along the time axis (axis=0 → time)
+        profile = np.nanmean(da.values, axis=0)
+        range_coord = da["range"].values
+
+        ax.plot(profile, range_coord, label=label, linewidth=1.2)
+
+    ax.set_ylabel("Range (km)")
+    ax.set_xlabel(f"{first_name}  [{next(iter(variables.values())).attrs.get('units', '')}]")
+    ax.legend(fontsize=8, loc="upper right")
+    ax.grid(True, linestyle="--", alpha=0.4)
+
+    if log_x:
+        ax.set_xscale("log")
+    if ylim is not None:
+        ax.set_ylim(ylim)
+
+    fig.tight_layout()
+
+    if output_path is not None:
+        fig.savefig(output_path, dpi=150, bbox_inches="tight")
+
+    return fig, ax
+
+
+def plot_cfad(
+    cfad_data_dict: "dict[str, tuple[np.ndarray, np.ndarray, np.ndarray]]",
+    product: str = "",
+    cmap: str = "Blues",
+    ylim: "tuple[float, float] | None" = None,
+    fig_width: float = 4.5,
+    panel_height: float = 6.0,
+    output_path: "str | Path | None" = None,
+    **kwargs,
+) -> tuple[plt.Figure, np.ndarray]:
+    """Plot Contoured Frequency by Altitude Diagrams (CFADs) for multiple instruments.
+
+    One subplot is created per instrument.  Each panel shows relative
+    frequency as a color fill, with range on the y-axis and value on the
+    x-axis.  Backscatter/extinction products use a logarithmic x-axis.
+
+    Parameters
+    ----------
+    cfad_data_dict : dict[str, (freq_2d, range_centers, value_centers)]
+        Pre-computed CFAD data per label, as returned by
+        :func:`real_prof_analysis_utils.compute_cfad_data`.
+    product : str, optional
+        Product name used solely for the x-axis label.
+    cmap : str, optional
+        Colormap for the 2-D frequency fill.  Default ``"Blues"``.
+    ylim : tuple of float, optional
+        ``(ymin, ymax)`` range axis limits in km.
+    fig_width : float, optional
+        Width of each subplot panel in inches.  Default ``4.5``.
+    panel_height : float, optional
+        Height of each subplot panel in inches.  Default ``6``.
+    output_path : str or Path, optional
+        If provided, the figure is saved to this path.
+    **kwargs
+        Forwarded to :func:`matplotlib.pyplot.subplots`.
+
+    Returns
+    -------
+    fig : matplotlib.figure.Figure
+    axes : np.ndarray of matplotlib.axes.Axes
+    """
+    labels = list(cfad_data_dict.keys())
+    n = len(labels)
+    log_x = _is_log_var(product)
+
+    fig, _axes = plt.subplots(
+        1, n, figsize=(fig_width * n, panel_height),
+        squeeze=False, sharey=True, **kwargs,
+    )
+    axes: np.ndarray = _axes.ravel()
+
+    for ax, label in zip(axes, labels):
+        freq_2d, range_centers, value_centers = cfad_data_dict[label]
+
+        # pcolormesh expects edges, not centers — compute from center spacing
+        def _edges(centers: np.ndarray) -> np.ndarray:
+            half = np.diff(centers) / 2
+            return np.concatenate([
+                [centers[0] - half[0]],
+                centers[:-1] + half,
+                [centers[-1] + half[-1]],
+            ])
+
+        v_edges = _edges(value_centers)
+        r_edges = _edges(range_centers)
+
+        # Mask zero-frequency cells so they render as white (not color)
+        data = np.ma.masked_where(~np.isfinite(freq_2d) | (freq_2d == 0), freq_2d)
+
+        mesh = ax.pcolormesh(
+            v_edges, r_edges, data,
+            cmap=cmap, vmin=0, shading="flat",
+        )
+        plt.colorbar(mesh, ax=ax, pad=0.02, label="Relative frequency")
+
+        ax.set_title(label, fontsize=9)
+        ax.set_xlabel(product)
+        ax.set_ylabel("Range (km)")
+        ax.grid(True, linestyle="--", alpha=0.3)
+
+        if log_x and value_centers[value_centers > 0].size > 0:
+            ax.set_xscale("log")
+        if ylim is not None:
+            ax.set_ylim(ylim)
+
+    fig.tight_layout()
+
+    if output_path is not None:
+        fig.savefig(output_path, dpi=150, bbox_inches="tight")
+
+    return fig, axes
+
+
 # ---------------------------------------------------------------------------
 # Example / quick-look
 # ---------------------------------------------------------------------------
diff --git a/profiles_eval/real_prof_utils.py b/profiles_eval/real_prof_utils.py
index 6cbc052..9252e99 100644
--- a/profiles_eval/real_prof_utils.py
+++ b/profiles_eval/real_prof_utils.py
@@ -114,7 +114,7 @@ def interpolate_data(
     step_ns = int(time_step / np.timedelta64(1, "ns"))
     t_min_ns = time_min.astype("datetime64[ns]").astype(np.int64)
     t_max_ns = time_max.astype("datetime64[ns]").astype(np.int64)
-    out_time = (np.arange(t_min_ns, t_max_ns + step_ns, step_ns)
+    out_time = (np.arange(t_min_ns, t_max_ns, step_ns)
                   .astype("datetime64[ns]"))
 
     # ------------------------------------------------------------------