diff --git a/cloudnetpy/cli.py b/cloudnetpy/cli.py
index 2f6fab50..db40214e 100644
--- a/cloudnetpy/cli.py
+++ b/cloudnetpy/cli.py
@@ -20,7 +20,7 @@
 
 from cloudnetpy import concat_lib, instruments
 from cloudnetpy.categorize import CategorizeInput, generate_categorize
-from cloudnetpy.exceptions import PlottingError
+from cloudnetpy.exceptions import ModelDataError, PlottingError
 from cloudnetpy.plotting import PlotParameters, generate_figure
 
 if TYPE_CHECKING:
@@ -29,6 +29,17 @@
 
 cloudnet_api_url: Final = "https://cloudnet.fmi.fi/api/"
 
+# Model-evaluation L3 products and the observation product they downsample.
+# The API does not declare these source products, so they are mapped here.
+L3_SOURCE_PRODUCTS: Final = {
+    "l3-cf": "categorize",
+    "l3-iwc": "iwc",
+    "l3-lwc": "lwc",
+}
+
+# Model used for L3 products when none is given with --model.
+DEFAULT_L3_MODEL: Final = "ecmwf"
+
 
 def run(args: argparse.Namespace, tmpdir: str, client: APIClient) -> None:
     cat_files = {}
@@ -96,6 +107,16 @@ def run(args: argparse.Namespace, tmpdir: str, client: APIClient) -> None:
         if epsilon_filepath is not None:
             _plot(epsilon_filepath, "epsilon-radar", args)
 
+    # Model evaluation L3 products (e.g. l3-cf)
+    for product in args.products:
+        base_product = _parse_instrument(product)[0]
+        if base_product in L3_SOURCE_PRODUCTS:
+            model = args.model or DEFAULT_L3_MODEL
+            if args.model is None:
+                logging.info("No model specified, using default '%s'", model)
+            l3_filepath = _process_l3_product(base_product, args, client, model)
+            _plot_l3(l3_filepath, base_product, args)
+
 
 def _process_epsilon_radar(
     cat_files: dict, args: argparse.Namespace, client: APIClient
@@ -118,6 +139,47 @@ def _process_epsilon_radar(
     return str(output_file)
 
 
+def _process_l3_product(
+    product: str, args: argparse.Namespace, client: APIClient, model: str
+) -> str | None:
+    obs = product.removeprefix("l3-")
+    source_product = L3_SOURCE_PRODUCTS[product]
+    product_file = _fetch_product(args, source_product, client)
+    if product_file is None:
+        logging.info("No %s data available for %s", source_product, product)
+        return None
+    model_file = _fetch_model(args, client, model)
+    if model_file is None:
+        logging.info("No model data available for %s", product)
+        return None
+    filename = f"{args.date.replace('-', '')}_{args.site}_{model}_{product}.nc"
+    output_file = _create_output_folder("evaluation", args) / filename
+    model_name = next((m.name for m in client.models() if m.id == model), model)
+    site_name = next(
+        (s.human_readable_name for s in client.sites() if s.id == args.site),
+        args.site,
+    )
+    module = importlib.import_module(
+        "cloudnetpy.model_evaluation.products.product_resampling"
+    )
+    try:
+        module.process_L3_day_product(
+            model,
+            obs,
+            model_file,
+            product_file,
+            str(output_file),
+            model_name=model_name,
+            site_name=site_name,
+            overwrite=True,
+        )
+    except ModelDataError as e:
+        logging.info("Failed to process %s: %s", product, e)
+        return None
+    logging.info("Processed %s: %s", product, output_file)
+    return str(output_file)
+
+
 def _process_categorize(
     input_files: dict,
     instrument_prefs: dict[str, list[str]],
@@ -536,9 +598,14 @@ def _fetch_product(
     return _download_product_file(meta, folder, client, force=args.force_download)
 
 
-def _fetch_model(args: argparse.Namespace, client: APIClient) -> str | None:
+def _fetch_model(
+    args: argparse.Namespace, client: APIClient, model_id: str | None = None
+) -> str | None:
     files = client.files(
-        product_id="model", model_id=args.model, date=args.date, site_id=args.site
+        product_id="model",
+        model_id=model_id or args.model,
+        date=args.date,
+        site_id=args.site,
     )
     if not files:
         logging.info("No model data available for this date")
@@ -623,6 +690,33 @@ def _plot(
         logging.info("Plotted %s: %s", product, image_name)
 
 
+def _plot_l3(
+    filepath: PathLike | str | None,
+    product: str,
+    args: argparse.Namespace,
+) -> None:
+    if filepath is None or (not args.plot and not args.show):
+        return
+    obs = product.removeprefix("l3-")
+    save_path = f"{Path(filepath).parent}/" if args.plot else None
+    var_list = args.variables.split(",") if args.variables is not None else None
+    module = importlib.import_module("cloudnetpy.model_evaluation.plotting.plotting")
+    try:
+        module.generate_L3_day_plots(
+            str(filepath),
+            obs,
+            var_list=var_list,
+            save_path=save_path,
+            show=args.show,
+            include_advection=False,
+        )
+    except (PlottingError, ValueError, KeyError) as e:
+        logging.info("Failed to plot %s: %s", product, e)
+        return
+    if args.plot:
+        logging.info("Plotted %s to %s", product, save_path)
+
+
 def _process_cat_product(product: str, categorize_file: str) -> str:
     output_file = categorize_file.replace("categorize", product)
     module = importlib.import_module("cloudnetpy.products")
@@ -648,7 +742,7 @@ def _parse_products(product_argument: str, client: APIClient) -> list[str]:
     valid_products = []
     for product in products:
         prod, _ = _parse_instrument(product)
-        if prod in valid_options:
+        if prod in valid_options or prod in L3_SOURCE_PRODUCTS:
             valid_products.append(product)
     return valid_products
 
@@ -704,7 +798,7 @@ def main() -> None:
         "-m",
         "--model",
         type=lambda arg: _parse_model(arg, client),
-        help="Model to use in categorize.",
+        help="Model to use in categorize and model evaluation (l3-*) products.",
     )
     parser.add_argument(
         "-i",
diff --git a/cloudnetpy/model_evaluation/file_handler.py b/cloudnetpy/model_evaluation/file_handler.py
index 92fdf638..debf5e50 100644
--- a/cloudnetpy/model_evaluation/file_handler.py
+++ b/cloudnetpy/model_evaluation/file_handler.py
@@ -1,12 +1,13 @@
 import os
 from datetime import datetime
 from os import PathLike
+from typing import TYPE_CHECKING
 from uuid import UUID
 
 import netCDF4
-import numpy.typing as npt
 
-from cloudnetpy import output, utils
+from cloudnetpy import output
+from cloudnetpy.model_evaluation.model_metadata import MODEL_PREFIX, PRODUCT_NAMES
 
 from .metadata import (
     CYCLE_ATTRIBUTES,
@@ -14,147 +15,101 @@
     MODEL_L3_ATTRIBUTES,
     REGRID_PRODUCT_ATTRIBUTES,
 )
-from .products.model_products import ModelManager
+
+if TYPE_CHECKING:
+    from cloudnetpy.model_evaluation.products.model_products import ModelManager
+    from cloudnetpy.model_evaluation.products.observation_products import (
+        ObservationManager,
+    )
 
 
 def update_attributes(model_downsample_variables: dict, attributes: dict) -> None:
-    """Overrides existing Cloudnet-ME Array-attributes.
-        Overrides existing attributes using hard-coded values.
-        New attributes are added.
+    """Sets variable attributes for the L3 downsampled file.
+
+    Model (simulated) fields are prefixed with ``model_``; observation fields
+    downsampled to the model grid use the bare product key.
 
     Args:
         model_downsample_variables (dict): Array instances.
-        attributes (dict): Product-specific attributes.
+        attributes (dict): Product-specific attributes (e.g. time units).
     """
-    for key in model_downsample_variables:
-        x = len(key.split("_")) - 1
-        key_parts = key.split("_", x)
-        if key in list(attributes.keys()):
-            model_downsample_variables[key].set_attributes(attributes[key])
+    for key, variable in model_downsample_variables.items():
+        if key in attributes:
+            variable.set_attributes(attributes[key])
         if key in MODEL_ATTRIBUTES:
-            model_downsample_variables[key].set_attributes(MODEL_ATTRIBUTES[key])
-        elif "_".join(key_parts[0:-1]) in REGRID_PRODUCT_ATTRIBUTES:
-            model_downsample_variables[key].set_attributes(
-                REGRID_PRODUCT_ATTRIBUTES["_".join(key_parts[0:-1])],
-            )
-        elif "_".join(key_parts[0:-2]) in REGRID_PRODUCT_ATTRIBUTES:
-            model_downsample_variables[key].set_attributes(
-                REGRID_PRODUCT_ATTRIBUTES["_".join(key_parts[0:-2])],
-            )
-        elif (
-            "_".join(key_parts[1:]) in MODEL_L3_ATTRIBUTES
-            or "_".join(key_parts[2:]) in MODEL_L3_ATTRIBUTES
-        ):
-            try:
-                model_downsample_variables[key].set_attributes(
-                    MODEL_L3_ATTRIBUTES["_".join(key_parts[1:])],
-                )
-            except KeyError:
-                model_downsample_variables[key].set_attributes(
-                    MODEL_L3_ATTRIBUTES["_".join(key_parts[2:])],
-                )
-        elif "_".join(key_parts[1:]) in CYCLE_ATTRIBUTES:
-            model_downsample_variables[key].set_attributes(
-                CYCLE_ATTRIBUTES["_".join(key_parts[1:])],
-            )
-        elif "_".join(key_parts[2:]) in CYCLE_ATTRIBUTES:
-            model_downsample_variables[key].set_attributes(
-                CYCLE_ATTRIBUTES["_".join(key_parts[2:])],
-            )
+            variable.set_attributes(MODEL_ATTRIBUTES[key])
+        elif key.startswith(MODEL_PREFIX):
+            base = key.removeprefix(MODEL_PREFIX)
+            if base in MODEL_L3_ATTRIBUTES:
+                variable.set_attributes(MODEL_L3_ATTRIBUTES[base])
+            elif base in CYCLE_ATTRIBUTES:
+                variable.set_attributes(CYCLE_ATTRIBUTES[base])
+        elif key in REGRID_PRODUCT_ATTRIBUTES:
+            variable.set_attributes(REGRID_PRODUCT_ATTRIBUTES[key])
 
 
 def save_downsampled_file(
-    id_mark: str,
+    product: str,
     file_name: str | PathLike,
-    objects: tuple,
-    files: tuple[list[str | PathLike], str | PathLike],
+    model_obj: "ModelManager",
+    obs_obj: "ObservationManager",
+    model_files: list[str | PathLike],
+    product_file: str | PathLike,
     uuid: UUID,
+    model_name: str | None = None,
+    site_name: str | None = None,
 ) -> None:
     """Saves a standard downsampled day product file.
 
     Args:
-        id_mark (str): File identifier, format "(product name)_(model name)"
-        file_name (str): Name of the output file to be generated
-        objects (tuple): Include two objects: The :class:'ModelManager' and
-                      The :class:'ObservationManager.
-        files (tuple): Includes two sourcefile group: List of model file(s) used
-                       for processing output file and Cloudnet L2 product file
-        keep_uuid (bool): If True, keeps the UUID of the old file, if that exists.
-                          Default is False when new UUID is generated.
+        product (str): Product name, e.g. "cf".
+        file_name (str): Name of the output file to be generated.
+        model_obj (ModelManager): Model fields downsampled to the model grid.
+        obs_obj (ObservationManager): Cloudnet L2 observation product.
+        model_files (list): Model file(s) used for processing the output file.
+        product_file (str): Cloudnet L2 product file.
         uuid (str): Set specific UUID for the file.
+        model_name (str): Human-readable model name for plot titles. Falls back
+                       to the model id when not given.
+        site_name (str): Human-readable site name for the location attribute and
+                       plot subtitle. Falls back to the source file's location
+                       when not given.
     """
-    obj = objects[0]
-    dimensions = {"time": len(obj.time), "level": len(obj.data["level"][:])}
-    with output.init_file(file_name, dimensions, obj.data, uuid) as root_group:
+    n_levels = model_obj.data[model_obj.keys["height"]][:].shape[-1]
+    dimensions = {"time": len(model_obj.time), "level": n_levels}
+    location = site_name or model_obj.dataset.location
+    with output.init_file(file_name, dimensions, model_obj.data, uuid) as root_group:
         _augment_global_attributes(root_group)
-        root_group.cloudnet_file_type = "l3-" + id_mark.split("_", maxsplit=1)[0]
-        root_group.title = (
-            f"Downsampled {id_mark.capitalize().replace('_', ' of ')} "
-            f"from {obj.dataset.location}"
-        )
-        _add_source(root_group, objects, files)
-        output.copy_global(
-            obj.dataset, root_group, ("location", "day", "month", "year")
-        )
-        if not hasattr(obj.dataset, "day"):
-            root_group.year, root_group.month, root_group.day = obj.date
-        output.merge_history(root_group, id_mark, obj)
-
-
-def add_var2ncfile(obj: ModelManager, file_name: str | PathLike) -> None:
-    with netCDF4.Dataset(file_name, "r+", format="NETCDF4_CLASSIC") as nc_file:
-        _write_vars2nc(nc_file, obj.data)
-
-
-def _write_vars2nc(rootgrp: netCDF4.Dataset, cloudnet_variables: dict) -> None:
-    """Iterates over Cloudnet-ME instances and write to given rootgrp."""
-
-    def _get_dimensions(array: npt.NDArray) -> tuple:
-        """Finds correct dimensions for a variable."""
-        if utils.isscalar(array):
-            return ()
-        variable_size: tuple = ()
-        file_dims = rootgrp.dimensions
-        array_dims = array.shape
-        for length in array_dims:
-            dim = [key for key in file_dims if file_dims[key].size == length][0]  # noqa: RUF015
-            variable_size = (*variable_size, dim)
-        return variable_size
-
-    for key in cloudnet_variables:
-        obj = cloudnet_variables[key]
-        size = _get_dimensions(obj.data)
-        try:
-            nc_variable = rootgrp.createVariable(
-                obj.name,
-                obj.data_type,
-                size,
-                zlib=True,
-            )
-            nc_variable[:] = obj.data
-            for attr in obj.fetch_attributes():
-                setattr(nc_variable, attr, getattr(obj, attr))
-        except RuntimeError:
-            continue
+        root_group.cloudnet_file_type = "l3-" + product
+        product_name = PRODUCT_NAMES.get(product, product)
+        root_group.title = f"Observed and modeled {product_name} over {location}"
+        root_group.model_id = model_obj.model
+        root_group.model_name = model_name or model_obj.model
+        _add_source(root_group, model_obj, obs_obj, model_files, product_file)
+        output.copy_global(model_obj.dataset, root_group, ("day", "month", "year"))
+        root_group.location = location
+        if not hasattr(model_obj.dataset, "day"):
+            root_group.year, root_group.month, root_group.day = model_obj.date
+        output.merge_history(root_group, f"L3 {product_name}", model_obj)
 
 
 def _augment_global_attributes(root_group: netCDF4.Dataset) -> None:
     root_group.Conventions = "CF-1.8"
 
 
-def _add_source(root_ground: netCDF4.Dataset, objects: tuple, files: tuple) -> None:
+def _add_source(
+    root_group: netCDF4.Dataset,
+    model_obj: "ModelManager",
+    obs_obj: "ObservationManager",
+    model_files: list[str | PathLike],
+    product_file: str | PathLike,
+) -> None:
     """Generates source info for multiple files."""
-    model, obs = objects
-    model_files, obs_file = files
-    source = f"Observation file: {os.path.basename(obs_file)}"
-    source += "\n"
-    source += f"{model.model} file(s): "
-    for i, f in enumerate(model_files):
-        source += f"{os.path.basename(f)}"
-        if i < len(model_files) - 1:
-            source += "\n"
-    root_ground.source = source
-    root_ground.source_file_uuids = output.get_source_uuids([model, obs])
+    filenames = [os.path.basename(product_file)] + [
+        os.path.basename(f) for f in model_files
+    ]
+    root_group.source = "\n".join(filenames)
+    root_group.source_file_uuids = output.get_source_uuids([model_obj, obs_obj])
 
 
 def add_time_attribute(date: datetime) -> dict:
diff --git a/cloudnetpy/model_evaluation/metadata.py b/cloudnetpy/model_evaluation/metadata.py
index 1e3c8c82..d991f47d 100644
--- a/cloudnetpy/model_evaluation/metadata.py
+++ b/cloudnetpy/model_evaluation/metadata.py
@@ -17,13 +17,14 @@
     "horizontal_resolution": MetaData(
         long_name="Horizontal resolution of model",
         units="km",
-        comment="Distance between two grid point",
+        comment="Distance between two grid points",
         dimensions=None,
     ),
     "level": MetaData(
         long_name="Model level",
         units="1",
-        comment="Level 1 describes the highest height from ground.",
+        comment="Level 1 is the topmost model level; the level number increases "
+        "towards the surface.",
         axis="Z",
         positive="down",
         dimensions=("level",),
@@ -45,7 +46,7 @@
         long_name="Height above ground",
         units="m",
         comment=(
-            "Height have been calculated using pressure, temperature and\n"
+            "Height has been calculated using pressure, temperature and\n"
             "specific humidity."
         ),
         positive="up",
@@ -77,8 +78,8 @@
     ),
     "omega": MetaData(
         long_name="Vertical wind in pressure coordinates",
-        units="PA s-1",
-        standard_name="omega",
+        units="Pa s-1",
+        standard_name="lagrangian_tendency_of_air_pressure",
         dimensions=None,
     ),
     "q": MetaData(long_name="Specific humidity", units="1", dimensions=None),
@@ -102,22 +103,22 @@
     "cf_cirrus": MetaData(
         long_name="Cloud fraction of model grid point with filtered cirrus fraction",
         units="1",
-        comment="High level cirrus fraction is reduce do to lack if a radar\n"
-        "capability to observe correctly particles small and high.",
+        comment="High-level cirrus fraction is reduced due to the lack of radar\n"
+        "capability to correctly observe small particles at high altitudes.",
         dimensions=("time", "level"),
     ),
     "iwc": MetaData(
         long_name="Ice water content of model grid point",
         units="kg m-3",
-        comment="Calculated using model ice water mixing ration, pressure\n"
-        "and temperature: qi*P/287*T",
+        comment="Calculated using model ice water mixing ratio, pressure\n"
+        "and temperature: qi*P/(287*T)",
         dimensions=("time", "level"),
     ),
     "lwc": MetaData(
         long_name="Liquid water content of model grid point",
         units="kg m-3",
-        comment="Calculated using model liquid water mixing ration, pressure\n"
-        "and temperature:  ql*P/287*T",
+        comment="Calculated using model liquid water mixing ratio, pressure\n"
+        "and temperature: ql*P/(287*T)",
         dimensions=("time", "level"),
     ),
 }
@@ -127,9 +128,9 @@
         long_name="Observed cloud fraction by volume",
         units="1",
         comment=(
-            "Cloud fraction generated from observations and by volume,\n"
-            "averaged onto the models grid with height and time. Volume is\n"
-            "space within four grid points."
+            "Cloud fraction generated from observations by volume, averaged onto\n"
+            "the model grid in height and time. The volume value is the mean over\n"
+            "all observation pixels within the model grid cell."
         ),
         dimensions=("time", "level"),
     ),
@@ -137,9 +138,9 @@
         long_name="Observed cloud fraction by area",
         units="1",
         comment=(
-            "Cloud fraction generated from observation  and by area,\n"
-            "averaged onto the models grid with height and time. Area is\n"
-            "sum of time columns with any cloud fraction."
+            "Cloud fraction generated from observations and by area,\n"
+            "averaged onto the model grid with height and time. Area is\n"
+            "the fraction of time columns containing any cloud."
         ),
         dimensions=("time", "level"),
     ),
@@ -147,9 +148,9 @@
         long_name="Observed cloud fraction by advection volume",
         units="1",
         comment=(
-            "This variable is the same as the observed cloud fraction by volume, cf_V\n"
-            "except that model winds were used to estimate the time taken to advect\n"
-            "airflow a distance equivalent to the models horizontal resolution."
+            "This variable is the same as the observed cloud fraction by volume,\n"
+            "cf_V, except that model winds were used to estimate the time taken to\n"
+            "advect airflow a distance equivalent to the model's horizontal resolution."
         ),
         dimensions=("time", "level"),
     ),
@@ -157,46 +158,23 @@
         long_name="Observed cloud fraction by advection area",
         units="1",
         comment=(
-            "This variable is the same as the observed cloud fraction by area, cf_A\n"
+            "This variable is the same as the observed cloud fraction by area, cf_A,\n"
             "except that model winds were used to estimate the time taken to advect\n"
-            "airflow a distance equivalent to the models horizontal resolution."
+            "airflow a distance equivalent to the model's horizontal resolution."
         ),
         dimensions=("time", "level"),
     ),
     "iwc": MetaData(
-        long_name="Observed ice water content reshaped to model dimensions by"
-        " averaging",
+        long_name="Observed ice water content reshaped to model grid by averaging",
         units="kg m-3",
         comment=(
             "This variable is the observed mean ice water content derived from radar\n"
             "reflectivity factor averaged onto the model grid with height and time.\n"
-            "The formula  has been applied where the categorization data has\n"
+            "The formula has been applied where the categorization data has\n"
             "diagnosed that the radar echo is due to ice."
         ),
         dimensions=("time", "level"),
     ),
-    "iwc_att": MetaData(
-        long_name="Observed ice water content with attenuation reshaped to model grid"
-        " by averaging",
-        units="kg m-3",
-        comment=(
-            "This variable is the same as the observed mean ice water content, iwc,\n"
-            "except that profiles with uncorrected attenuation of the radar\n"
-            "reflectivity were included."
-        ),
-        dimensions=("time", "level"),
-    ),
-    "iwc_rain": MetaData(
-        long_name="Observed ice water content with raining reshaped to model grid"
-        " by averaging",
-        units="kg m-3",
-        comment=(
-            "This variable is the same as the observed mean ice water content\n"
-            "including attenuation, iwc_att, except that profiles with rain\n"
-            "at the surface were also included."
-        ),
-        dimensions=("time", "level"),
-    ),
     "iwc_adv": MetaData(
         long_name="Observed ice water content reshaped to model advection grid"
         " by averaging",
@@ -204,29 +182,7 @@
         comment=(
             "This variable is the same as the observed mean ice water content, iwc,\n"
             "except that model winds were used to estimate the time taken to advect\n"
-            "the flow a distance equivalent to the models horizontal resolution."
-        ),
-        dimensions=("time", "level"),
-    ),
-    "iwc_att_adv": MetaData(
-        long_name="Observed ice water content with attenuation reshaped to model"
-        " advection grid by averaging",
-        units="kg m-3",
-        comment=(
-            "This variable is the same as the observed mean ice water content,\n"
-            "iwc_adv,  except that profiles with uncorrected attenuation of the radar\n"
-            "reflectivity were included."
-        ),
-        dimensions=("time", "level"),
-    ),
-    "iwc_rain_adv": MetaData(
-        long_name="Observed ice water content with raining reshaped to model"
-        " advection rid by averaging",
-        units="kg m-3",
-        comment=(
-            "This variable is the same as the observed mean ice water content\n"
-            "including attenuation, iwc_att_adv, except that profiles with rain\n"
-            "at the surface were also included."
+            "the flow a distance equivalent to the model's horizontal resolution."
         ),
         dimensions=("time", "level"),
     ),
@@ -235,7 +191,7 @@
         units="kg m-3",
         comment=(
             "This variable is the observed mean liquid water content estimated for\n"
-            "pixels where the categorization data has diagnosed that liquid water \n"
+            "pixels where the categorization data has diagnosed that liquid water\n"
             "is present, averaged onto the model grid with height and time."
         ),
         dimensions=("time", "level"),
@@ -247,7 +203,8 @@
         comment=(
             "This variable is the same as the observed mean liquid water content,\n"
             "lwc, except that model winds were used to estimate the time taken to\n"
-            "advect flow over the models grid points."
+            "advect the flow a distance equivalent to the model's horizontal "
+            "resolution."
         ),
         dimensions=("time", "level"),
     ),
diff --git a/cloudnetpy/model_evaluation/model_metadata.py b/cloudnetpy/model_evaluation/model_metadata.py
index 8dfee90b..e9faa41e 100644
--- a/cloudnetpy/model_evaluation/model_metadata.py
+++ b/cloudnetpy/model_evaluation/model_metadata.py
@@ -1,50 +1,60 @@
-from typing import NamedTuple
-
-
-class ModelMetaData(NamedTuple):
-    long_name: str | None = None
-    cycle_var: str | None = None
-    common_var: str | None = None
-    cycle: str | None = None
-    level: int | None = None
-    model_name: str | None = None
-
-
-MODELS = {
-    "ecmwf": ModelMetaData(
-        model_name="ECMWF",
-        long_name="European Centre for Medium-Range Weather Forecasts",
-        level=88,
-    ),
-    "icon": ModelMetaData(
-        model_name="ICON-Iglo",
-        long_name="Icosahedral Nonhydrostatic Model",
-        level=62,
-        cycle="12-23, 24-35, 36-47",
-    ),
-    "era5": ModelMetaData(
-        model_name="ERA5",
-        long_name="Earth Re-Analysis System",
-        level=88,
-        cycle="1-12, 7-18",
-    ),
-    "harmonie-fmi-6-11": ModelMetaData(
-        model_name="HARMONIE-AROME",
-        long_name="the HIRLAM–ALADIN Research on Mesoscale Operational NWP in Euromed",
-        level=65,
-        cycle="6-11",
-    ),
+"""Model-agnostic configuration for the model evaluation subpackage.
+
+All model files consumed here are expected to be harmonized to a common set of
+variable names and units (temperature [K], pressure [Pa], qi/ql/
+cloud_fraction [1], uwind/vwind [m s-1], height [m]). There is therefore no
+per-model variable mapping: the only thing that differs between models is the
+number of vertical levels, which is handled dynamically via `ALTITUDE_LIMIT`.
+"""
+
+# Mapping from internal product keys to the harmonized model-file variable
+# names. These names are expected to be identical for every model.
+MODEL_VARIABLE_NAMES = {
+    "T": "temperature",
+    "p": "pressure",
+    "h": "height",
+    "iwc": "qi",
+    "lwc": "ql",
+    "cf": "cloud_fraction",
 }
 
-VARIABLES = {
-    "variables": ModelMetaData(
-        common_var="time, level, latitude, longitude, horizontal_resolution",
-        cycle_var="forecast_time, height",
-    ),
-    "T": ModelMetaData(long_name="temperature"),
-    "p": ModelMetaData(long_name="pressure"),
-    "h": ModelMetaData(long_name="height"),
-    "iwc": ModelMetaData(long_name="qi"),
-    "lwc": ModelMetaData(long_name="ql"),
-    "cf": ModelMetaData(long_name="cloud_fraction"),
+# Optional frozen-condensate mixing-ratio variables summed into the model ice
+# water content alongside the required cloud-ice field `qi`. Some models (e.g.
+# ECMWF) carry all grid-mean ice in `qi`, while others (e.g. HARMONIE) leave
+# `qi` empty and store the frozen mass in the snow and graupel categories. The
+# observed IWC includes this precipitating ice, so it must be included in the
+# model IWC too. Components absent from a given model file are simply skipped.
+OPTIONAL_IWC_VARIABLES = ("qs", "qg")
+
+# Prefix for the model's own (simulated) fields in the L3 output. Observation
+# fields downsampled to the model grid are stored without a prefix. The model
+# identity itself is stored in the file's global attributes, not in variable
+# names, so one L3 file holds exactly one model run.
+MODEL_PREFIX = "model_"
+
+# Name of the vertical dimension in the harmonized model files. The vertical
+# axis is expected under this dimension for every model.
+LEVEL_DIMENSION = "level"
+
+# Human-readable names for the observation products, used in file titles.
+PRODUCT_NAMES = {
+    "cf": "cloud fraction",
+    "iwc": "ice water content",
+    "lwc": "liquid water content",
 }
+
+# Coordinate / metadata variables that are identical between forecast cycles.
+COMMON_VARIABLES = (
+    "time",
+    "level",
+    "latitude",
+    "longitude",
+    "horizontal_resolution",
+)
+
+# Variables that may differ between forecast cycles.
+CYCLE_VARIABLES = ("forecast_time", "height")
+
+# Vertical levels above this altitude (m) are dropped: the radar cannot observe
+# them and the model grid above is not evaluated.
+ALTITUDE_LIMIT = 22000.0
diff --git a/cloudnetpy/model_evaluation/plotting/plot_tools.py b/cloudnetpy/model_evaluation/plotting/plot_tools.py
index dae0cd8d..f75d8d78 100644
--- a/cloudnetpy/model_evaluation/plotting/plot_tools.py
+++ b/cloudnetpy/model_evaluation/plotting/plot_tools.py
@@ -4,23 +4,42 @@
 from matplotlib.axes import Axes
 from numpy import ma
 
-from cloudnetpy.model_evaluation.model_metadata import MODELS
+from cloudnetpy.model_evaluation.model_metadata import MODEL_PREFIX
+
+
+def read_model_id(nc_file: str) -> str:
+    """Returns the model identifier (e.g. "ecmwf").
+
+    Stored as the `model_id` global attribute when the L3 file is created.
+    """
+    with netCDF4.Dataset(nc_file) as nc:
+        return nc.model_id
+
+
+def read_model_name(nc_file: str, model: str) -> str:
+    """Returns a human-readable model name for plot titles.
+
+    The model description is stored as the `model_name` global attribute when
+    the file is created (from the model file's `source`). Falls back to the
+    model identifier if the attribute is missing.
+    """
+    with netCDF4.Dataset(nc_file) as nc:
+        return getattr(nc, "model_name", model)
 
 
 def parse_wanted_names(
     nc_file: str,
     name: str,
-    model: str,
     variables: list | None = None,
     *,
     advance: bool = False,
 ) -> tuple[list, list]:
     """Returns standard and advection lists of product types to plot."""
-    names = variables or parse_dataset_keys(nc_file, name, advance=advance, model=model)
+    names = variables or parse_dataset_keys(nc_file, name, advance=advance)
     standard_n = [n for n in names if name in n and "adv" not in n]
-    standard_n = sort_model2first_element(standard_n, model)
+    standard_n = sort_model2first_element(standard_n)
     advection_n = [n for n in names if name in n and "adv" in n]
-    model_names = [n for n in names if f"{model}_" in n and f"_{model}_" not in n]
+    model_names = [n for n in names if n.startswith(MODEL_PREFIX)]
     for i, model_n in enumerate(model_names):
         advection_n.insert(0 + i, model_n)
     if len(advection_n) < len(standard_n):
@@ -30,68 +49,34 @@ def parse_wanted_names(
     return standard_n, advection_n
 
 
-def parse_dataset_keys(
-    nc_file: str,
-    product: str,
-    *,
-    advance: bool,
-    model: str,
-) -> list:
-    names = list(netCDF4.Dataset(nc_file).variables.keys())
-    a_names = ["cirrus", "snow"]
-    model_vars = []
-    for n in names:
-        if model not in n or (model in n and product not in n):
-            model_vars.append(n)
+def parse_dataset_keys(nc_file: str, product: str, *, advance: bool) -> list:
+    with netCDF4.Dataset(nc_file) as nc:
+        names = [n for n in nc.variables if product in n]
     if not advance:
-        for a in a_names:
-            for n in names:
-                if a in n:
-                    model_vars.append(n)
-    for m in model_vars:
-        names.remove(m)
+        names = [n for n in names if "cirrus" not in n and "snow" not in n]
     return names
 
 
-def sort_model2first_element(a: list, model: str) -> list:
-    mm = [n for n in a if f"{model}_" in n and f"_{model}_" not in n]
+def sort_model2first_element(a: list) -> list:
+    mm = [n for n in a if n.startswith(MODEL_PREFIX)]
     for i, m in enumerate(mm):
         a.remove(m)
         a.insert(0 + i, m)
     return a
 
 
-def sort_cycles(names: list, model: str) -> tuple[list, list]:
-    model_info = MODELS[model]
-    cycles = model_info.cycle
-    if cycles is None:
-        raise AttributeError
-    cycles_split = [x.strip() for x in cycles.split(",")]
-    cycles_names = [[name for name in names if cycle in name] for cycle in cycles_split]
-    cycles_names.sort()
-    cycles_names = [c for c in cycles_names if c]
-    cycles_new = [c for c in cycles_split for name in cycles_names if c in name[0]]
-    return cycles_names, cycles_new
-
-
-def read_data_characters(nc_file: str, name: str, model: str) -> tuple:
+def read_data_characters(nc_file: str, name: str) -> tuple:
     """Gets dimensions and data for plotting."""
-    nc = netCDF4.Dataset(nc_file)
-    data = nc.variables[name][:]
-    data = mask_small_values(data, name)
-    x = nc.variables["time"][:]
-    x = reshape_1d2nd(x, data)
-    try:
-        y = nc.variables[f"{model}_height"][:]
-    except KeyError as err:
-        model_info = MODELS[model]
-        cycles = model_info.cycle
-        if cycles is None:
-            msg = f"Invalid model: {model}"
+    with netCDF4.Dataset(nc_file) as nc:
+        data = nc.variables[name][:]
+        data = mask_small_values(data, name)
+        x = nc.variables["time"][:]
+        x = reshape_1d2nd(x, data)
+        try:
+            y = nc.variables[f"{MODEL_PREFIX}height"][:]
+        except KeyError as err:
+            msg = f"Missing variable {MODEL_PREFIX}height"
             raise RuntimeError(msg) from err
-        cycles_split = [x.strip() for x in cycles.split(",")]
-        cycle = [cycle for cycle in cycles_split if cycle in name]
-        y = nc.variables[f"{model}_{cycle[0]}_height"][:]
     y = y / 1000
     try:
         mask = y.mask
diff --git a/cloudnetpy/model_evaluation/plotting/plotting.py b/cloudnetpy/model_evaluation/plotting/plotting.py
index 6bd971e1..a68df614 100644
--- a/cloudnetpy/model_evaluation/plotting/plotting.py
+++ b/cloudnetpy/model_evaluation/plotting/plotting.py
@@ -1,7 +1,7 @@
 import logging
-import os
-import sys
+from dataclasses import dataclass
 from datetime import date
+from typing import TYPE_CHECKING
 
 import matplotlib as mpl
 import matplotlib.pyplot as plt
@@ -18,18 +18,44 @@
 from numpy import ma
 
 import cloudnetpy.model_evaluation.plotting.plot_tools as p_tools
-from cloudnetpy.model_evaluation.model_metadata import MODELS
 from cloudnetpy.model_evaluation.plotting.plot_meta import ATTRIBUTES, PlotMeta
 from cloudnetpy.model_evaluation.statistics.statistical_methods import DayStatistics
-from cloudnetpy.plotting.plotting import Dimensions, get_log_cbar_tick_labels, lin2log
+from cloudnetpy.plotting.plotting import (
+    Dimensions,
+    add_subtitle,
+    get_log_cbar_tick_labels,
+    get_time_tick_labels,
+    lin2log,
+)
 
-sys.path.append(os.path.dirname(os.path.abspath(__file__)))
+if TYPE_CHECKING:
+    from collections.abc import Callable
+
+
+@dataclass
+class L3PlotConfig:
+    """Shared context passed to the ``get_*_plots`` helpers.
+
+    Bundles the fields that are constant across a single
+    ``generate_L3_day_plots`` call so the helpers take just this config and the
+    per-figure ``names``.
+    """
+
+    product: str
+    nc_file: str
+    model: str
+    model_name: str
+    save_path: str | None
+    image_name: str | None
+    stats: tuple | None = None
+    show: bool = False
+    title: bool = True
+    include_xlimits: bool = False
 
 
 def generate_L3_day_plots(
     nc_file: str,
     product: str,
-    model: str,
     *,
     title: bool = True,
     var_list: list | None = None,
@@ -38,26 +64,26 @@ def generate_L3_day_plots(
     save_path: str | None = None,
     image_name: str | None = None,
     show: bool | None = False,
+    include_advection: bool = True,
 ) -> list[Dimensions]:
     """Generate visualizations for level 3 dayscale products.
-    With figure type visualizations can be subplot in group, pair, single or
-    statistic of given product. In group fig_type all different methods are plot
-    in same figure but standard timegrid is separated from advection timegrid.
-    In pair fig_type upper subplot is always model product and below one is
-    product method. All product method in given file will be plotted in loop.
-    Single fig_type will plot each product variable in a own figure.
-    Statistic fig_type will plot select statistical method of all product method
-    in same fig.
+
+    Two kinds of figure are produced. The default 'group' figure draws the
+    requested variables (the model field and its downsampling methods) as
+    colormesh subplots in a single figure - one subplot if a single variable is
+    requested via ``var_list``, many otherwise. The 'statistic' figure instead
+    plots a model-vs-observation comparison (relative error, area, histogram or
+    vertical profile). In both cases the standard and advection timegrids are
+    drawn into separate figures.
 
     Args:
         nc_file (str): Path to source file
         product (str): Name of product wanted to plot
-        model (str): Name of model which downsampling was done with
 
         title (bool): Show title and additional labels. Operational processing
              sets this to False. Default is True.
-        fig_type (str, optional): Type of figure wanted to produce. Options
-            are 'group', 'pair', 'single' and 'statistic'. Default value is 'group'
+        fig_type (str, optional): Type of figure to produce, 'group' (default)
+            or 'statistic'.
         stats (list, optional): List of statistical methods to visualize in
             'statistic' fig_type generation. Default is all types, but methods can
             be called individually.
@@ -68,285 +94,116 @@ def generate_L3_day_plots(
 
         image_name (str, optional): Saving name of generated fig
         show (bool, optional): If True, shows visualization
+        include_advection (bool, optional): If True (default), also plot the
+            wind-advection time grid in its own figure(s). If False, plot only
+            the standard time grid.
+
     Notes:
         In case of 'group' and 'statistic' fig_type advection timegrid is
         separated from standard timegrid to their own figures.
-        In case of model cycles, cycles are visualized in their on figures same
-        way as an individual model run would be visualized in its own in a group
-        figure.
 
     Examples:
         >>> from cloudnetpy.model_evaluation.plotting.plotting
         import generate_L3_day_plots
         >>> l3_day_file = 'cf_ecmwf.nc'
         >>> product = 'cf'
-        >>> model = 'ecmwf'
-        >>> generate_L3_day_plots(l3_day_file, product, model)
-        >>> l3_day_file = 'cf_ecmwf.nc'
-        >>> product = 'cf'
-        >>> model = 'ecmwf'
-        >>> generate_L3_day_plots(l3_day_file, product, model,
+        >>> generate_L3_day_plots(l3_day_file, product)
+        >>> generate_L3_day_plots(l3_day_file, product,
         >>>                       fig_type='statistic', stats=['error'])
     """
-
-    def _check_cycle_names() -> None:
-        if not c_names:
-            raise AttributeError
-
-    cls = __import__("plotting")
-    model_info = MODELS[model]
-    model_name = model_info.model_name
-    name_set = p_tools.parse_wanted_names(nc_file, product, model, var_list)
-    unique_tuples = {tuple(lst) for lst in name_set}
-    name_set_unique = tuple(list(tup) for tup in unique_tuples)
+    plotters: dict[str, Callable[[L3PlotConfig, list], tuple]] = {
+        "group": get_group_plots,
+        "statistic": get_statistic_plots,
+    }
+    if fig_type not in plotters:
+        msg = f"Invalid fig_type: {fig_type}"
+        raise ValueError(msg)
+    plot_function = plotters[fig_type]
+    model = p_tools.read_model_id(nc_file)
+    config = L3PlotConfig(
+        product=product,
+        nc_file=nc_file,
+        model=model,
+        model_name=p_tools.read_model_name(nc_file, model),
+        save_path=save_path,
+        image_name=image_name,
+        stats=stats,
+        show=bool(show),
+        title=title,
+    )
+    standard_names, advection_names = p_tools.parse_wanted_names(
+        nc_file, product, var_list
+    )
+    name_set = (
+        [standard_names, advection_names] if include_advection else [standard_names]
+    )
+    name_set_unique: list[list] = []
+    for names in name_set:
+        if names not in name_set_unique:
+            name_set_unique.append(names)
 
     dimensions = []
     for names in name_set_unique:
-        if len(names) > 0:
-            try:
-                cycle_names, cycles = p_tools.sort_cycles(names, model)
-                for i, c_names in enumerate(cycle_names):
-                    _check_cycle_names()
-                    params = [
-                        product,
-                        c_names,
-                        nc_file,
-                        model,
-                        model_name,
-                        save_path,
-                        image_name,
-                    ]
-                    kwargs = {"show": show, "title": title, "cycle": cycles[i]}
-                    if fig_type == "statistic":
-                        params = [
-                            product,
-                            c_names,
-                            nc_file,
-                            model,
-                            model_name,
-                            stats,
-                            save_path,
-                            image_name,
-                        ]
-                    figs, axes = getattr(cls, f"get_{fig_type}_plots")(
-                        *params, **kwargs
-                    )
-                    for fig, ax in zip(figs, axes, strict=False):
-                        dimensions.append(Dimensions(fig, ax))
-            except AttributeError:
-                params = [
-                    product,
-                    names,
-                    nc_file,
-                    model,
-                    model_name,
-                    save_path,
-                    image_name,
-                ]
-                kwargs = {"show": show, "title": title}
-                if fig_type == "statistic":
-                    params = [
-                        product,
-                        names,
-                        nc_file,
-                        model,
-                        model_name,
-                        stats,
-                        save_path,
-                        image_name,
-                    ]
-                figs, axes = getattr(cls, f"get_{fig_type}_plots")(*params, **kwargs)
-                for fig, ax in zip(figs, axes, strict=False):
-                    dimensions.append(Dimensions(fig, ax))
+        if len(names) == 0:
+            continue
+        figs, axes = plot_function(config, names)
+        if isinstance(figs, Figure):  # group plots return a single figure
+            dimensions.append(Dimensions(figs, axes))
+        else:
+            for fig, ax in zip(figs, axes, strict=False):
+                dimensions.append(Dimensions(fig, ax))
     return dimensions
 
 
-def get_group_plots(
-    product: str,
-    names: list,
-    nc_file: str,
-    model: str,
-    model_name: str,
-    save_path: str,
-    image_name: str,
-    *,
-    show: bool,
-    cycle: str = "",
-    title: bool = True,
-    include_xlimits: bool = False,
-) -> tuple:
-    """Group subplot visualization for both standard and advection downsampling.
-    Generates group subplot figure for product with model and all different
-    downsampling methods. Generates separated figures for standard and advection
-    timegrids. All model cycles if any will be generated to their own figures.
+def get_group_plots(config: L3PlotConfig, names: list) -> tuple:
+    """Draws the requested variables as colormesh subplots in a single figure.
+
+    Produces one subplot per name (a single subplot if only one variable was
+    requested), with the model field in the first subplot. Standard and
+    advection timegrids are drawn into separate figures by the caller.
 
     Args:
-        product (str): Name of the product
-        names (list): List of variables to be visualized to same fig
-        nc_file (str): Path to a source file
-        model (str): Name of used model in a downsampling process
-        model_name (str): Correct name of a model
-        save_path (str): Path for saving figures
-        image_name (str, optional): Saving name of generated fig
-        show (bool): Show figure before saving if True
-        cycle (str): Name of cycle if exists
-        title (bool): True or False if wanted to add title to subfig
-        include_xlimits (bool): Show labels at the ends of x-axis
+        config (L3PlotConfig): Shared plotting context
+        names (list): Variables to draw as subplots in the same figure
     """
+    product = config.product
     fig, ax = initialize_figure(len(names))
-    model_run = model
     name = ""
     j = 0
     for j, name in enumerate(names):
         variable_info = ATTRIBUTES[product]
         set_yax(ax[j], 12, ylabel=None)
-        set_xax(ax[j], include_xlimits=include_xlimits)
-        if title:
+        set_xax(ax[j], include_xlimits=config.include_xlimits)
+        if config.title:
             _set_title(ax[j], name, product, variable_info)
-        if j == 0 and title:
-            _set_title(ax[j], model, product, variable_info, model_name)
-        data, x, y = p_tools.read_data_characters(nc_file, name, model)
+        if j == 0 and config.title:
+            _set_title(ax[j], config.model, product, variable_info, config.model_name)
+        data, x, y = p_tools.read_data_characters(config.nc_file, name)
         plot_colormesh(ax[j], data, (x, y), variable_info)
-    casedate = set_labels(fig, ax[j], nc_file)
+    casedate = set_labels(fig, ax[j], config.nc_file)
     if "adv" in name:
         product = product + "_adv"
-    if len(cycle) > 1:
-        fig.text(0.64, 0.885, f"Cycle: {cycle}", fontsize=13)
-        model_run = f"{model}_{cycle}"
     handle_saving(
-        image_name,
-        save_path,
+        config.image_name,
+        config.save_path,
         casedate,
-        [product, model_run, "group"],
-        show=show,
+        [product, config.model, "group"],
+        show=config.show,
     )
     return fig, ax
 
 
-def get_pair_plots(
-    product: str,
-    names: list,
-    nc_file: str,
-    model: str,
-    model_name: str,
-    save_path: str,
-    image_name: str,
-    cycle: str = "",
-    *,
-    show: bool = False,
-    title: bool = True,
-    include_xlimits: bool = False,
-) -> None:
-    """Pair subplots of model and product method.
-    In upper subplot is model product and lower subplot one of the
-    downsampled method of select product. Function generates all product methods
-    in a given nc-file in loop.
-
-    Args:
-        product (str): Name of the product
-        names (list): List of variables to be visualized to same fig
-        nc_file (str): Path to a source file
-        model (str): Name of used model in a downsampling process
-        model_name (str): Correct name of a model
-        save_path (str): Path for saving figures
-        image_name (str, optional): Saving name of generated fig
-        show (bool): Show figure before saving if True
-        cycle (str): Name of cycle if exists
-        title (bool): True or False if wanted add title to subfig
-        include_xlimits (bool): Show labels at the ends of x-axis
-    """
-    variable_info = ATTRIBUTES[product]
-    model_ax = names[0]
-    for i, name in enumerate(names):
-        if i == 0:
-            continue
-        fig, ax = initialize_figure(2)
-        set_yax(ax[0], 12, ylabel=None)
-        set_yax(ax[-1], 12, ylabel=None)
-        set_xax(ax[0], include_xlimits=include_xlimits)
-        set_xax(ax[-1], include_xlimits=include_xlimits)
-
-        if title:
-            _set_title(ax[0], model, product, variable_info, model_name)
-            _set_title(ax[-1], name, product, variable_info)
-        model_data, mx, my = p_tools.read_data_characters(nc_file, model_ax, model)
-        data, x, y = p_tools.read_data_characters(nc_file, name, model)
-        plot_colormesh(ax[0], model_data, (mx, my), variable_info)
-        plot_colormesh(ax[-1], data, (x, y), variable_info)
-        casedate = set_labels(fig, ax[-1], nc_file)
-        if len(cycle) > 1:
-            fig.text(0.64, 0.889, f"Cycle: {cycle}", fontsize=13)
-        handle_saving(
-            image_name,
-            save_path,
-            casedate,
-            [name, "pair"],
-            show=show,
-        )
-
-
-def get_single_plots(
-    product: str,
-    names: list,
-    nc_file: str,
-    model: str,
-    model_name: str,
-    save_path: str,
-    image_name: str,
-    *,
-    show: bool,
-    cycle: str = "",
-    title: bool = True,
-    include_xlimits: bool = False,
-) -> tuple[list, list]:
-    """Generates figures of each product variable from given file in loop.
-
-    Args:
-        product (str): Name of the product
-        names (list): List of variables to be visualized to same fig
-        nc_file (str): Path to a source file
-        model (str): Name of used model in a downsampling process
-        model_name (str): Correct name of a model
-        save_path (str): Path for saving figures
-        image_name (str, optional): Saving name of generated fig
-        show (bool): Show figure before saving if True
-        cycle (str): Name of cycle if exists
-        title (bool): True or False if wanted to add title to subfig
-        include_xlimits (bool): Show labels at the ends of x-axis
-    """
-    figs = []
-    axes = []
-    variable_info = ATTRIBUTES[product]
-    for name in names:
-        fig, ax = initialize_figure(1)
-        figs.append(fig)
-        axes.append(ax)
-        set_yax(ax[0], 12, ylabel=None)
-        set_xax(ax[0], include_xlimits=include_xlimits)
-
-        if title:
-            _set_title(ax[0], name, product, variable_info)
-        data, x, y = p_tools.read_data_characters(nc_file, name, model)
-        plot_colormesh(ax[0], data, (x, y), variable_info)
-        casedate = set_labels(fig, ax[0], nc_file, sub_title=title)
-        if title:
-            if len(cycle) > 1:
-                fig.text(0.64, 0.9, f"{model_name} cycle: {cycle}", fontsize=13)
-            else:
-                fig.text(0.64, 0.9, f"{model_name}", fontsize=13)
-        handle_saving(
-            image_name,
-            save_path,
-            casedate,
-            [name, "single"],
-            show=show,
-        )
-    return figs, axes
-
-
 def plot_colormesh(
     ax: Axes, data: npt.NDArray, axes: tuple, variable_info: PlotMeta
 ) -> None:
+    x, y = axes
+    if getattr(y, "ndim", 1) > 1:
+        # Collapse the 2D (time, level) model grid to 1D monotonic cell centers
+        # (mean height over time) so pcolormesh does not warn about non-monotonic
+        # coordinates.
+        x = x[:, 0]
+        y = ma.mean(y, axis=0)
+        data = data.T
     vmin, vmax = variable_info.plot_range
     if variable_info.clabel == "\u00b0C":
         data = data - 273.15
@@ -354,7 +211,7 @@ def plot_colormesh(
         data, vmin, vmax = lin2log(data, vmin, vmax)
     cmap = plt.get_cmap(variable_info.cbar, 22)
     data[data < vmin] = ma.masked
-    pl = ax.pcolormesh(*axes, data, vmin=vmin, vmax=vmax, cmap=cmap)
+    pl = ax.pcolormesh(x, y, data, vmin=vmin, vmax=vmax, cmap=cmap)
     colorbar = init_colorbar(pl, ax)
     if variable_info.plot_scale == "logarithmic":
         tick_labels = get_log_cbar_tick_labels(vmin, vmax)
@@ -365,43 +222,20 @@ def plot_colormesh(
         colorbar.set_label(variable_info.clabel, fontsize=13)
 
 
-def get_statistic_plots(
-    product: str,
-    names: list,
-    nc_file: str,
-    model: str,
-    model_name: str,
-    stats: list,
-    save_path: str,
-    image_name: str,
-    *,
-    show: bool,
-    cycle: str = "",
-    title: bool = True,
-) -> tuple:
+def get_statistic_plots(config: L3PlotConfig, names: list) -> tuple:
     """Statistical subplots for day scale products.
     Statistical analysis can be done by day scale with relative error ('error'),
     total data area analysis ('area'), histogram ('hist') or vertical profiles
     ('vertical'). Each given stats are looped through and generated as one figure
     per statistical method for a select product. All different downsampled method
-    are in a same fig. Standard and advection timegrids are separated to own figs
-    as well as different cycle runs.
+    are in a same fig. Standard and advection timegrids are separated to own figs.
 
     Args:
-        product (str): Name of the product
+        config (L3PlotConfig): Shared plotting context, including the list of
+            statistical methods to process in ``config.stats``
         names (list): List of variables to be visualized to same fig
-        nc_file (str): Path to a source file
-        model (str): Name of used model in a downsampling process
-        model_name (str): Correct name of a model
-        stats (list): List of statistical method to process analysis with.
-                      Options are ['error', 'area', 'hist', 'vertical']
-        save_path (str): Path for saving figures
-        image_name (str, optional): Saving name of generated fig
-        show (bool): Show figure before saving if True
-        cycle (str): Name of cycle if exists
-        title (bool): True or False if wanted to add title to subfig
     """
-    model_run = model
+    stats = config.stats or ()
     name = ""
     j = 0
 
@@ -410,11 +244,11 @@ def _check_data() -> None:
             _raise()
 
     def _check_data2() -> None:
-        if "error" in stat and np.all(day_stat.model_stat.mask is True):
+        if "error" in stat and np.all(ma.getmaskarray(day_stat.model_stat)):
             _raise()
 
     def _raise() -> None:
-        err_msg = f"No data in {model_name} or observation"
+        err_msg = f"No data in {config.model_name} or observation"
         raise ValueError(err_msg)
 
     figs = []
@@ -423,25 +257,27 @@ def _raise() -> None:
         try:
             obs_missing = False
             model_missing = False
-            variable_info = ATTRIBUTES[product]
+            variable_info = ATTRIBUTES[config.product]
             fig, ax = initialize_figure(len(names) - 1, stat)
             figs.append(fig)
             axes.append(ax)
-            model_data, _, _ = p_tools.read_data_characters(nc_file, names[0], model)
-            if np.all(model_data.mask is True):
+            model_data, _, _ = p_tools.read_data_characters(config.nc_file, names[0])
+            if np.all(ma.getmaskarray(model_data)):
                 model_missing = True
             for j, name in enumerate(names):
-                data, x, y = p_tools.read_data_characters(nc_file, name, model)
-                if np.all(data.mask is True):
+                data, x, y = p_tools.read_data_characters(config.nc_file, name)
+                if np.all(ma.getmaskarray(data)):
                     obs_missing = True
                 _check_data()
-                statistics = "aerror" if product == "cf" and stat == "error" else stat
+                statistics = (
+                    "aerror" if config.product == "cf" and stat == "error" else stat
+                )
                 if j > 0:
                     name_new = ""
-                    name_new = _get_stat_titles(name_new, product, variable_info)
+                    name_new = _get_stat_titles(name_new, config.product, variable_info)
                     day_stat = DayStatistics(
                         statistics,
-                        [product, model_name, name_new],
+                        [config.product, config.model_name, name_new],
                         model_data,
                         data,
                     )
@@ -456,26 +292,25 @@ def _raise() -> None:
                         data,
                         (x, y),
                         variable_info,
-                        title=title,
+                        title=config.title,
                     )
         except ValueError:
             logging.exception("Exception occurred")
         if stat not in ("hist", "vertical"):
-            casedate = set_labels(fig, ax[j - 1], nc_file, sub_title=title)
+            casedate = set_labels(
+                fig, ax[j - 1], config.nc_file, sub_title=config.title
+            )
         else:
-            casedate = read_date(nc_file)
-            _name = read_location(nc_file)
-            if title:
+            casedate = read_date(config.nc_file)
+            _name = read_location(config.nc_file)
+            if config.title:
                 add_subtitle(fig, casedate, _name.capitalize())
-        if len(cycle) > 1:
-            fig.text(0.64, 0.885, f"Cycle: {cycle}", fontsize=13)
-            model_run = f"{model}_{cycle}"
         handle_saving(
-            image_name,
-            save_path,
+            config.image_name,
+            config.save_path,
             casedate,
-            [name, stat, model_run],
-            show=show,
+            [name, stat, config.model],
+            show=config.show,
         )
     return figs, axes
 
@@ -497,7 +332,7 @@ def initialize_statistic_plots(
     if method in ("error", "aerror"):
         plot_relative_error(ax, day_stat.model_stat.T, args, method)
         if title:
-            ax.set_title(day_stat.title, fontsize=14)
+            ax.set_title(f"{day_stat.title}", fontsize=14)
         set_yax(ax, 12, ylabel=None)
         set_xax(ax, include_xlimits=include_xlimits)
 
@@ -691,7 +526,7 @@ def plot_vertical_profile(
         "-",
         color="orange",
         lw=2,
-        label="Mean of {day_stat.title[0]}",
+        label=f"Mean of {day_stat.title[0]}",
     )
     ax.plot(orm, axes, "-", color="green", lw=2, label="Mean of observation")
 
@@ -770,8 +605,6 @@ def _set_title(
         title = _get_product_title(variable_info)
         if product == "cf":
             title = _get_cf_title(field_name, variable_info)
-        if product == "iwc":
-            title = _get_iwc_title(field_name, variable_info)
         if "adv" in field_name:
             adv = " Downsampled using advection time"
             ax.text(0.9, -0.13, adv, size=12, ha="center", transform=ax.transAxes)
@@ -779,7 +612,7 @@ def _set_title(
     else:
         name = variable_info.name
         if len(model_name) > 1:
-            ax.set_title(f"{name} of {model_name}", fontsize=14)
+            ax.set_title(f"{name}, {model_name}", fontsize=14)
         else:
             ax.set_title(f"Simulated {name}")
 
@@ -791,17 +624,6 @@ def _get_cf_title(field_name: str, variable_info: PlotMeta) -> str:
     return title
 
 
-def _get_iwc_title(field_name: str, variable_info: PlotMeta) -> str:
-    name = variable_info.name
-    if "att" in field_name:
-        title = f"{name} with good attenuation"
-    elif "rain" in field_name:
-        title = f"{name} with rain"
-    else:
-        title = f"{name}"
-    return title
-
-
 def _get_product_title(variable_info: PlotMeta) -> str:
     return f"{variable_info.name}"
 
@@ -810,8 +632,6 @@ def _get_stat_titles(field_name: str, product: str, variable_info: PlotMeta) ->
     title = _get_product_title_stat(variable_info)
     if product == "cf":
         title = _get_cf_title_stat(field_name, variable_info)
-    if product == "iwc":
-        title = _get_iwc_title_stat(field_name, variable_info)
     if "adv" in field_name:
         adv = " (Advection time)"
         return f"{title}{adv}"
@@ -826,17 +646,6 @@ def _get_cf_title_stat(field_name: str, variable_info: PlotMeta) -> str:
     return title
 
 
-def _get_iwc_title_stat(field_name: str, variable_info: PlotMeta) -> str:
-    name = variable_info.name
-    if "att" in field_name:
-        title = f"{name} with good attenuation"
-    elif "rain" in field_name:
-        title = f"{name} with rain"
-    else:
-        title = f"{name}"
-    return title
-
-
 def _get_product_title_stat(variable_info: PlotMeta) -> str:
     name = variable_info.name
     return f"{name}"
@@ -852,29 +661,12 @@ def set_yax(ax: Axes, max_y: float, ylabel: str | None, min_y: float = 0.0) -> N
 
 def set_xax(ax: Axes, *, include_xlimits: bool = False) -> None:
     """Sets xticks and xtick labels for plt.imshow()."""
-    ticks_x_labels = _get_standard_time_ticks(include_xlimits=include_xlimits)
+    ticks_x_labels = get_time_tick_labels(edge_tick_labels=include_xlimits)
     ax.set_xticks(np.arange(0, 25, 4, dtype=int))
     ax.set_xticklabels(ticks_x_labels, fontsize=12)
     ax.set_xlim(0, 24)
 
 
-def _get_standard_time_ticks(
-    resolution: int = 4,
-    *,
-    include_xlimits: bool = False,
-) -> list:
-    """Returns typical ticks / labels for a time vector between 0-24h."""
-    if include_xlimits:
-        return [
-            f"{int(i):02d}:00" if 24 >= i >= 0 else ""
-            for i in np.arange(0, 24.01, resolution)
-        ]
-    return [
-        f"{int(i):02d}:00" if 24 > i > 0 else ""
-        for i in np.arange(0, 24.01, resolution)
-    ]
-
-
 def set_labels(fig: Figure, ax: Axes, nc_file: str, *, sub_title: bool = True) -> date:
     ax.set_xlabel("Time (UTC)", fontsize=13)
     case_date = read_date(nc_file)
@@ -896,24 +688,6 @@ def read_date(nc_file: str) -> date:
         return date(int(nc.year), int(nc.month), int(nc.day))
 
 
-def add_subtitle(fig: Figure, case_date: date, site_name: str) -> None:
-    """Adds subtitle into figure."""
-    text = _get_subtitle_text(case_date, site_name)
-    fig.suptitle(
-        text,
-        fontsize=13,
-        y=0.885,
-        x=0.07,
-        horizontalalignment="left",
-        verticalalignment="bottom",
-    )
-
-
-def _get_subtitle_text(case_date: date, site_name: str) -> str:
-    site_name = site_name.replace("-", " ")
-    return f"{site_name}, {case_date.strftime('%d %b %Y').lstrip('0')}"
-
-
 def _create_save_name(
     save_path: str,
     case_date: date,
diff --git a/cloudnetpy/model_evaluation/products/advance_methods.py b/cloudnetpy/model_evaluation/products/advance_methods.py
index 28473c60..906a1cc4 100644
--- a/cloudnetpy/model_evaluation/products/advance_methods.py
+++ b/cloudnetpy/model_evaluation/products/advance_methods.py
@@ -1,7 +1,3 @@
-import importlib
-import logging
-from os import PathLike
-
 import numpy as np
 import numpy.typing as npt
 import scipy.special
@@ -9,12 +5,17 @@
 from numpy import ma
 
 import cloudnetpy.utils as cl_tools
-from cloudnetpy.datasource import DataSource
+from cloudnetpy.model_evaluation.model_metadata import MODEL_PREFIX
 from cloudnetpy.model_evaluation.products.model_products import ModelManager
 from cloudnetpy.model_evaluation.products.observation_products import ObservationManager
+from cloudnetpy.products.product_tools import (
+    IceCoefficients,
+    get_ice_coefficients,
+    z_to_iwc,
+)
 
 
-class AdvanceProductMethods(DataSource):
+class AdvanceProductMethods:
     """Class that adds advance methods of product to nc-file.
     Different methods could be filtering or adding info by making
     assumptions of model or observation data.
@@ -27,38 +28,30 @@ class AdvanceProductMethods(DataSource):
     def __init__(
         self,
         model_obj: ModelManager,
-        model_file: str | PathLike,
         obs_obj: ObservationManager,
     ) -> None:
-        super().__init__(model_file)
         self._obs_obj = obs_obj
-        self.product = obs_obj.obs
+        self.product = obs_obj.product
         self._date = obs_obj.date
         self._obs_height = obs_obj.data["height"][:]
-        self._obs_data = obs_obj.data[obs_obj.obs][:]
+        self._obs_data = obs_obj.data[obs_obj.product][:]
         self._model_obj = model_obj
         self._model_time = model_obj.time
         self._model_height = model_obj.data[model_obj.keys["height"]][:]
         self.generate_products()
 
     def generate_products(self) -> None:
-        cls = importlib.import_module(__name__).AdvanceProductMethods
-        try:
-            name = f"get_advance_{self.product}"
-            getattr(cls, name)(self)
-        except AttributeError as error:
-            logging.debug("No advance method for %s: %s", self.product, error)
-
-    def get_advance_cf(self) -> None:
-        self.cf_cirrus_filter()
+        # Cloud fraction is the only product with an advance (cirrus) step.
+        if self.product == "cf":
+            self.cf_cirrus_filter()
 
     def cf_cirrus_filter(self) -> None:
         cf = self.getvar_from_object("cf")
         h = self.getvar_from_object("h")
-        temperature = self.getvar("temperature")
+        temperature = self._model_obj.getvar("temperature")
         t_screened = self.remove_extra_levels(temperature - 273.15)
         iwc, lwc = (self._model_obj.get_water_content(var) for var in ["iwc", "lwc"])
-        tZT, tT, tZ, t = self.set_frequency_parameters()
+        coeffs = self.set_frequency_parameters()
         z_sen = self.fit_z_sensitivity(h)
         cf_filtered = self.filter_high_iwc_low_cf(cf, iwc, lwc)
         cloud_iwc, ice_ind = self.find_ice_in_clouds(cf_filtered, iwc, lwc)
@@ -82,39 +75,27 @@ def cf_cirrus_filter(self) -> None:
                 cf_filtered[ind] = ma.masked
                 continue
 
-            min_iwc = 10 ** (
-                tZT * z_sen[ind] * t_screened[ind]
-                + tT * t_screened[ind]
-                + tZ * z_sen[ind]
-                + t
-            )
+            min_iwc = z_to_iwc(coeffs, z_sen[ind], t_screened[ind])
             obs_index = iwc_dist > min_iwc
             cf_filtered[ind] = self.filter_cirrus(p_iwc, obs_index, cf_filtered[ind])
 
         cf_filtered[cf_filtered < 0.05] = ma.masked
 
-        self._model_obj.append_data(
-            cf_filtered,
-            f"{self._model_obj.model}{self._model_obj.cycle}_cf_cirrus",
-        )
+        self._model_obj.append_data(cf_filtered, f"{MODEL_PREFIX}cf_cirrus")
 
     def getvar_from_object(self, arg: str) -> npt.NDArray:
         v_name = arg if arg == "cf" else self._model_obj.get_model_var_names((arg,))[0]
-        key = f"{self._model_obj.model}{self._model_obj.cycle}_{v_name}"
-        return self._model_obj.data[key][:]
+        return self._model_obj.data[f"{MODEL_PREFIX}{v_name}"][:]
 
     def remove_extra_levels(self, arg: npt.NDArray) -> npt.NDArray:
         return self._model_obj.cut_off_extra_levels(arg)
 
-    def set_frequency_parameters(self) -> tuple:
+    def set_frequency_parameters(self) -> IceCoefficients:
         if self._obs_obj.radar_freq is None:
             msg = "No radar frequency in observation file"
             raise ValueError(msg)
-        if 30 <= self._obs_obj.radar_freq <= 40:
-            return 0.000242, -0.0186, 0.0699, -1.63
-        if 90 <= float(self._obs_obj.radar_freq) <= 100:
-            return 0.00058, -0.00706, 0.0923, -0.992
-        raise ValueError
+        wl_band = cl_tools.get_wl_band(float(self._obs_obj.radar_freq))
+        return get_ice_coefficients("iwc", wl_band)
 
     def fit_z_sensitivity(self, h: npt.NDArray) -> npt.NDArray:
         if self._obs_obj.z_sensitivity is None:
@@ -172,8 +153,8 @@ def get_ice_indices(
         return tuple(np.where((cf_filtered > 0) & (iwc > 0) & (lwc < iwc / 10)))
 
     def iwc_variance(self, height: npt.NDArray, ice_ind: tuple) -> npt.NDArray:
-        u = self.getvar("uwind")
-        v = self.getvar("vwind")
+        u = self._model_obj.getvar("uwind")
+        v = self._model_obj.getvar("vwind")
         u = self.remove_extra_levels(u)
         v = self.remove_extra_levels(v)
         w_shear = self.calculate_wind_shear(self._model_obj.wind, u, v, height)
@@ -199,9 +180,7 @@ def calculate_wind_shear(
         for w in (wind, u, v):
             grad_w = np.zeros(w.shape)
             grad_w[0, :] = (w[1, :] - w[0, :]) / (height[1, :] - height[0, :])
-            grad_w[1:-2, :] = (w[2:-1, :] - 2 * w[1:-2, :] + w[1:-2, :]) / (
-                height[2:-1, :] - height[1:-2, :]
-            )
+            grad_w[1:-1, :] = (w[2:, :] - w[:-2, :]) / (height[2:, :] - height[:-2, :])
             grad_w[-1, :] = (w[-1, :] - w[-2, :]) / (height[-1, :] - height[-2, :])
             grand_winds.append(grad_w)
 
@@ -222,7 +201,7 @@ def calculate_iwc_distribution(
         iwc_dist = np.arange(0, finish, finish / (n_dist - 1))
         if cloud_iwc < iwc_dist[2]:
             finish = cloud_iwc * 10
-            iwc_dist = np.arange(0, finish, finish / n_dist - 1)
+            iwc_dist = np.arange(0, finish, finish / (n_dist - 1))
         return iwc_dist
 
     @staticmethod
@@ -240,3 +219,11 @@ def filter_cirrus(
         cf_filtered: npt.NDArray,
     ) -> npt.NDArray:
         return (ma.sum(p_iwc * obs_index) / ma.sum(p_iwc)) * cf_filtered
+
+
+def apply_cirrus_filter(model_obj: ModelManager, obs_obj: ObservationManager) -> None:
+    """Add the cirrus-filtered model cloud fraction to `model_obj` in place.
+
+    Raises ValueError if the model has no ice clouds to filter.
+    """
+    AdvanceProductMethods(model_obj, obs_obj)
diff --git a/cloudnetpy/model_evaluation/products/grid_methods.py b/cloudnetpy/model_evaluation/products/grid_methods.py
index db0e7949..364209d0 100644
--- a/cloudnetpy/model_evaluation/products/grid_methods.py
+++ b/cloudnetpy/model_evaluation/products/grid_methods.py
@@ -25,12 +25,12 @@ def __init__(self, model_obj: ModelManager, obs_obj: ObservationManager) -> None
         self._date = obs_obj.date
         self._obs_time = tl.time2datetime(obs_obj.time, self._date)
         self._obs_height = obs_obj.data["height"][:]
-        self._obs_data = obs_obj.data[obs_obj.obs][:]
+        self._obs_data = obs_obj.data[obs_obj.product][:]
         self.model_obj = model_obj
         self._model_time = model_obj.time
         self._model_height = model_obj.data[model_obj.keys["height"]][:]
         self._time_adv = tl.calculate_advection_time(
-            int(model_obj.resolution_h),
+            model_obj.resolution_h,
             ma.array(model_obj.wind),
             1,
         )
@@ -49,12 +49,6 @@ def _generate_downsample_product(self) -> None:
                 (self._time_steps[i], self._time_steps[i + 1]),
                 self._obs_time,
             )
-            if self._obs_obj.obs == "iwc":
-                x_ind_no_rain = tl.get_1d_indices(
-                    (self._time_steps[i], self._time_steps[i + 1]),
-                    self._obs_time,
-                    mask=self._obs_obj.data["iwc_rain"][:],
-                )
             y_steps = tl.rebin_edges(self._model_height[i])
             for j in range(len(y_steps) - 1):
                 x_ind_adv = tl.get_adv_indices(
@@ -68,39 +62,18 @@ def _generate_downsample_product(self) -> None:
                 )
                 ind = np.outer(x_ind, y_ind)
                 ind_avd = np.outer(x_ind_adv, y_ind)
-                if self._obs_obj.obs == "cf":
+                if self._obs_obj.product == "cf":
                     data = self._reshape_data_to_window(ind, x_ind, y_ind)
                     if data is None:
                         continue
                     product_dict = self._regrid_cf(product_dict, i, j, data)
                     data_adv = self._reshape_data_to_window(ind_avd, x_ind_adv, y_ind)
-                    if data_adv is None:
-                        msg = "No data for advection"
-                        raise RuntimeError(msg)
-                    product_adv_dict = self._regrid_cf(product_adv_dict, i, j, data_adv)
-                elif self._obs_obj.obs == "iwc":
-                    x_ind_no_rain_adv = tl.get_adv_indices(
-                        model_t[i],
-                        self._time_adv[i, j],
-                        self._obs_time,
-                        mask=self._obs_obj.data["iwc_rain"][:],
-                    )
-                    ind_no_rain = np.outer(x_ind_no_rain, y_ind)
-                    ind_no_rain_adv = np.outer(x_ind_no_rain_adv, y_ind)
-                    product_dict = self._regrid_iwc(
-                        product_dict,
-                        i,
-                        j,
-                        ind,
-                        ind_no_rain,
-                    )
-                    product_adv_dict = self._regrid_iwc(
-                        product_adv_dict,
-                        i,
-                        j,
-                        ind_avd,
-                        ind_no_rain_adv,
-                    )
+                    # No advection window data (e.g. masked model wind collapses
+                    # the window); leave the cell unfilled like the standard grid.
+                    if data_adv is not None:
+                        product_adv_dict = self._regrid_cf(
+                            product_adv_dict, i, j, data_adv
+                        )
                 else:
                     product_dict = self._regrid_product(product_dict, i, j, ind)
                     product_adv_dict = self._regrid_product(
@@ -112,10 +85,8 @@ def _generate_downsample_product(self) -> None:
         self._append_data2object([product_dict, product_adv_dict])
 
     def _get_method_storage(self) -> tuple[dict, dict]:
-        if self._obs_obj.obs == "cf":
+        if self._obs_obj.product == "cf":
             return self._cf_method_storage()
-        if self._obs_obj.obs == "iwc":
-            return self._iwc_method_storage()
         return self._product_method_storage()
 
     def _cf_method_storage(self) -> tuple[dict, dict]:
@@ -129,23 +100,10 @@ def _cf_method_storage(self) -> tuple[dict, dict]:
         }
         return cf_dict, cf_adv_dict
 
-    def _iwc_method_storage(self) -> tuple[dict, dict]:
-        iwc_dict = {
-            "iwc": ma.zeros(self._model_height.shape),
-            "iwc_att": ma.zeros(self._model_height.shape),
-            "iwc_rain": ma.zeros(self._model_height.shape),
-        }
-        iwc_adv_dict = {
-            "iwc_adv": ma.zeros(self._model_height.shape),
-            "iwc_att_adv": ma.zeros(self._model_height.shape),
-            "iwc_rain_adv": ma.zeros(self._model_height.shape),
-        }
-        return iwc_dict, iwc_adv_dict
-
     def _product_method_storage(self) -> tuple[dict, dict]:
-        product_dict = {f"{self._obs_obj.obs}": ma.zeros(self._model_height.shape)}
+        product_dict = {f"{self._obs_obj.product}": ma.zeros(self._model_height.shape)}
         product_adv_dict = {
-            f"{self._obs_obj.obs}_adv": ma.zeros(self._model_height.shape),
+            f"{self._obs_obj.product}_adv": ma.zeros(self._model_height.shape),
         }
         return product_dict, product_adv_dict
 
@@ -168,36 +126,10 @@ def _regrid_cf(storage: dict, i: int, j: int, data: npt.NDArray) -> dict:
         for key, downsample in storage.items():
             downsample[i, j] = ma.mean(data_ma)
             if "_A" in key and not data_ma.mask.all():
-                downsample[i, j] = tl.average_column_sum(data_ma)
+                downsample[i, j] = tl.fraction_of_cloudy_columns(data_ma)
             storage[key] = downsample
         return storage
 
-    def _regrid_iwc(
-        self,
-        storage: dict,
-        i: int,
-        j: int,
-        ind_rain: npt.NDArray,
-        ind_no_rain: npt.NDArray,
-    ) -> dict:
-        """Calculates average iwc value for each grid point."""
-        for key, down_sample in storage.items():
-            down_sample[i, j] = ma.masked
-            if "rain" not in key:
-                no_rain_data = self._obs_data[ind_no_rain]
-                if ind_no_rain.any() and not no_rain_data.mask.all():
-                    down_sample[i, j] = ma.mean(no_rain_data)
-            if "rain" in key:
-                rain_data = self._obs_data[ind_rain]
-                if ind_rain.any() and not rain_data.mask.all():
-                    down_sample[i, j] = ma.mean(rain_data)
-            if "att" in key:
-                no_rain_att_data = self._obs_obj.data["iwc_att"][ind_no_rain]
-                if ind_no_rain.any() and not no_rain_att_data.mask.all():
-                    down_sample[i, j] = ma.mean(no_rain_att_data)
-            storage[key] = down_sample
-        return storage
-
     def _regrid_product(self, storage: dict, i: int, j: int, ind: npt.NDArray) -> dict:
         """Calculates average of standard product value for each grid point."""
         for key, down_sample in storage.items():
@@ -213,8 +145,14 @@ def _regrid_product(self, storage: dict, i: int, j: int, ind: npt.NDArray) -> di
     def _append_data2object(self, data_storage: list) -> None:
         for storage in data_storage:
             for key in storage:
-                down_sample = storage[key]
-                self.model_obj.append_data(
-                    down_sample,
-                    f"{key}_{self.model_obj.model}{self.model_obj.cycle}",
-                )
+                self.model_obj.append_data(storage[key], key)
+
+
+def downsample_to_model_grid(
+    model_obj: ModelManager, obs_obj: ObservationManager
+) -> None:
+    """Downsample the observation product onto the model grid.
+
+    The downsampled fields are added to `model_obj` in place.
+    """
+    ProductGrid(model_obj, obs_obj)
diff --git a/cloudnetpy/model_evaluation/products/model_products.py b/cloudnetpy/model_evaluation/products/model_products.py
index 69bf32a5..2e934763 100644
--- a/cloudnetpy/model_evaluation/products/model_products.py
+++ b/cloudnetpy/model_evaluation/products/model_products.py
@@ -1,16 +1,22 @@
-import importlib
 import logging
-import os.path
 from os import PathLike
 
 import numpy as np
 import numpy.typing as npt
 from numpy import ma
 
+from cloudnetpy.constants import RS
 from cloudnetpy.datasource import DataSource
-from cloudnetpy.model_evaluation.model_metadata import MODELS, VARIABLES
-from cloudnetpy.model_evaluation.utils import file_exists
-from cloudnetpy.utils import isscalar
+from cloudnetpy.exceptions import ModelDataError
+from cloudnetpy.model_evaluation.model_metadata import (
+    ALTITUDE_LIMIT,
+    COMMON_VARIABLES,
+    CYCLE_VARIABLES,
+    LEVEL_DIMENSION,
+    MODEL_PREFIX,
+    MODEL_VARIABLE_NAMES,
+    OPTIONAL_IWC_VARIABLES,
+)
 
 
 class ModelManager(DataSource):
@@ -19,16 +25,17 @@ class ModelManager(DataSource):
     Args:
         model_file (str): Path to source model file.
         model (str): Name of model
-        output_file (str): name of output file name and path to save data
         product (str): name of product to generate
 
     Notes:
-        For this class to work, needed information of model in use should be found in
-        model_metadata.py
+        Model files are expected to be harmonized, so all models share the
+        same variable names and units. The only model-specific quantity is the
+        number of vertical levels, which is derived from the model height and
+        the shared `ALTITUDE_LIMIT`.
 
-        Output_file is given for saving all cycles to same nc-file. Some variables
-        are same in control run and cycles so checking existence of output-file
-        prevents duplicates as well as unnecessary processing.
+        One :class:`ModelManager` (and one output file) corresponds to exactly
+        one model run. The model's own fields are stored with the `model_`
+        prefix; the model identity is kept in the file's global attributes.
 
         Class inherits DataSource interface from CloudnetPy.
     """
@@ -37,142 +44,147 @@ def __init__(
         self,
         model_file: str | PathLike,
         model: str,
-        output_file: str | PathLike,
         product: str,
-        *,
-        check_file: bool = True,
     ) -> None:
         super().__init__(model_file)
         self.model = model
-        self.model_info = MODELS[model]
-        self.model_vars = VARIABLES["variables"]
         self._product = product
         self.keys: dict = {}
-        self._is_file = file_exists(output_file) if check_file else False
-        self.cycle = self._read_cycle_name(model_file)
+        self._n_levels = self._read_number_of_levels()
         self._add_variables()
         self._generate_products()
         self.date: list = []
         self.wind = self._calculate_wind_speed()
         self.resolution_h = self._get_horizontal_resolution()
 
-    def _read_cycle_name(self, model_file: str | PathLike) -> str:
-        """Get cycle name from model_metadata.py for saving variable name(s)."""
-        basename = os.path.basename(model_file)
-        try:
-            cycles = self.model_info.cycle
-            if cycles is None:
-                return ""
-            cycles_split = [x.strip() for x in cycles.split(",")]
-            for cycle in cycles_split:
-                if cycle in basename:
-                    return f"_{cycle}"
-        except AttributeError:
-            return ""
-        return ""
+    def _read_number_of_levels(self) -> int:
+        """Number of vertical levels below the altitude limit.
+
+        Model heights are ground-first and increase with level index, and the
+        number of levels below the limit is constant in time, so a single level
+        count can be used to drop the unused upper levels for any model.
+        """
+        if "height" not in self.dataset.variables:
+            msg = (
+                f"Model '{self.model}' is missing the 'height' variable. "
+                "It needs to be added to the model file."
+            )
+            raise ModelDataError(msg)
+        height_var = self.dataset.variables["height"]
+        if LEVEL_DIMENSION not in height_var.dimensions:
+            msg = (
+                f"Model '{self.model}' height is missing the "
+                f"'{LEVEL_DIMENSION}' dimension. It needs to be fixed in the "
+                "model file."
+            )
+            raise ModelDataError(msg)
+        height = self.to_m(height_var)
+        below_limit = (
+            np.any(height < ALTITUDE_LIMIT, axis=0)
+            if height.ndim > 1
+            else height < ALTITUDE_LIMIT
+        )
+        return int(np.sum(below_limit))
 
     def _generate_products(self) -> None:
         """Process needed data of model to a ModelManager object."""
-        cls = importlib.import_module(__name__).ModelManager
+        methods = {"cf": self._get_cf, "iwc": self._get_iwc, "lwc": self._get_lwc}
         try:
-            name = f"_get_{self._product}"
-            getattr(cls, name)(self)
-        except AttributeError:
+            method = methods[self._product]
+        except KeyError:
             msg = f"Invalid product name: {self._product}"
             logging.exception(msg)
             raise
+        method()
 
     def _get_cf(self) -> None:
         """Collect cloud fraction straight from model file."""
-        cf_name = self.get_model_var_names(("cf",))[0]
-        cf = self.getvar(cf_name)
+        cf = self._getvar_checked("cf")
         cf = self.cut_off_extra_levels(cf)
         cf[cf < 0.05] = ma.masked
-        self.append_data(cf, f"{self.model}{self.cycle}_cf")
-        self.keys[self._product] = f"{self.model}{self.cycle}_cf"
+        self.append_data(cf, f"{MODEL_PREFIX}cf")
+        self.keys[self._product] = f"{MODEL_PREFIX}cf"
 
     def _get_iwc(self) -> None:
         iwc = self.get_water_content("iwc")
         iwc[iwc < 1e-7] = ma.masked
-        self.append_data(iwc, f"{self.model}{self.cycle}_iwc")
-        self.keys[self._product] = f"{self.model}{self.cycle}_iwc"
+        self.append_data(iwc, f"{MODEL_PREFIX}iwc")
+        self.keys[self._product] = f"{MODEL_PREFIX}iwc"
 
     def _get_lwc(self) -> None:
         lwc = self.get_water_content("lwc")
         lwc[lwc < 1e-5] = ma.masked
-        self.append_data(lwc, f"{self.model}{self.cycle}_lwc")
-        self.keys[self._product] = f"{self.model}{self.cycle}_lwc"
+        self.append_data(lwc, f"{MODEL_PREFIX}lwc")
+        self.keys[self._product] = f"{MODEL_PREFIX}lwc"
 
     @staticmethod
     def get_model_var_names(args: tuple) -> list:
-        var = []
-        for arg in args:
-            var.append(VARIABLES[arg].long_name)
-        return var
+        return [MODEL_VARIABLE_NAMES[arg] for arg in args]
+
+    def _getvar_checked(self, *internal_keys: str) -> npt.NDArray:
+        """Fetch a model variable, raising a clear error if it is missing."""
+        names = [MODEL_VARIABLE_NAMES[key] for key in internal_keys]
+        try:
+            return self.getvar(*names)
+        except KeyError as err:
+            msg = (
+                f"Model '{self.model}' is missing variable '{names[0]}' "
+                f"required for product '{self._product}'."
+            )
+            raise ModelDataError(msg) from err
 
     def get_water_content(self, var: str) -> npt.NDArray:
-        p_name = self.get_model_var_names(("p",))[0]
-        t_name = self.get_model_var_names(("T",))[0]
-        lwc_name = self.get_model_var_names((var,))[0]
-        p = self.getvar(p_name)
-        t = self.getvar(t_name)
-        q = self.getvar(lwc_name)
+        p = self._getvar_checked("p")
+        t = self._getvar_checked("T")
+        q = self._get_mixing_ratio(var)
         wc = self._calc_water_content(q, p, t)
         wc = self.cut_off_extra_levels(wc)
         wc[wc < 0.0] = ma.masked
         return wc
 
+    def _get_mixing_ratio(self, var: str) -> npt.NDArray:
+        """Mixing ratio for a water-content product.
+
+        For ice (`iwc`) the total frozen-condensate mixing ratio is returned:
+        the required cloud-ice field plus any snow/graupel categories present in
+        the model file. This keeps the model IWC consistent with the observed
+        IWC, which also includes precipitating ice.
+        """
+        q = self._getvar_checked(var)
+        if var == "iwc":
+            for name in OPTIONAL_IWC_VARIABLES:
+                if name in self.dataset.variables:
+                    q = q + self.getvar(name)
+        return q
+
     @staticmethod
     def _calc_water_content(
         q: npt.NDArray, p: npt.NDArray, t: npt.NDArray
     ) -> npt.NDArray:
-        return q * p / (287 * t)
+        return q * p / (RS * t)
 
     def _add_variables(self) -> None:
-        """Add basic variables off model and cycle."""
-
-        def _add_common_variables() -> None:
-            """Model variables that are always the same within cycles."""
-            wanted_vars = self.model_vars.common_var
-            if wanted_vars is None:
-                msg = f"Model {self.model} has no common variables"
-                raise ValueError(msg)
-            wanted_vars_split = [x.strip() for x in wanted_vars.split(",")]
-            for var in wanted_vars_split:
-                if var in self.dataset.variables:
-                    data = self.dataset.variables[var][:]
-                    if not isscalar(data) and len(data) > 25:
-                        data = self.cut_off_extra_levels(self.dataset.variables[var][:])
-                    self.append_data(data, f"{var}")
-
-        def _add_cycle_variables() -> None:
-            """Add cycle depending variables."""
-            wanted_vars = self.model_vars.cycle_var
-            if wanted_vars is None:
-                msg = f"Model {self.model} has no cycle variables"
-                raise ValueError(msg)
-            wanted_vars_split = [x.strip() for x in wanted_vars.split(",")]
-            for var in wanted_vars_split:
-                if var in self.dataset.variables:
-                    data = self.dataset.variables[var][:]
-                    if data.ndim > 1 or len(data) > 25:
-                        data = self.cut_off_extra_levels(self.dataset.variables[var][:])
-                    self.append_data(data, f"{self.model}{self.cycle}_{var}")
-                if var == "height":
-                    self.keys["height"] = f"{self.model}{self.cycle}_{var}"
-
-        if not self._is_file:
-            _add_common_variables()
-        _add_cycle_variables()
+        """Add common coordinate variables and the model's own fields."""
+
+        def _add_variable(var: str, key: str) -> None:
+            ncvar = self.dataset.variables[var]
+            data = ncvar[:]
+            if LEVEL_DIMENSION in ncvar.dimensions:
+                data = self.cut_off_extra_levels(data)
+            self.append_data(data, key)
+
+        for var in COMMON_VARIABLES:
+            if var in self.dataset.variables:
+                _add_variable(var, var)
+        for var in CYCLE_VARIABLES:
+            if var in self.dataset.variables:
+                _add_variable(var, f"{MODEL_PREFIX}{var}")
+            if var == "height":
+                self.keys["height"] = f"{MODEL_PREFIX}{var}"
 
     def cut_off_extra_levels(self, data: npt.NDArray) -> npt.NDArray:
-        """Remove unused levels (over 22km) from model data."""
-        try:
-            level = self.model_info.level
-        except KeyError:
-            return data
-
-        return data[:, :level] if data.ndim > 1 else data[:level]
+        """Remove unused levels (above the altitude limit) from model data."""
+        return data[:, : self._n_levels] if data.ndim > 1 else data[: self._n_levels]
 
     def _calculate_wind_speed(self) -> npt.NDArray:
         """Real wind from x- and y-components."""
@@ -183,5 +195,20 @@ def _calculate_wind_speed(self) -> npt.NDArray:
         return np.sqrt(ma.power(u.data, 2) + ma.power(v.data, 2))
 
     def _get_horizontal_resolution(self) -> float:
-        h_res = self.getvar("horizontal_resolution")
-        return float(np.unique(h_res.data)[0])
+        try:
+            h_res = self.getvar("horizontal_resolution")
+        except KeyError as err:
+            msg = (
+                f"Model '{self.model}' is missing 'horizontal_resolution'. "
+                "It needs to be added to the model file."
+            )
+            raise ModelDataError(msg) from err
+        unique = np.unique(ma.masked_invalid(h_res).compressed())
+        if unique.size != 1 or unique[0] <= 0:
+            msg = (
+                f"Model '{self.model}' has invalid horizontal_resolution "
+                f"({unique.tolist()}); expected a single positive value. "
+                "It needs to be fixed in the model file."
+            )
+            raise ModelDataError(msg)
+        return float(unique[0])
diff --git a/cloudnetpy/model_evaluation/products/observation_products.py b/cloudnetpy/model_evaluation/products/observation_products.py
index fbf8fd87..92ec8489 100644
--- a/cloudnetpy/model_evaluation/products/observation_products.py
+++ b/cloudnetpy/model_evaluation/products/observation_products.py
@@ -6,7 +6,6 @@
 import numpy.typing as npt
 from numpy import ma
 
-from cloudnetpy import utils
 from cloudnetpy.datasource import DataSource
 from cloudnetpy.products.product_tools import CategorizeBits, CategoryBits
 
@@ -15,7 +14,7 @@ class ObservationManager(DataSource):
     """Class to collect and manage observations for downsampling.
 
     Args:
-        obs (str): Name of observation product
+        product (str): Name of observation product
         obs_file (str): Path to source observation file
 
     Notes:
@@ -26,9 +25,9 @@ class ObservationManager(DataSource):
         should be processed using CloudnetPy for this class to work properly.
     """
 
-    def __init__(self, obs: str, obs_file: str | PathLike) -> None:
+    def __init__(self, product: str, obs_file: str | PathLike) -> None:
         super().__init__(obs_file)
-        self.obs = obs
+        self.product = product
         self._file = obs_file
         self.date = self._get_date()
         self.radar_freq = self._get_radar_frequency()
@@ -62,15 +61,15 @@ def _get_z_sensitivity(self) -> npt.NDArray | None:
     def _generate_product(self) -> None:
         """Process needed data of observation to a ObservationManager object."""
         try:
-            if self.obs == "cf":
+            if self.product == "cf":
                 self.append_data(self._generate_cf(), "cf")
             else:
-                self.append_data(self.getvar(self.obs), self.obs)
-                if self.obs == "iwc":
-                    self._generate_iwc_masks()
+                self.append_data(self.getvar(self.product), self.product)
+                if self.product == "iwc":
+                    self._mask_iwc()
             self.append_data(self.getvar("height"), "height")
         except (KeyError, RuntimeError):
-            msg = f"Failed to read {self.obs} from {self._file}"
+            msg = f"Failed to read {self.product} from {self._file}"
             logging.exception(msg)
             raise
 
@@ -100,55 +99,14 @@ def _mask_cloud_bits(cloud_mask: npt.NDArray) -> npt.NDArray:
             cloud_mask[cloud_mask == i] = 0
         return cloud_mask
 
-    def _check_rainrate(self) -> bool:
-        """Check if rainrate in file."""
-        try:
-            self.getvar("rainrate")
-        except KeyError:
-            return False
-        return True
-
-    def _get_rainrate_threshold(self) -> int:
-        if self.radar_freq is None:
-            msg = "Radar frequency not found from file"
-            raise RuntimeError(msg)
-        wband = utils.get_wl_band(float(self.radar_freq))
-        rainrate_threshold = 8
-        if wband == "W":
-            rainrate_threshold = 2
-        return rainrate_threshold
-
-    def _rain_index(self) -> npt.NDArray:
-        rainrate = self.getvar("rainrate")
-        rainrate_threshold = self._get_rainrate_threshold()
-        return rainrate > rainrate_threshold
-
-    def _generate_iwc_masks(self) -> None:
-        """Generates ice water content variables with different masks."""
-        # TODO: Differences with CloudnetPy (status=2) and Legacy data (status=3)
-        iwc = self.getvar(self.obs)
-        iwc_status = self.getvar("iwc_retrieval_status")
-        self._mask_iwc_att(iwc, iwc_status)
-        self._get_rain_iwc(iwc_status)
-        self._mask_iwc(iwc, iwc_status)
-
-    def _mask_iwc(self, iwc: npt.NDArray, iwc_status: npt.NDArray) -> None:
-        """Leaves only reliable data and corrected liquid attenuation."""
-        iwc_mask = ma.copy(iwc)
-        iwc_mask[np.bitwise_and(iwc_status != 1, iwc_status != 2)] = ma.masked
-        self.append_data(iwc_mask, "iwc")
-
-    def _mask_iwc_att(self, iwc: npt.NDArray, iwc_status: npt.NDArray) -> None:
-        """Leaves only where reliable data, corrected liquid attenuation
-        and uncorrected liquid attenuation.
+    def _mask_iwc(self) -> None:
+        """Keeps only reliable ice water content retrievals.
+
+        Status 1 is a reliable retrieval and status 3 is a retrieval with the
+        radar corrected for liquid, rain and melting attenuation; everything
+        else (uncorrected attenuation, lidar-only, rain) is masked out.
         """
-        iwc_att = ma.copy(iwc)
-        iwc_att[iwc_status > 3] = ma.masked
-        self.append_data(iwc_att, "iwc_att")
-
-    def _get_rain_iwc(self, iwc_status: npt.NDArray) -> None:
-        """Finds columns where is rain, return boolean of x-axis shape."""
-        iwc_rain = np.zeros(iwc_status.shape, dtype=bool)
-        iwc_rain[iwc_status == 5] = 1
-        iwc_rain = np.any(iwc_rain, axis=1)
-        self.append_data(iwc_rain, "iwc_rain")
+        iwc = self.getvar("iwc")
+        iwc_status = self.getvar("iwc_retrieval_status")
+        iwc[~np.isin(iwc_status, (1, 3))] = ma.masked
+        self.append_data(iwc, "iwc")
diff --git a/cloudnetpy/model_evaluation/products/product_resampling.py b/cloudnetpy/model_evaluation/products/product_resampling.py
index eba6513c..f2b0f3b7 100644
--- a/cloudnetpy/model_evaluation/products/product_resampling.py
+++ b/cloudnetpy/model_evaluation/products/product_resampling.py
@@ -5,12 +5,11 @@
 import cloudnetpy.model_evaluation.products.tools as tl
 from cloudnetpy.model_evaluation.file_handler import (
     add_time_attribute,
-    add_var2ncfile,
     save_downsampled_file,
     update_attributes,
 )
-from cloudnetpy.model_evaluation.products.advance_methods import AdvanceProductMethods
-from cloudnetpy.model_evaluation.products.grid_methods import ProductGrid
+from cloudnetpy.model_evaluation.products.advance_methods import apply_cirrus_filter
+from cloudnetpy.model_evaluation.products.grid_methods import downsample_to_model_grid
 from cloudnetpy.model_evaluation.products.model_products import ModelManager
 from cloudnetpy.model_evaluation.products.observation_products import ObservationManager
 from cloudnetpy.model_evaluation.utils import file_exists
@@ -19,80 +18,75 @@
 
 def process_L3_day_product(
     model: str,
-    obs: str,
-    model_files: list[str | PathLike],
+    product: str,
+    model_file: str | PathLike,
     product_file: str | PathLike,
     output_file: str | PathLike,
     uuid: str | UUID | None = None,
+    model_name: str | None = None,
+    site_name: str | None = None,
     *,
     overwrite: bool = False,
 ) -> UUID:
-    """Main function to generate downsample of observations to match model grid.
+    """Generate a downsampling of observations to match the model grid.
 
-    This function will generate a L3 product nc-file. It includes the information of
-    downsampled observation products for each model cycles and model products
-    and other variables of each cycles.
+    Generates an L3 product nc-file for a single model run. It contains the
+    model's own fields (prefixed with ``model_``) and the observation products
+    downsampled onto the model grid. The model identity is stored in the file's
+    global attributes, not in the variable names.
 
     Args:
         model (str): Name of model
-        obs (str): Name of product to generate
-        model_files (list): List of model + cycles file path(s) to be generated
+        product (str): Name of product to generate
+        model_file (str): Model file path
         product_file (str): Source file path of L2 observation product
         output_file (str): Path and name of L3 day scale product output file
-        keep_uuid (bool): If True, keeps the UUID of the old file, if that exists.
-                          Default is False when new UUID is generated.
         uuid (str): Set specific UUID for the file.
-        overwrite (bool): If file exists, but still want to recreate it then True,
-                          default False
+        model_name (str): Human-readable model name for plot titles. Falls back
+                       to the model id when not given.
+        site_name (str): Human-readable site name for the location attribute and
+                       plot subtitle. Falls back to the source file's location
+                       when not given.
+        overwrite (bool): Recreate the output file if it already exists.
 
     Raises:
         RuntimeError: Failed to create the L3 product file.
         ValueError (Warning): No ice clouds in model data
-
-    Notes:
-        Model file(s) are given as a list to make all different cycles to be at same
-        nc-file. If list includes more than one model file, nc-file is created within
-        the first round. With rest of rounds, downsample observation and model data
-        is added to same L3 day nc-file.
+        FileExistsError: Output file exists and overwrite is False.
 
     Examples:
         >>> from cloudnetpy.model_evaluation.products.product_resampling import \
         process_L3_day_product
-        >>> product = 'cf'
-        >>> model = 'ecmwf'
-        >>> model_file = 'ecmwf.nc'
-        >>> input_file = 220190517_mace-head_categorize.nchead_categorize.nc
-        >>> output_file = 'cf_ecmwf.nc'
-        >>> process_L3_day_product(model, product, [model_file], input_file,
-        output_file)
+        >>> process_L3_day_product('ecmwf', 'cf', 'ecmwf.nc',
+        ... 'categorize.nc', 'l3-cf.nc')
     """
     uuid = get_uuid(uuid)
-    product_obj = ObservationManager(obs, product_file)
-    tl.check_model_file_list(model, model_files)
-    for m_file in model_files:
-        model_obj = ModelManager(
-            m_file,
-            model,
-            output_file,
-            obs,
-            check_file=not overwrite,
-        )
+    if file_exists(output_file) and not overwrite:
+        msg = f"Output file {output_file} exists, use overwrite=True to recreate"
+        raise FileExistsError(msg)
+    product_obj = ObservationManager(product, product_file)
+    model_obj = ModelManager(model_file, model, product)
+    try:
         try:
-            AdvanceProductMethods(model_obj, m_file, product_obj)
+            apply_cirrus_filter(model_obj, product_obj)
         except ValueError as e:
             logging.info(e)
-        ProductGrid(model_obj, product_obj)
+        downsample_to_model_grid(model_obj, product_obj)
         attributes = add_time_attribute(product_obj.date)
         update_attributes(model_obj.data, attributes)
-        if not file_exists(output_file) or overwrite:
-            tl.add_date(model_obj, product_obj)
-            save_downsampled_file(
-                f"{obs}_{model}",
-                output_file,
-                (model_obj, product_obj),
-                (model_files, product_file),
-                uuid,
-            )
-        else:
-            add_var2ncfile(model_obj, output_file)
+        tl.add_date(model_obj, product_obj)
+        save_downsampled_file(
+            product,
+            output_file,
+            model_obj,
+            product_obj,
+            [model_file],
+            product_file,
+            uuid,
+            model_name,
+            site_name,
+        )
+    finally:
+        model_obj.close()
+        product_obj.close()
     return uuid
diff --git a/cloudnetpy/model_evaluation/products/tools.py b/cloudnetpy/model_evaluation/products/tools.py
index ef4d1463..2aab2dca 100644
--- a/cloudnetpy/model_evaluation/products/tools.py
+++ b/cloudnetpy/model_evaluation/products/tools.py
@@ -1,9 +1,5 @@
 import datetime
-import logging
-import os.path
-from collections.abc import Sequence
 from datetime import timedelta
-from os import PathLike
 
 import numpy as np
 import numpy.typing as npt
@@ -13,15 +9,6 @@
 from cloudnetpy.model_evaluation.products.observation_products import ObservationManager
 
 
-def check_model_file_list(name: str, models: Sequence[str | PathLike]) -> None:
-    """Check that files in models are from same model and date."""
-    for m in models:
-        if name not in os.path.basename(m):
-            logging.error("Invalid model file set")
-            msg = f"{m} not from {name}"
-            raise AttributeError(msg)
-
-
 def time2datetime(time: npt.NDArray, date: datetime.datetime) -> npt.NDArray:
     return np.asarray([date + timedelta(hours=float(t)) for t in time])
 
@@ -29,13 +16,13 @@ def time2datetime(time: npt.NDArray, date: datetime.datetime) -> npt.NDArray:
 def rebin_edges(arr: npt.NDArray) -> npt.NDArray:
     """Rebins array bins by half and adds boundaries."""
     new_arr = [(arr[i] + arr[i + 1]) / 2 for i in range(len(arr) - 1)]
-    new_arr.insert(0, arr[0] - ((arr[0] + arr[1]) / 2))
+    new_arr.insert(0, arr[0] - (arr[1] - arr[0]))
     new_arr.insert(len(new_arr), arr[-1] + (arr[-1] - arr[-2]))
     return np.array(new_arr)
 
 
 def calculate_advection_time(
-    resolution: int,
+    resolution: float,
     wind: ma.MaskedArray,
     sampling: int,
 ) -> npt.NDArray:
@@ -81,7 +68,7 @@ def get_obs_window_size(ind_x: npt.NDArray, ind_y: npt.NDArray) -> tuple | None:
     """Returns shape (tuple) of window area, where values are True."""
     x = np.where(ind_x)[0]
     y = np.where(ind_y)[0]
-    if np.any(x) and np.any(y):
+    if x.size > 0 and y.size > 0:
         return x[-1] - x[0] + 1, y[-1] - y[0] + 1
     return None
 
@@ -91,6 +78,13 @@ def add_date(model_obj: ModelManager, obs_obj: ObservationManager) -> None:
         model_obj.date.append(getattr(obs_obj.dataset, a))
 
 
-def average_column_sum(data: npt.NDArray) -> npt.NDArray:
-    """Returns average sum of columns which have any data."""
-    return np.nanmean(np.nansum(data, 1) > 0)
+def fraction_of_cloudy_columns(data: ma.MaskedArray) -> np.float64:
+    """Fraction of columns (profiles) containing cloud at any height.
+
+    This is the "by area" cloud fraction: a column counts as cloudy if it has
+    cloud anywhere in the vertical, regardless of cloud depth. Masked
+    observation pixels are ignored, consistent with the "by volume" cloud
+    fraction (``ma.mean``).
+    """
+    column_has_cloud = ma.sum(data, axis=1) > 0
+    return ma.mean(column_has_cloud)
diff --git a/cloudnetpy/model_evaluation/statistics/statistical_methods.py b/cloudnetpy/model_evaluation/statistics/statistical_methods.py
index ab00c5c7..e29e1e33 100644
--- a/cloudnetpy/model_evaluation/statistics/statistical_methods.py
+++ b/cloudnetpy/model_evaluation/statistics/statistical_methods.py
@@ -1,13 +1,10 @@
 import logging
-import os
-import sys
+from collections.abc import Callable
 
 import numpy as np
 import numpy.typing as npt
 from numpy import ma
 
-sys.path.append(os.path.dirname(os.path.abspath(__file__)))
-
 
 class DayStatistics:
     """Class for calculating statistical analysis of day scale products.
@@ -54,33 +51,26 @@ def __init__(
         self.observation_data = observation
         self._generate_day_statistics()
 
-    def _get_method_attr(self) -> tuple[str, tuple]:
-        full_name = ""
+    def _get_method_attr(self) -> tuple[Callable, tuple]:
         params = (self.model_data, self.observation_data)
         if self.method == "error":
-            full_name = "relative_error"
+            return relative_error, params
         if self.method == "aerror":
-            full_name = "absolute_error"
+            return absolute_error, params
         if self.method == "area":
-            full_name = "calc_common_area_sum"
+            return calc_common_area_sum, params
         if self.method == "hist":
-            return "histogram", (
-                self.product,
-                self.model_data,
-                self.observation_data,
-            )
+            return histogram, (self.product, *params)
         if self.method == "vertical":
-            full_name = "vertical_profile"
-        return full_name, params
+            return vertical_profile, params
+        msg = f"Unknown statistical method: {self.method}"
+        raise RuntimeError(msg)
 
     def _generate_day_statistics(self) -> None:
-        full_name, params = self._get_method_attr()
-        cls = __import__("statistical_methods")
+        method, params = self._get_method_attr()
         try:
-            self.model_stat, self.observation_stat = getattr(cls, f"{full_name}")(
-                *params,
-            )
-            self.title = cls.day_stat_title(self.method, self.product)
+            self.model_stat, self.observation_stat = method(*params)
+            self.title = day_stat_title(self.method, self.product)
         except RuntimeError:
             msg = f"Failed to calculate {self.method} of {self.product[0]}"
             logging.exception(msg)
@@ -91,7 +81,10 @@ def relative_error(
     observation: ma.MaskedArray,
 ) -> tuple[float, str]:
     model, observation = combine_masked_indices(model, observation)
-    error = ((model - observation) / observation) * 100
+    # Very small observation values make the relative error overflow; the
+    # resulting huge values are expected, so ignore the warning.
+    with np.errstate(divide="ignore", over="ignore"):
+        error = ((model - observation) / observation) * 100
     return np.round(error, 2), ""
 
 
@@ -129,7 +122,10 @@ def _indices_of_mask_sum() -> float:
         model[model < np.min(observation)] = ma.masked
         combine_mask = model.mask + observation.mask
         common_mask = np.bitwise_and(model.mask, observation.mask)
-        the_match = np.sum(~combine_mask) / np.sum(~common_mask) * 100
+        valid_area = np.sum(~common_mask)
+        if valid_area == 0:
+            return 0.0
+        the_match = np.sum(~combine_mask) / valid_area * 100
         return np.round(the_match, 2)
 
     match = _indices_of_mask_sum()
@@ -160,9 +156,8 @@ def histogram(
 
 
 def vertical_profile(model: ma.MaskedArray, observation: ma.MaskedArray) -> tuple:
-    if model.shape[0] > 25:
-        model = model.T
-        observation = observation.T
+    # L3 fields are shaped (time, level); average over time (axis 0) to get the
+    # vertical profile over height.
     model_vertical = ma.mean(model, axis=0)
     obs_vertical = np.nanmean(observation, axis=0)
     return model_vertical, obs_vertical
diff --git a/cloudnetpy/model_evaluation/tests/e2e/process_cf/main.py b/cloudnetpy/model_evaluation/tests/e2e/process_cf/main.py
index 7a7effcb..90037e75 100755
--- a/cloudnetpy/model_evaluation/tests/e2e/process_cf/main.py
+++ b/cloudnetpy/model_evaluation/tests/e2e/process_cf/main.py
@@ -21,7 +21,7 @@ def _process() -> None:
     product_resampling.process_L3_day_product(
         "ecmwf",
         "cf",
-        [test_file_model],
+        test_file_model,
         test_file_product,
         temp_file,
     )
diff --git a/cloudnetpy/model_evaluation/tests/e2e/process_cf/tests.py b/cloudnetpy/model_evaluation/tests/e2e/process_cf/tests.py
index d0010c6a..4e8e65da 100644
--- a/cloudnetpy/model_evaluation/tests/e2e/process_cf/tests.py
+++ b/cloudnetpy/model_evaluation/tests/e2e/process_cf/tests.py
@@ -16,19 +16,18 @@ def test_that_has_correct_attributes(self) -> None:
         assert nc.year == "2019"
         assert nc.month == "05"
         assert nc.day == "17"
-        assert nc.title == "Downsampled Cf of ecmwf from Mace-Head"
+        assert nc.title == "Observed and modeled cloud fraction over Mace-Head"
         assert nc.cloudnet_file_type == "l3-cf"
         assert nc.Conventions == "CF-1.8"
         assert (
-            nc.source
-            == "Observation file: 20190517_mace-head_categorize.nc\necmwf file(s): 20190517_mace-head_ecmwf.nc"
+            nc.source == "20190517_mace-head_categorize.nc\n20190517_mace-head_ecmwf.nc"
         )
         nc.close()
 
     @pytest.mark.reprocess()
     @pytest.mark.parametrize(
         "key",
-        ["cf_V_ecmwf", "cf_A_ecmwf", "cf_V_adv_ecmwf", "cf_A_adv_ecmwf"],
+        ["cf_V", "cf_A", "cf_V_adv", "cf_A_adv"],
     )
     def test_that_has_correct_product_variables(self, key) -> None:
         nc = netCDF4.Dataset(self.full_path)
@@ -48,7 +47,7 @@ def test_that_has_correct_model_variables(self, key) -> None:
     @pytest.mark.reprocess()
     @pytest.mark.parametrize(
         "key",
-        ["ecmwf_forecast_time", "ecmwf_height", "ecmwf_cf", "ecmwf_cf_cirrus"],
+        ["model_forecast_time", "model_height", "model_cf", "model_cf_cirrus"],
     )
     def test_that_has_correct_cycle_variables(self, key) -> None:
         nc = netCDF4.Dataset(self.full_path)
diff --git a/cloudnetpy/model_evaluation/tests/e2e/process_iwc/main.py b/cloudnetpy/model_evaluation/tests/e2e/process_iwc/main.py
index 6b57a8b8..5418e609 100755
--- a/cloudnetpy/model_evaluation/tests/e2e/process_iwc/main.py
+++ b/cloudnetpy/model_evaluation/tests/e2e/process_iwc/main.py
@@ -25,7 +25,7 @@ def _process() -> None:
     product_resampling.process_L3_day_product(
         "ecmwf",
         "iwc",
-        [test_file_model],
+        test_file_model,
         test_file_product,
         temp_file,
     )
diff --git a/cloudnetpy/model_evaluation/tests/e2e/process_iwc/tests.py b/cloudnetpy/model_evaluation/tests/e2e/process_iwc/tests.py
index bf87fca0..3c05afd0 100644
--- a/cloudnetpy/model_evaluation/tests/e2e/process_iwc/tests.py
+++ b/cloudnetpy/model_evaluation/tests/e2e/process_iwc/tests.py
@@ -16,12 +16,12 @@ def test_that_has_correct_attributes(self) -> None:
         assert nc.year == "2019"
         assert nc.month == "05"
         assert nc.day == "17"
-        assert nc.title == "Downsampled Iwc of ecmwf from Mace-Head"
+        assert nc.title == "Observed and modeled ice water content over Mace-Head"
         assert nc.cloudnet_file_type == "l3-iwc"
         assert nc.Conventions == "CF-1.8"
         assert (
             nc.source
-            == "Observation file: 20190517_mace-head_iwc-Z-T-method.nc\necmwf file(s): 20190517_mace-head_ecmwf.nc"
+            == "20190517_mace-head_iwc-Z-T-method.nc\n20190517_mace-head_ecmwf.nc"
         )
         nc.close()
 
@@ -29,12 +29,8 @@ def test_that_has_correct_attributes(self) -> None:
     @pytest.mark.parametrize(
         "key",
         [
-            "iwc_ecmwf",
-            "iwc_att_ecmwf",
-            "iwc_rain_ecmwf",
-            "iwc_adv_ecmwf",
-            "iwc_att_adv_ecmwf",
-            "iwc_rain_adv_ecmwf",
+            "iwc",
+            "iwc_adv",
         ],
     )
     def test_that_has_correct_product_variables(self, key) -> None:
@@ -55,7 +51,7 @@ def test_that_has_correct_model_variables(self, key) -> None:
     @pytest.mark.reprocess()
     @pytest.mark.parametrize(
         "key",
-        ["ecmwf_forecast_time", "ecmwf_height", "ecmwf_iwc"],
+        ["model_forecast_time", "model_height", "model_iwc"],
     )
     def test_that_has_correct_cycle_variables(self, key) -> None:
         nc = netCDF4.Dataset(self.full_path)
diff --git a/cloudnetpy/model_evaluation/tests/e2e/process_lwc/main.py b/cloudnetpy/model_evaluation/tests/e2e/process_lwc/main.py
index f8705695..25f3e651 100755
--- a/cloudnetpy/model_evaluation/tests/e2e/process_lwc/main.py
+++ b/cloudnetpy/model_evaluation/tests/e2e/process_lwc/main.py
@@ -21,7 +21,7 @@ def _process() -> None:
     product_resampling.process_L3_day_product(
         "ecmwf",
         "lwc",
-        [test_file_model],
+        test_file_model,
         test_file_product,
         temp_file,
     )
diff --git a/cloudnetpy/model_evaluation/tests/e2e/process_lwc/tests.py b/cloudnetpy/model_evaluation/tests/e2e/process_lwc/tests.py
index 8f1c3bde..c54d9340 100644
--- a/cloudnetpy/model_evaluation/tests/e2e/process_lwc/tests.py
+++ b/cloudnetpy/model_evaluation/tests/e2e/process_lwc/tests.py
@@ -16,17 +16,17 @@ def test_that_has_correct_attributes(self) -> None:
         assert nc.year == "2019"
         assert nc.month == "05"
         assert nc.day == "17"
-        assert nc.title == "Downsampled Lwc of ecmwf from Mace-Head"
+        assert nc.title == "Observed and modeled liquid water content over Mace-Head"
         assert nc.cloudnet_file_type == "l3-lwc"
         assert nc.Conventions == "CF-1.8"
         assert (
             nc.source
-            == "Observation file: 20190517_mace-head_lwc-scaled-adiabatic.nc\necmwf file(s): 20190517_mace-head_ecmwf.nc"
+            == "20190517_mace-head_lwc-scaled-adiabatic.nc\n20190517_mace-head_ecmwf.nc"
         )
         nc.close()
 
     @pytest.mark.reprocess()
-    @pytest.mark.parametrize("key", ["lwc_ecmwf", "lwc_adv_ecmwf"])
+    @pytest.mark.parametrize("key", ["lwc", "lwc_adv"])
     def test_that_has_correct_product_variables(self, key) -> None:
         nc = netCDF4.Dataset(self.full_path)
         assert key in nc.variables
@@ -45,7 +45,7 @@ def test_that_has_correct_model_variables(self, key) -> None:
     @pytest.mark.reprocess()
     @pytest.mark.parametrize(
         "key",
-        ["ecmwf_forecast_time", "ecmwf_height", "ecmwf_lwc"],
+        ["model_forecast_time", "model_height", "model_lwc"],
     )
     def test_that_has_correct_cycle_variables(self, key) -> None:
         nc = netCDF4.Dataset(self.full_path)
diff --git a/cloudnetpy/model_evaluation/tests/unit/conftest.py b/cloudnetpy/model_evaluation/tests/unit/conftest.py
index 7e0c654e..afd25a54 100644
--- a/cloudnetpy/model_evaluation/tests/unit/conftest.py
+++ b/cloudnetpy/model_evaluation/tests/unit/conftest.py
@@ -61,6 +61,117 @@ def model_file(tmpdir_factory, file_metadata) -> str:
     return file_name
 
 
+def _write_model_file(
+    file_name,
+    file_metadata,
+    *,
+    n_levels: int,
+    heights: list,
+    horizontal_resolution: bool = True,
+    resolution: float = 9,
+    clouds: bool = True,
+    level_var: bool = True,
+) -> str:
+    """Writes a harmonized model file, optionally omitting certain variables.
+
+    Used to emulate the structural differences between models (e.g. arpege and
+    icon lack horizontal_resolution and a level coordinate; ecmwf-open lacks
+    cloud variables).
+    """
+    root_grp = netCDF4.Dataset(file_name, "w", format="NETCDF4_CLASSIC")
+    time = 3
+    root_grp.createDimension("time", time)
+    root_grp.createDimension("level", n_levels)
+    _create_global_attributes(root_grp, file_metadata)
+    root_grp.source = "Test Model (TM)"
+    var = root_grp.createVariable("time", "f8", "time")
+    var[:] = ma.array([2, 6, 10])
+    if level_var:
+        var = root_grp.createVariable("level", "f8", "level")
+        var[:] = np.arange(n_levels)
+    var = root_grp.createVariable("latitude", "f8")
+    var[:] = 1
+    var = root_grp.createVariable("longitude", "f8")
+    var[:] = 1
+    if horizontal_resolution:
+        var = root_grp.createVariable("horizontal_resolution", "f8")
+        var[:] = resolution
+    var = root_grp.createVariable("height", "f8", ("time", "level"))
+    var[:] = ma.array([heights, heights, heights])
+    var.units = "m"
+    var = root_grp.createVariable("forecast_time", "f8", "time")
+    var[:] = ma.array([1, 5, 10])
+    shape = (time, n_levels)
+    if clouds:
+        for name in ("cloud_fraction", "qi", "ql"):
+            var = root_grp.createVariable(name, "f8", ("time", "level"))
+            var[:] = ma.ones(shape) * 0.05
+    var = root_grp.createVariable("temperature", "f8", ("time", "level"))
+    var[:] = ma.ones(shape) * 300
+    var = root_grp.createVariable("pressure", "f8", ("time", "level"))
+    var[:] = ma.ones(shape) * 1000
+    var = root_grp.createVariable("uwind", "f8", ("time", "level"))
+    var[:] = ma.ones(shape) * 2
+    var = root_grp.createVariable("vwind", "f8", ("time", "level"))
+    var[:] = ma.ones(shape) * 3
+    root_grp.close()
+    return file_name
+
+
+@pytest.fixture(scope="session")
+def model_file_tall(tmpdir_factory, file_metadata) -> str:
+    """Model with 5 levels, of which only the lowest 3 are below 22 km."""
+    file_name = tmpdir_factory.mktemp("data").join("tall.nc")
+    return _write_model_file(
+        file_name,
+        file_metadata,
+        n_levels=5,
+        heights=[100, 5000, 15000, 25000, 40000],
+    )
+
+
+@pytest.fixture(scope="session")
+def model_file_no_hres(tmpdir_factory, file_metadata) -> str:
+    """Model without horizontal_resolution and level coordinate (icon-like)."""
+    file_name = tmpdir_factory.mktemp("data").join("no_hres.nc")
+    return _write_model_file(
+        file_name,
+        file_metadata,
+        n_levels=2,
+        heights=[100, 5000],
+        horizontal_resolution=False,
+        level_var=False,
+    )
+
+
+@pytest.fixture(scope="session")
+def model_file_zero_hres(tmpdir_factory, file_metadata) -> str:
+    """Model with an invalid (zero) horizontal_resolution."""
+    file_name = tmpdir_factory.mktemp("data").join("zero_hres.nc")
+    return _write_model_file(
+        file_name,
+        file_metadata,
+        n_levels=2,
+        heights=[100, 5000],
+        resolution=0,
+        level_var=False,
+    )
+
+
+@pytest.fixture(scope="session")
+def model_file_no_clouds(tmpdir_factory, file_metadata) -> str:
+    """Model without cloud variables (ecmwf-open-like)."""
+    file_name = tmpdir_factory.mktemp("data").join("no_clouds.nc")
+    return _write_model_file(
+        file_name,
+        file_metadata,
+        n_levels=2,
+        heights=[100, 5000],
+        clouds=False,
+        level_var=False,
+    )
+
+
 @pytest.fixture(scope="session")
 def obs_file(tmpdir_factory, file_metadata) -> str:
     file_name = tmpdir_factory.mktemp("data").join("file.nc")
@@ -177,19 +288,19 @@ def regrid_file(tmpdir_factory, file_metadata) -> str:
     var[:] = 1
     var = root_grp.createVariable("horizontal_resolution", "f8")
     var[:] = 9
-    var = root_grp.createVariable("ecmwf_height", "f8", ("time", "level"))
+    var = root_grp.createVariable("model_height", "f8", ("time", "level"))
     var[:] = ma.array([[10, 14], [8, 14], [9, 15]])
-    var = root_grp.createVariable("ecmwf_forecast_time", "f8", "time")
+    var = root_grp.createVariable("model_forecast_time", "f8", "time")
     var[:] = ma.array([1, 5, 10])
-    var = root_grp.createVariable("ecmwf_cf", "f8", ("time", "level"))
+    var = root_grp.createVariable("model_cf", "f8", ("time", "level"))
     var[:] = ma.array([[0, 2], [3, 6], [5, 8]])
-    var = root_grp.createVariable("ecmwf_cf_cirrus", "f8", ("time", "level"))
+    var = root_grp.createVariable("model_cf_cirrus", "f8", ("time", "level"))
     var[:] = ma.array([[0, 2], [3, 6], [5, 7]])
-    var = root_grp.createVariable("ecmwf_cf_snow", "f8", ("time", "level"))
+    var = root_grp.createVariable("model_cf_snow", "f8", ("time", "level"))
     var[:] = ma.array([[0, 2], [4, 6], [5, 8]])
-    var = root_grp.createVariable("cf_ecmwf", "f8", ("time", "level"))
+    var = root_grp.createVariable("cf_V", "f8", ("time", "level"))
     var[:] = ma.array([[0, 2], [3, 6], [5, 8]])
-    var = root_grp.createVariable("cf_adv_ecmwf", "f8", ("time", "level"))
+    var = root_grp.createVariable("cf_V_adv", "f8", ("time", "level"))
     var[:] = ma.array([[0, 2], [3, 6], [5, 8]])
     var = root_grp.createVariable("temperature", "f8", ("time", "level"))
     var[:] = ma.array([[300, 301], [302, 299], [305, 298]])
diff --git a/cloudnetpy/model_evaluation/tests/unit/test_advance_methods.py b/cloudnetpy/model_evaluation/tests/unit/test_advance_methods.py
index 12fcd26f..da8f0679 100644
--- a/cloudnetpy/model_evaluation/tests/unit/test_advance_methods.py
+++ b/cloudnetpy/model_evaluation/tests/unit/test_advance_methods.py
@@ -8,15 +8,14 @@
 from cloudnetpy.model_evaluation.products.observation_products import ObservationManager
 
 MODEL = "ecmwf"
-OUTPUT_FILE = ""
 PRODUCT = "cf"
 
 
-@pytest.mark.parametrize("name", ("ecmwf_cf_cirrus",))
+@pytest.mark.parametrize("name", ("model_cf_cirrus",))
 def test_cf_cirrus_filter(obs_file, model_file, name) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    AdvanceProductMethods(model, obs)
     assert name in model.data
 
 
@@ -35,19 +34,19 @@ def test_cf_cirrus_filter(obs_file, model_file, name) -> None:
 )
 def test_getvar_from_object(obs_file, model_file, name, data) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
-    x = adv_pro.getvar_from_object(name)
-    testing.assert_array_almost_equal(x, data)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
+    result = advance.getvar_from_object(name)
+    testing.assert_array_almost_equal(result, data)
 
 
 @pytest.mark.parametrize("name", ("T",))
 def test_getvar_from_object_None(obs_file, model_file, name) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     with pytest.raises(KeyError):
-        adv_pro.getvar_from_object(name)
+        advance.getvar_from_object(name)
 
 
 @pytest.mark.parametrize(
@@ -59,83 +58,85 @@ def test_getvar_from_object_None(obs_file, model_file, name) -> None:
 )
 def test_set_frequency_parameters(obs_file, model_file, radar_f, values) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     obs.radar_freq = radar_f
-    x = adv_pro.set_frequency_parameters()
-    assert x == values
+    result = advance.set_frequency_parameters()
+    assert (result.ZT, result.T, result.Z, result.c) == values
 
 
 def test_fit_z_sensitivity(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
-    h = np.array([[5000, 9000, 13000], [10000, 15000, 20000], [8000, 12000, 16000]])
-    compare = ma.masked_invalid([[np.nan, 0.15, 0.5], [0.1, 1, np.nan], [0.15, 0, 1]])
-    x = adv_pro.fit_z_sensitivity(h)
-    testing.assert_array_almost_equal(x, compare)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
+    height = np.array(
+        [[5000, 9000, 13000], [10000, 15000, 20000], [8000, 12000, 16000]]
+    )
+    expected = ma.masked_invalid([[np.nan, 0.15, 0.5], [0.1, 1, np.nan], [0.15, 0, 1]])
+    result = advance.fit_z_sensitivity(height)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_filter_high_iwc_low_cf(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     cf = ma.array([0.0001, 0.0002, 0, 0.0001, 1, 0.0006])
     iwc = np.array([0.0, 0, 0, 0.2, 0.4, 0])
     lwc = np.array([0.0, 0.02, 0.01, 0, 0.01, 0.01])
-    compare = ma.array(
+    expected = ma.array(
         [0.0001, 0.0002, 0, -99, 1, 0.0006],
         mask=[False, False, False, True, False, False],
     )
-    x = adv_pro.filter_high_iwc_low_cf(cf, iwc, lwc)
-    testing.assert_array_almost_equal(x, compare)
+    result = advance.filter_high_iwc_low_cf(cf, iwc, lwc)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_filter_high_iwc_low_cf_no_ice(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     cf = ma.array([0.0001, 0.0002, 0, 0, 0, 0.0006])
     iwc = np.array([0.0, 0, 0, 0.2, 0.4, 0])
     lwc = np.array([0.0, 0.02, 0.01, 0, 0.01, 0.01])
     with pytest.raises(ValueError):
-        adv_pro.filter_high_iwc_low_cf(cf, iwc, lwc)
+        advance.filter_high_iwc_low_cf(cf, iwc, lwc)
 
 
 def test_mask_weird_indices(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     cf = ma.array([0.0001, 0.0002, 0, 0.0001, 1, 0.0006])
-    compare = ma.copy(cf)
+    expected = ma.copy(cf)
     iwc = np.array([0.0, 0, 0, 0.2, 0.4, 0])
     lwc = np.array([0.0, 0.02, 0.01, 0, 0.01, 0.01])
     ind = (iwc / cf > 0.5e-3) & (cf < 0.001)
     ind = ind | (iwc == 0) & (lwc == 0) & (cf == 0)
-    compare[ind] = ma.masked
-    x = adv_pro.mask_weird_indices(cf, iwc, lwc)
-    testing.assert_array_almost_equal(x, compare)
+    expected[ind] = ma.masked
+    result = advance.mask_weird_indices(cf, iwc, lwc)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_mask_weird_indices_values(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     cf = ma.array([0.0001, 0.0002, 0, 0.0001, 1, 0.0006])
     iwc = np.array([0.0, 0, 0, 0.2, 0.4, 0])
     lwc = np.array([0.0, 0.02, 0.01, 0, 0.01, 0.01])
-    compare = ma.array(
+    expected = ma.array(
         [0.0001, 0.0002, 0, -99, 1, 0.0006],
         mask=[False, False, False, True, False, False],
     )
-    x = adv_pro.mask_weird_indices(cf, iwc, lwc)
-    testing.assert_array_almost_equal(x, compare)
+    result = advance.mask_weird_indices(cf, iwc, lwc)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_find_ice_in_clouds(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     cf_f = np.array(
         [
             [-1, 0, 2, 3],
@@ -155,14 +156,14 @@ def test_find_ice_in_clouds(obs_file, model_file) -> None:
         ],
     )
     expected = np.array([31 / 3 * 1000, 200 / 2 * 1000])
-    x, _ = adv_pro.find_ice_in_clouds(cf_f, iwc, lwc)
-    testing.assert_array_almost_equal(x, expected)
+    result, _ = advance.find_ice_in_clouds(cf_f, iwc, lwc)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_get_ice_indices(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     cf_f = np.array(
         [
             [-1, 0, 2, 3],
@@ -181,92 +182,92 @@ def test_get_ice_indices(obs_file, model_file) -> None:
             [0, 1, 2, 3],
         ],
     )
-    x = np.array([0, 1])
-    y = np.array([3, 2])
-    expected = (x, y)
-    result = adv_pro.get_ice_indices(cf_f, iwc, lwc)
+    rows = np.array([0, 1])
+    cols = np.array([3, 2])
+    expected = (rows, cols)
+    result = advance.get_ice_indices(cf_f, iwc, lwc)
     testing.assert_array_almost_equal(result, expected)
 
 
 def test_iwc_variance(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
-    x = np.array([0, 1, 2])
-    y = np.array([0, 0, 1])
-    ind = (x, y)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
+    rows = np.array([0, 1, 2])
+    cols = np.array([0, 0, 1])
+    ind = (rows, cols)
     height = np.array([[1, 5], [2, 6], [3, 7]])
-    x = adv_pro.iwc_variance(height, ind)
-    assert len(x) == 3
+    result = advance.iwc_variance(height, ind)
+    assert len(result) == 3
 
 
 def test_calculate_variance_iwc(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     shear = np.array([[1, 1, 2, 1], [2, 2, 1, 0], [0, 0, 1, 0]])
     ind_arr = np.array([[0, 0, 1, 1], [0, 1, 0, 0], [1, 1, 0, 0]])
     ind = np.where(ind_arr > -1)
-    compare = 10 ** (0.3 * np.log10(model.resolution_h) - 0.04 * shear[ind] - 1.03)
-    x = adv_pro.calculate_variance_iwc(shear, tuple(ind))
-    testing.assert_array_almost_equal(x, compare)
+    expected = 10 ** (0.3 * np.log10(model.resolution_h) - 0.04 * shear[ind] - 1.03)
+    result = advance.calculate_variance_iwc(shear, tuple(ind))
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_calculate_wind_shear(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     u = np.array([[1, 2, 0, 1], [-1, 0, 1, -1], [1, 0, 1, -1]])
     v = np.array([[1, 0, 1, -1], [1, 2, -1, 0], [1, 2, 0, 1]])
     wind = np.sqrt(np.power(u, 2) + np.power(v, 2))
-    h = np.array([[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3]])
-    compare = np.array([[2, 2.83, 2.24, -2.24], [0, 0, 0, 0], [2, 0, -1, 1]])
-    x = adv_pro.calculate_wind_shear(wind, u, v, h)
-    testing.assert_array_almost_equal(np.round(x, 2), compare)
+    height = np.array([[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3]])
+    expected = np.array([[2, 2.83, 2.24, -2.24], [0, 1.41, 0.71, 1.41], [2, 0, -1, 1]])
+    result = advance.calculate_wind_shear(wind, u, v, height)
+    testing.assert_array_almost_equal(np.round(result, 2), expected)
 
 
 def test_calculate_iwc_distribution(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     n_std = 5
     n_dist = 250
     f_variance_iwc = 0.1
     cloud_iwc = 0.2
     finish = cloud_iwc + n_std * (np.sqrt(f_variance_iwc) * cloud_iwc)
-    compare = np.arange(0, finish, finish / (n_dist - 1))
-    x = adv_pro.calculate_iwc_distribution(cloud_iwc, f_variance_iwc)
-    testing.assert_array_almost_equal(x, compare)
+    expected = np.arange(0, finish, finish / (n_dist - 1))
+    result = advance.calculate_iwc_distribution(cloud_iwc, f_variance_iwc)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_gamma_distribution(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     iwc_dist = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6])
-    compare = np.zeros(iwc_dist.shape)
+    expected = np.zeros(iwc_dist.shape)
     f_variance_iwc = 0.1
     cloud_iwc = 0.2
     alpha = 1 / f_variance_iwc
     for i in range(len(iwc_dist)):
-        compare[i] = (
+        expected[i] = (
             1
             / gamma(alpha)
             * (alpha / cloud_iwc) ** alpha
             * iwc_dist[i] ** (alpha - 1)
             * ma.exp(-(alpha * iwc_dist[i] / cloud_iwc))
         )
-    x = adv_pro.gamma_distribution(iwc_dist, f_variance_iwc, cloud_iwc)
-    testing.assert_array_almost_equal(x, compare)
+    result = advance.gamma_distribution(iwc_dist, f_variance_iwc, cloud_iwc)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_filter_cirrus(obs_file, model_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    adv_pro = AdvanceProductMethods(model, str(model_file), obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    advance = AdvanceProductMethods(model, obs)
     cf_f = 0.7
-    p = np.array([1, 2, 3, 4, 5, 6])
+    p_iwc = np.array([1, 2, 3, 4, 5, 6])
     ind = np.array([0, 0, 1, 1, 0, 1], dtype=bool)
-    compare = (np.sum(p * ind) / np.sum(p)) * cf_f
-    x = adv_pro.filter_cirrus(p, ind, np.array([cf_f]))
-    testing.assert_almost_equal(x, compare)
+    expected = (np.sum(p_iwc * ind) / np.sum(p_iwc)) * cf_f
+    result = advance.filter_cirrus(p_iwc, ind, np.array([cf_f]))
+    testing.assert_almost_equal(result, expected)
diff --git a/cloudnetpy/model_evaluation/tests/unit/test_grid_methods.py b/cloudnetpy/model_evaluation/tests/unit/test_grid_methods.py
index 3ead2956..02d1a76c 100644
--- a/cloudnetpy/model_evaluation/tests/unit/test_grid_methods.py
+++ b/cloudnetpy/model_evaluation/tests/unit/test_grid_methods.py
@@ -7,7 +7,6 @@
 from cloudnetpy.model_evaluation.products.observation_products import ObservationManager
 
 MODEL = "ecmwf"
-OUTPUT_FILE = ""
 PRODUCT = "iwc"
 
 
@@ -17,157 +16,135 @@
         (
             "iwc",
             (
-                "ecmwf_iwc",
-                "iwc_ecmwf",
-                "iwc_att_ecmwf",
-                "iwc_rain_ecmwf",
-                "iwc_adv_ecmwf",
-                "iwc_att_adv_ecmwf",
-                "iwc_rain_adv_ecmwf",
+                "model_iwc",
+                "iwc",
+                "iwc_adv",
             ),
         ),
         (
             "cf",
             (
-                "ecmwf_cf",
-                "cf_A_ecmwf",
-                "cf_V_ecmwf",
-                "cf_A_adv_ecmwf",
-                "cf_V_adv_ecmwf",
+                "model_cf",
+                "cf_A",
+                "cf_V",
+                "cf_A_adv",
+                "cf_V_adv",
             ),
         ),
-        ("lwc", ("lwc_ecmwf", "lwc_ecmwf", "lwc_adv_ecmwf")),
+        ("lwc", ("model_lwc", "lwc", "lwc_adv")),
     ],
 )
 def test_generate_regrid_product(model_file, obs_file, product, variables) -> None:
     obs = ObservationManager(product, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, product)
+    model = ModelManager(str(model_file), MODEL, product)
     ProductGrid(model, obs)
     for var in variables:
         assert var in model.data
 
 
-@pytest.mark.parametrize("key, value", [("iwc", 3), ("lwc", 1), ("cf", 2)])
+@pytest.mark.parametrize("key, value", [("iwc", 1), ("lwc", 1), ("cf", 2)])
 def test_get_method_storage(key, value, model_file, obs_file) -> None:
     obs = ObservationManager(key, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, key)
-    obj = ProductGrid(model, obs)
-    x, y = obj._get_method_storage()
-    assert len(x.keys()) == value
+    model = ModelManager(str(model_file), MODEL, key)
+    grid = ProductGrid(model, obs)
+    storage, adv_storage = grid._get_method_storage()
+    assert len(storage.keys()) == value
 
 
-@pytest.mark.parametrize("key, value", [("iwc", 3), ("lwc", 1), ("cf", 2)])
+@pytest.mark.parametrize("key, value", [("iwc", 1), ("lwc", 1), ("cf", 2)])
 def test_get_method_storage_adv(key, value, model_file, obs_file) -> None:
     obs = ObservationManager(key, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, key)
-    obj = ProductGrid(model, obs)
-    x, y = obj._get_method_storage()
-    assert len(y.keys()) == value
+    model = ModelManager(str(model_file), MODEL, key)
+    grid = ProductGrid(model, obs)
+    storage, adv_storage = grid._get_method_storage()
+    assert len(adv_storage.keys()) == value
 
 
 @pytest.mark.parametrize("name", ["cf_V", "cf_A"])
 def test_cf_method_storage(name, model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    x, y = obj._cf_method_storage()
-    assert name in x
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
+    storage, adv_storage = grid._cf_method_storage()
+    assert name in storage
 
 
 @pytest.mark.parametrize("name", ["cf_V_adv", "cf_A_adv"])
 def test_cf_method_storage_adv(name, model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    x, y = obj._cf_method_storage()
-    assert name in y
-
-
-@pytest.mark.parametrize("name", ["iwc", "iwc_att", "iwc_rain"])
-def test_iwc_method_storage(name, model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    x, y = obj._iwc_method_storage()
-    assert name in x
-
-
-@pytest.mark.parametrize("name", ["iwc_adv", "iwc_att_adv", "iwc_rain_adv"])
-def test_iwc_method_storage_adv(name, model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    x, y = obj._iwc_method_storage()
-    assert name in y
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
+    storage, adv_storage = grid._cf_method_storage()
+    assert name in adv_storage
 
 
 def test_product_method_storage(model_file, obs_file) -> None:
     obs = ObservationManager("lwc", str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "lwc")
-    obj = ProductGrid(model, obs)
-    x, y = obj._product_method_storage()
-    assert "lwc" in x
+    model = ModelManager(str(model_file), MODEL, "lwc")
+    grid = ProductGrid(model, obs)
+    storage, adv_storage = grid._product_method_storage()
+    assert "lwc" in storage
 
 
 def test_product_method_storage_adv(model_file, obs_file) -> None:
     obs = ObservationManager("lwc", str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "lwc")
-    obj = ProductGrid(model, obs)
-    _, y = obj._product_method_storage()
-    assert "lwc_adv" in y
+    model = ModelManager(str(model_file), MODEL, "lwc")
+    grid = ProductGrid(model, obs)
+    _, adv_storage = grid._product_method_storage()
+    assert "lwc_adv" in adv_storage
 
 
 def test_regrid_cf_area(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
     data = ma.array([[1, 1, 1], [0, 1, 1], [0, 0, 1], [0, 0, 0]])
-    d = {"cf_A": ma.zeros((1, 1))}
-    d = obj._regrid_cf(d, 0, 0, data)
-    x = d["cf_A"]
-    assert x[0, 0] == 0.75
+    storage = {"cf_A": ma.zeros((1, 1))}
+    storage = grid._regrid_cf(storage, 0, 0, data)
+    result = storage["cf_A"]
+    assert result[0, 0] == 0.75
 
 
 def test_regrid_cf_area_masked(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
     data = ma.array([[1, 1, 1], [0, 1, 1], [0, 0, 1], [0, 0, 0]])
     data[1, :] = ma.masked
-    d = {"cf_A": ma.zeros((1, 1))}
-    d = obj._regrid_cf(d, 0, 0, data)
-    x = d["cf_A"]
-    assert round(x[0, 0], 3) == 0.667
+    storage = {"cf_A": ma.zeros((1, 1))}
+    storage = grid._regrid_cf(storage, 0, 0, data)
+    result = storage["cf_A"]
+    assert round(result[0, 0], 3) == 0.667
 
 
 def test_regrid_cf_area_all_masked(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
     data = ma.array([[1, 1, 1], [0, 1, 1], [0, 0, 1], [0, 0, 0]])
     data[:, :] = ma.masked
-    d = {"cf_A": ma.zeros((1, 1))}
-    d = obj._regrid_cf(d, 0, 0, data)
-    x = d["cf_A"]
-    testing.assert_equal(x, ma.masked)
+    storage = {"cf_A": ma.zeros((1, 1))}
+    storage = grid._regrid_cf(storage, 0, 0, data)
+    result = storage["cf_A"]
+    testing.assert_equal(result, ma.masked)
 
 
 def test_regrid_cf_area_nan(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
     data = ma.array([[1, 99, 1], [0, 1, 1], [99, 0, 1], [0, 0, 0]])
     data = ma.masked_where(data == 99, data)
-    d = {"cf_A": ma.zeros((1, 1))}
-    d = obj._regrid_cf(d, 0, 0, data)
-    x = d["cf_A"]
-    assert x[0, 0] == 0.75
+    storage = {"cf_A": ma.zeros((1, 1))}
+    storage = grid._regrid_cf(storage, 0, 0, data)
+    result = storage["cf_A"]
+    assert result[0, 0] == 0.75
 
 
 def test_regrid_cf_area_all_nan(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
     data = ma.array(
         [
             [1, 1, 1],
@@ -177,39 +154,39 @@ def test_regrid_cf_area_all_nan(model_file, obs_file) -> None:
         ],
         mask=True,
     )
-    d = {"cf_A": ma.zeros((1, 1))}
-    d = obj._regrid_cf(d, 0, 0, data)
-    x = d["cf_A"]
-    testing.assert_equal(x, ma.masked)
+    storage = {"cf_A": ma.zeros((1, 1))}
+    storage = grid._regrid_cf(storage, 0, 0, data)
+    result = storage["cf_A"]
+    testing.assert_equal(result, ma.masked)
 
 
 def test_regrid_cf_volume(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
     data = ma.array([[1, 1, 1], [0, 1, 1], [0, 0, 1], [0, 0, 0]])
-    d = {"cf_V": ma.zeros((1, 1))}
-    d = obj._regrid_cf(d, 0, 0, data)
-    x = d["cf_V"]
-    assert x[0, 0] == 0.5
+    storage = {"cf_V": ma.zeros((1, 1))}
+    storage = grid._regrid_cf(storage, 0, 0, data)
+    result = storage["cf_V"]
+    assert result[0, 0] == 0.5
 
 
 def test_regrid_cf_volume_nan(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
     data = ma.array([[1, 99, 1], [0, 1, 1], [99, 0, 1], [0, 0, 0]])
     data = ma.masked_where(data == 99, data)
-    d = {"cf_V": ma.zeros((1, 1))}
-    d = obj._regrid_cf(d, 0, 0, data)
-    x = d["cf_V"]
-    assert x[0, 0] == 0.5
+    storage = {"cf_V": ma.zeros((1, 1))}
+    storage = grid._regrid_cf(storage, 0, 0, data)
+    result = storage["cf_V"]
+    assert result[0, 0] == 0.5
 
 
 def test_regrid_cf_volume_all_nan(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
     data = ma.array(
         [
             [99, 99, 99],
@@ -219,42 +196,42 @@ def test_regrid_cf_volume_all_nan(model_file, obs_file) -> None:
         ],
     )
     data = ma.masked_where(data == 99, data)
-    d = {"cf_V": ma.zeros((1, 1))}
-    d = obj._regrid_cf(d, 0, 0, data)
-    x = d["cf_V"]
-    testing.assert_equal(x, ma.masked)
+    storage = {"cf_V": ma.zeros((1, 1))}
+    storage = grid._regrid_cf(storage, 0, 0, data)
+    result = storage["cf_V"]
+    testing.assert_equal(result, ma.masked)
 
 
 def test_regrid_cf_volume_masked(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
     data = ma.array([[1, 1, 1], [0, 1, 1], [0, 0, 1], [0, 0, 0]])
     data[1, :] = ma.masked
-    d = {"cf_V": ma.zeros((1, 1))}
-    d = obj._regrid_cf(d, 0, 0, data)
-    x = d["cf_V"]
-    assert round(x[0, 0], 3) == 0.444
+    storage = {"cf_V": ma.zeros((1, 1))}
+    storage = grid._regrid_cf(storage, 0, 0, data)
+    result = storage["cf_V"]
+    assert round(result[0, 0], 3) == 0.444
 
 
 def test_regrid_cf_volume_all_masked(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
     data = ma.array([[1, 1, 1], [0, 1, 1], [0, 0, 1], [0, 0, 0]])
     data[:, :] = ma.masked
-    d = {"cf_V": ma.zeros((1, 1))}
-    d = obj._regrid_cf(d, 0, 0, data)
-    x = d["cf_V"]
-    testing.assert_equal(x, ma.masked)
+    storage = {"cf_V": ma.zeros((1, 1))}
+    storage = grid._regrid_cf(storage, 0, 0, data)
+    result = storage["cf_V"]
+    testing.assert_equal(result, ma.masked)
 
 
 def test_reshape_data_to_window(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    xnd = np.array([1, 1, 1, 0, 0, 0])
-    ynd = np.array([1, 1, 0, 0])
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
+    x_ind = np.array([1, 1, 1, 0, 0, 0])
+    y_ind = np.array([1, 1, 0, 0])
     ind = np.array(
         [
             [1, 1, 0, 0],
@@ -266,7 +243,7 @@ def test_reshape_data_to_window(model_file, obs_file) -> None:
         ],
         dtype=bool,
     )
-    obj._obs_data = np.array(
+    grid._obs_data = np.array(
         [
             [1, 2, 3, 4],
             [11, 22, 33, 44],
@@ -276,18 +253,18 @@ def test_reshape_data_to_window(model_file, obs_file) -> None:
             [555, 666, 777, 888],
         ],
     )
-    x = obj._reshape_data_to_window(ind, xnd, ynd)
-    compare = np.array([[1, 2], [11, 22], [111, 222]])
-    assert x is not None
-    testing.assert_array_almost_equal(x, compare)
+    result = grid._reshape_data_to_window(ind, x_ind, y_ind)
+    expected = np.array([[1, 2], [11, 22], [111, 222]])
+    assert result is not None
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_reshape_data_to_window_middle(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    xnd = np.array([0, 0, 1, 1, 1, 0])
-    ynd = np.array([0, 1, 1, 0])
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
+    x_ind = np.array([0, 0, 1, 1, 1, 0])
+    y_ind = np.array([0, 1, 1, 0])
     ind = np.array(
         [
             [0, 0, 0, 0],
@@ -299,7 +276,7 @@ def test_reshape_data_to_window_middle(model_file, obs_file) -> None:
         ],
         dtype=bool,
     )
-    obj._obs_data = np.array(
+    grid._obs_data = np.array(
         [
             [1, 2, 3, 4],
             [11, 22, 33, 44],
@@ -309,17 +286,17 @@ def test_reshape_data_to_window_middle(model_file, obs_file) -> None:
             [555, 666, 777, 888],
         ],
     )
-    x = obj._reshape_data_to_window(ind, xnd, ynd)
-    compare = np.array([[222, 333], [6, 7], [66, 77]])
-    assert x is not None
-    testing.assert_array_almost_equal(x, compare)
+    result = grid._reshape_data_to_window(ind, x_ind, y_ind)
+    expected = np.array([[222, 333], [6, 7], [66, 77]])
+    assert result is not None
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_reshape_data_to_window_empty(model_file, obs_file) -> None:
     obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    xnd = np.array(
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    grid = ProductGrid(model, obs)
+    x_ind = np.array(
         [
             1,
             1,
@@ -329,357 +306,31 @@ def test_reshape_data_to_window_empty(model_file, obs_file) -> None:
             0,
         ],
     )
-    ynd = np.array([0, 0, 0, 0])
+    y_ind = np.array([0, 0, 0, 0])
     ind = np.array([1, 1, 0, 0], dtype=bool)
-    x = obj._reshape_data_to_window(ind, xnd, ynd)
-    assert x is None
-
-
-def test_regrid_iwc(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array([[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3], [4, 4, 4, 3]])
-    d = {"iwc": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array([[0, 1, 1, 1]], dtype=bool)
-    no_rain: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 0], [0, 0, 0, 0]],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc"]
-    testing.assert_almost_equal(x[0, 0], 1.4)
-
-
-def test_regrid_iwc_nan(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array(
-        [[1, 1, 99, 1], [2, 99, 2, 2], [3, 3, 3, 3], [4, 4, 4, 99]],
-    )
-    obj._obs_data = ma.masked_where(obj._obs_data == 99, obj._obs_data)
-    d = {"iwc": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array([[0, 1, 1, 1]], dtype=bool)
-    no_rain: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 0], [0, 0, 0, 0]],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc"]
-    testing.assert_almost_equal(x[0, 0], 1.5)
-
-
-def test_regrid_iwc_all_nan(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array(
-        [
-            [99, 99, 99, 99],
-            [99, 99, 99, 99],
-            [99, 99, 99, 99],
-            [99, 99, 99, 99],
-        ],
-    )
-    obj._obs_data = ma.masked_where(obj._obs_data == 99, obj._obs_data)
-    d = {"iwc": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array([[0, 1, 1, 1]], dtype=bool)
-    no_rain: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 0], [0, 0, 0, 0]],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    assert d["iwc"][0, 0].mask == True
-
-
-def test_regrid_iwc_masked(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array([[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
-    obj._obs_data[1, :] = ma.masked
-    d = {"iwc": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array([[0, 1, 1, 1]], dtype=bool)
-    no_rain: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 0], [0, 0, 0, 0]],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc"]
-    testing.assert_almost_equal(x[0, 0], 1.0)
-
-
-def test_regrid_iwc_all_masked(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array(
-        [[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]], mask=True
-    )
-    d = {"iwc": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array([[0, 1, 1, 1]], dtype=bool)
-    no_rain: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 0], [0, 0, 0, 0]],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    assert d["iwc"][0, 0].mask == True
-
-
-def test_regrid_iwc_none(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array([[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
-    d = {"iwc": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array([[0, 1, 1, 1]], dtype=bool)
-    no_rain: ma.MaskedArray = ma.array(
-        [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    assert d["iwc"][0, 0].mask == True
-
-
-def test_regrid_iwc_att(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    d = {"iwc_att": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array([[0, 1, 1, 1]], dtype=bool)
-    no_rain: ma.MaskedArray = ma.array(
-        [
-            [0, 1, 1, 1],
-            [0, 0, 1, 1],
-            [0, 1, 1, 1],
-            [0, 0, 1, 1],
-            [0, 0, 0, 0],
-            [0, 0, 0, 0],
-        ],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc_att"]
-    testing.assert_almost_equal(x[0, 0], 0.018)
-
-
-def test_regrid_iwc_att_masked(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_obj.data["iwc_att"][:].mask = ma.array(
-        [
-            [1, 0, 1, 0],
-            [0, 1, 1, 0],
-            [1, 0, 0, 1],
-            [0, 1, 1, 1],
-            [1, 1, 0, 0],
-            [0, 1, 0, 1],
-        ],
-        dtype=bool,
-    )
-    d = {"iwc_att": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array([[0, 1, 1, 1]], dtype=bool)
-    no_rain: ma.MaskedArray = ma.array(
-        [
-            [0, 1, 1, 1],
-            [0, 0, 1, 1],
-            [0, 1, 1, 1],
-            [0, 0, 1, 1],
-            [0, 0, 0, 0],
-            [0, 0, 0, 0],
-        ],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc_att"]
-    testing.assert_almost_equal(x[0, 0], 0.018)
-
-
-def test_regrid_iwc_att_all_masked(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_obj.data["iwc_att"][:].mask = ma.array(
-        [
-            [1, 1, 1, 1],
-            [1, 1, 1, 1],
-            [1, 1, 1, 1],
-            [1, 1, 1, 1],
-            [1, 1, 1, 1],
-            [1, 1, 1, 1],
-        ],
-        dtype=bool,
-    )
-    d = {"iwc_att": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array([[0, 1, 1, 1]], dtype=bool)
-    no_rain: ma.MaskedArray = ma.array(
-        [
-            [0, 1, 1, 1],
-            [0, 0, 1, 1],
-            [0, 1, 1, 1],
-            [0, 0, 1, 1],
-            [0, 0, 0, 0],
-            [0, 0, 0, 0],
-        ],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc_att"]
-    # Simo: not sure if this should be masked or not
-    assert x[0, 0] == 0.01
-
-
-def test_regrid_iwc_att_none(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    d = {"iwc_att": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array([[0, 1, 1, 1]], dtype=bool)
-    no_rain: ma.MaskedArray = ma.array(
-        [
-            [0, 0, 0, 0],
-            [0, 0, 0, 0],
-            [0, 0, 0, 0],
-            [0, 0, 0, 0],
-            [0, 0, 0, 0],
-            [0, 0, 0, 0],
-        ],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc_att"]
-    assert x[0, 0].mask == True
-
-
-def test_regrid_iwc_rain(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array([[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3], [4, 4, 4, 3]])
-    d = {"iwc_rain": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 1, 1, 1], [0, 0, 1, 1]], dtype=bool
-    )
-    no_rain: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 1, 1, 1], [0, 0, 1, 1]],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc_rain"]
-    testing.assert_almost_equal(x[0, 0], 2.3)
-
-
-def test_regrid_iwc_rain_nan(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array(
-        [
-            [1, 99, 1, 1],
-            [2, 2, 2, 99],
-            [3, 3, 3, 3],
-            [99, 4, 4, 99],
-        ],
-    )
-    obj._obs_data = ma.masked_where(obj._obs_data == 99, obj._obs_data)
-    d = {"iwc_rain": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 1, 1, 1], [0, 0, 1, 1]], dtype=bool
-    )
-    no_rain: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 1, 1, 1], [0, 0, 1, 1]],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc_rain"]
-    testing.assert_almost_equal(round(x[0, 0], 3), 2.429)
-
-
-def test_regrid_iwc_rain_all_nan(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array(
-        [
-            [99, 99, 99, 99],
-            [99, 99, 99, 99],
-            [99, 99, 99, 99],
-            [99, 99, 99, 99],
-        ],
-    )
-    obj._obs_data = ma.masked_where(obj._obs_data == 99, obj._obs_data)
-    d = {"iwc_rain": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 1, 1, 1], [0, 0, 1, 1]], dtype=bool
-    )
-    no_rain: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 1, 1, 1], [0, 0, 1, 1]],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc_rain"]
-    assert x[0, 0].mask == True
-
-
-def test_regrid_iwc_rain_masked(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array([[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
-    obj._obs_data[2, :] = ma.masked
-    d = {"iwc_rain": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 1, 1, 1], [0, 0, 1, 1]], dtype=bool
-    )
-    no_rain: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 1, 1, 1], [0, 0, 1, 1]],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc_rain"]
-    testing.assert_almost_equal(round(x[0, 0], 3), 2.143)
-
-
-def test_regrid_iwc_rain_all_masked(model_file, obs_file) -> None:
-    obs = ObservationManager(PRODUCT, str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array([[1, 3, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
-    obj._obs_data[:, :] = ma.masked
-    d = {"iwc_rain": ma.zeros((1, 1))}
-    ind: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 1, 1, 1], [0, 0, 1, 1]], dtype=bool
-    )
-    no_rain: ma.MaskedArray = ma.array(
-        [[0, 1, 1, 1], [0, 0, 1, 1], [0, 1, 1, 1], [0, 0, 1, 1]],
-        dtype=bool,
-    )
-    d = obj._regrid_iwc(d, 0, 0, ind, no_rain)
-    x = d["iwc_rain"]
-    assert x[0, 0].mask == True
+    result = grid._reshape_data_to_window(ind, x_ind, y_ind)
+    assert result is None
 
 
 def test_regrid_product(model_file, obs_file) -> None:
     obs = ObservationManager("lwc", str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "lwc")
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array([[1, 1, 1, 1], [2, 1, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
-    d = {"lwc": ma.zeros((1, 1))}
+    model = ModelManager(str(model_file), MODEL, "lwc")
+    grid = ProductGrid(model, obs)
+    grid._obs_data = ma.array([[1, 1, 1, 1], [2, 1, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
+    storage = {"lwc": ma.zeros((1, 1))}
     ind: ma.MaskedArray = ma.array(
         [[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 0], [0, 0, 0, 0]], dtype=bool
     )
-    d = obj._regrid_product(d, 0, 0, ind)
-    x = d["lwc"]
-    testing.assert_almost_equal(x[0, 0], 1.4)
+    storage = grid._regrid_product(storage, 0, 0, ind)
+    result = storage["lwc"]
+    testing.assert_almost_equal(result[0, 0], 1.4)
 
 
 def test_regrid_product_nan(model_file, obs_file) -> None:
     obs = ObservationManager("lwc", str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "lwc")
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array(
+    model = ModelManager(str(model_file), MODEL, "lwc")
+    grid = ProductGrid(model, obs)
+    grid._obs_data = ma.array(
         [
             [1, 99, 1, 1],
             [99, 1, 2, 2],
@@ -687,20 +338,20 @@ def test_regrid_product_nan(model_file, obs_file) -> None:
             [4, 99, 4, 4],
         ],
     )
-    obj._obs_data = ma.masked_where(obj._obs_data == 99, obj._obs_data)
-    d = {"lwc": ma.zeros((1, 1))}
+    grid._obs_data = ma.masked_where(grid._obs_data == 99, grid._obs_data)
+    storage = {"lwc": ma.zeros((1, 1))}
     ind: ma.MaskedArray = ma.array(
         [[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 0], [0, 0, 0, 0]], dtype=bool
     )
-    d = obj._regrid_product(d, 0, 0, ind)
-    x = d["lwc"]
-    testing.assert_almost_equal(x[0, 0], 1.5)
+    storage = grid._regrid_product(storage, 0, 0, ind)
+    result = storage["lwc"]
+    testing.assert_almost_equal(result[0, 0], 1.5)
 
 
 def test_regrid_product_all_nan(model_file, obs_file) -> None:
     obs = ObservationManager("lwc", str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "lwc")
-    obj = ProductGrid(model, obs)
+    model = ModelManager(str(model_file), MODEL, "lwc")
+    grid = ProductGrid(model, obs)
     data = ma.array(
         [
             [99, 99, 99, 99],
@@ -709,73 +360,70 @@ def test_regrid_product_all_nan(model_file, obs_file) -> None:
             [99, 99, 99, 99],
         ],
     )
-    obj._obs_data = ma.masked_where(data == 99, data)
-    d = {"lwc": ma.zeros((1, 1))}
+    grid._obs_data = ma.masked_where(data == 99, data)
+    storage = {"lwc": ma.zeros((1, 1))}
     ind = np.array([[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 0], [0, 0, 0, 0]], dtype=bool)
-    d = obj._regrid_product(d, 0, 0, ind)
-    assert d["lwc"][0, 0].mask == True
+    storage = grid._regrid_product(storage, 0, 0, ind)
+    assert storage["lwc"][0, 0].mask == True
 
 
 def test_regrid_product_masked(model_file, obs_file) -> None:
     obs = ObservationManager("lwc", str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "lwc")
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array([[1, 1, 1, 1], [2, 1, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
-    obj._obs_data[2, :] = ma.masked
-    d = {"lwc": ma.zeros((1, 1))}
+    model = ModelManager(str(model_file), MODEL, "lwc")
+    grid = ProductGrid(model, obs)
+    grid._obs_data = ma.array([[1, 1, 1, 1], [2, 1, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
+    grid._obs_data[2, :] = ma.masked
+    storage = {"lwc": ma.zeros((1, 1))}
     ind = np.array([[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 0], [0, 0, 0, 0]], dtype=bool)
-    d = obj._regrid_product(d, 0, 0, ind)
-    x = d["lwc"]
-    testing.assert_almost_equal(x[0, 0], 1.4)
+    storage = grid._regrid_product(storage, 0, 0, ind)
+    result = storage["lwc"]
+    testing.assert_almost_equal(result[0, 0], 1.4)
 
 
 def test_regrid_product_all_masked(model_file, obs_file) -> None:
     obs = ObservationManager("lwc", str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "lwc")
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array([[1, 1, 1, 1], [2, 1, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
-    obj._obs_data[:, :] = ma.masked
-    d = {"lwc": ma.zeros((1, 1))}
+    model = ModelManager(str(model_file), MODEL, "lwc")
+    grid = ProductGrid(model, obs)
+    grid._obs_data = ma.array([[1, 1, 1, 1], [2, 1, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
+    grid._obs_data[:, :] = ma.masked
+    storage = {"lwc": ma.zeros((1, 1))}
     ind = np.array([[0, 1, 1, 1], [0, 0, 1, 1], [0, 0, 0, 0], [0, 0, 0, 0]], dtype=bool)
-    d = obj._regrid_product(d, 0, 0, ind)
-    x = d["lwc"]
-    testing.assert_almost_equal(x, ma.masked)
+    storage = grid._regrid_product(storage, 0, 0, ind)
+    result = storage["lwc"]
+    testing.assert_almost_equal(result, ma.masked)
 
 
 def test_regrid_product_none(model_file, obs_file) -> None:
     obs = ObservationManager("lwc", str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "lwc")
-    obj = ProductGrid(model, obs)
-    obj._obs_data = ma.array([[1, 1, 1, 1], [2, 1, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
-    d = {"lwc": ma.zeros((1, 1))}
+    model = ModelManager(str(model_file), MODEL, "lwc")
+    grid = ProductGrid(model, obs)
+    grid._obs_data = ma.array([[1, 1, 1, 1], [2, 1, 2, 2], [3, 3, 3, 3], [4, 4, 4, 4]])
+    storage = {"lwc": ma.zeros((1, 1))}
     ind = np.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], dtype=bool)
-    d = obj._regrid_product(d, 0, 0, ind)
-    x = d["lwc"]
-    assert x[0, 0].mask == True
+    storage = grid._regrid_product(storage, 0, 0, ind)
+    result = storage["lwc"]
+    assert result[0, 0].mask == True
 
 
 @pytest.mark.parametrize("product", ["cf_A", "cf_V", "cf_A_adv", "cf_V_adv"])
 def test_append_data2object_cf(product, model_file, obs_file) -> None:
     obs = ObservationManager("cf", str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "cf")
+    model = ModelManager(str(model_file), MODEL, "cf")
     ProductGrid(model, obs)
-    assert product + "_" + MODEL in model.data
+    assert product in model.data
 
 
-@pytest.mark.parametrize(
-    "product",
-    ["iwc", "iwc_att", "iwc_rain", "iwc_adv", "iwc_att_adv", "iwc_rain_adv"],
-)
+@pytest.mark.parametrize("product", ["iwc", "iwc_adv"])
 def test_append_data2object_iwc(product, model_file, obs_file) -> None:
     obs = ObservationManager("iwc", str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "iwc")
+    model = ModelManager(str(model_file), MODEL, "iwc")
     ProductGrid(model, obs)
-    assert product + "_" + MODEL in model.data
+    assert product in model.data
 
 
 @pytest.mark.parametrize("product", ["lwc", "lwc_adv"])
 def test_append_data2object_lwc(product, model_file, obs_file) -> None:
     obs = ObservationManager("lwc", str(obs_file))
-    model = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "lwc")
+    model = ModelManager(str(model_file), MODEL, "lwc")
     ProductGrid(model, obs)
-    assert product + "_" + MODEL in model.data
+    assert product in model.data
diff --git a/cloudnetpy/model_evaluation/tests/unit/test_model_products.py b/cloudnetpy/model_evaluation/tests/unit/test_model_products.py
index de3a242e..c62353ad 100644
--- a/cloudnetpy/model_evaluation/tests/unit/test_model_products.py
+++ b/cloudnetpy/model_evaluation/tests/unit/test_model_products.py
@@ -3,106 +3,106 @@
 import pytest
 from numpy import testing
 
+from cloudnetpy.constants import RS
+from cloudnetpy.exceptions import ModelDataError
 from cloudnetpy.model_evaluation.products.model_products import ModelManager
 
 MODEL = "ecmwf"
-OUTPUT_FILE = ""
 PRODUCT = "iwc"
 
 
-@pytest.mark.parametrize(
-    "cycle, model, answer",
-    [("test_file_12-23", "icon", "_12-23"), ("test_file", "ecmwf", "")],
-)
-def test_read_cycle_name(cycle, model, answer, model_file) -> None:
-    obj = ModelManager(str(model_file), model, OUTPUT_FILE, PRODUCT)
-    x = obj._read_cycle_name(cycle)
-    assert x == answer
-
-
 def test_get_cf(model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj._get_cf()
-    assert f"{MODEL}_cf" in obj.data
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    model._get_cf()
+    assert "model_cf" in model.data
 
 
 def test_get_iwc(model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj._get_iwc()
-    assert f"{MODEL}_iwc" in obj.data
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    model._get_iwc()
+    assert "model_iwc" in model.data
 
 
 def test_get_lwc(model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, "lwc")
-    obj._get_lwc()
-    assert f"{MODEL}_lwc" in obj.data
+    model = ModelManager(str(model_file), MODEL, "lwc")
+    model._get_lwc()
+    assert "model_lwc" in model.data
 
 
 def test_read_config(model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    var = obj.get_model_var_names(("p",))
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    var = model.get_model_var_names(("p",))
     assert "pressure" in var
-    var = obj.get_model_var_names(("T",))
+    var = model.get_model_var_names(("T",))
     assert "temperature" in var
 
 
 @pytest.mark.parametrize("key", ["pressure", "temperature"])
 def test_set_variables(key, model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    var = obj.getvar(key)
-    x = netCDF4.Dataset(model_file).variables[key]
-    testing.assert_almost_equal(x, var)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    actual = model.getvar(key)
+    expected = netCDF4.Dataset(model_file).variables[key]
+    testing.assert_almost_equal(expected, actual)
 
 
 @pytest.mark.parametrize("p, T, q", [(1, 2, 3), (20, 40, 80), (0.3, 0.6, 0.9)])
 def test_calc_water_content(p, T, q, model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    x = q * p / (287 * T)
-    testing.assert_almost_equal(x, obj._calc_water_content(q, p, T))
-
-
-@pytest.mark.parametrize(
-    "key",
-    ["time", "level", "horizontal_resolution", "latitude", "longitude"],
-)
-def test_add_common_variables_false(key, model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj._is_file = False
-    obj._add_variables()
-    assert key in obj.data
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    result = q * p / (RS * T)
+    testing.assert_almost_equal(result, model._calc_water_content(q, p, T))
 
 
 @pytest.mark.parametrize(
     "key",
     ["time", "level", "horizontal_resolution", "latitude", "longitude"],
 )
-def test_add_common_variables_true(key, model_file, regrid_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, regrid_file, PRODUCT)
-    obj._is_file = True
-    obj._add_variables()
-    assert key not in obj.data
+def test_add_common_variables(key, model_file) -> None:
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    assert key in model.data
 
 
 @pytest.mark.parametrize("key", ["height", "forecast_time"])
-def test_add_cycle_variables_no_products(key, model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    obj._is_file = False
-    obj._add_variables()
-    assert f"{MODEL}_{key}" in obj.data
+def test_add_model_variables(key, model_file) -> None:
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    assert f"model_{key}" in model.data
 
 
 def test_cut_off_extra_levels(model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
+    # The fixture has 2 levels, both below the altitude limit, so nothing is cut.
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
     data = np.array([np.arange(100), np.arange(100)])
-    compare = np.array([np.arange(88), np.arange(88)])
-    x = obj.cut_off_extra_levels(data)
-    testing.assert_array_almost_equal(x, compare)
+    result = model.cut_off_extra_levels(data)
+    testing.assert_array_almost_equal(result, data[:, :2])
+
+
+def test_cut_off_extra_levels_drops_high_levels(model_file_tall) -> None:
+    # The tall fixture has 5 levels, of which 3 are below the altitude limit.
+    model = ModelManager(str(model_file_tall), MODEL, PRODUCT)
+    assert model._n_levels == 3
+    data = np.arange(5 * 3).reshape(3, 5)
+    result = model.cut_off_extra_levels(data)
+    testing.assert_array_almost_equal(result, data[:, :3])
+
+
+def test_missing_horizontal_resolution_raises(model_file_no_hres) -> None:
+    with pytest.raises(ModelDataError, match="horizontal_resolution"):
+        ModelManager(str(model_file_no_hres), MODEL, "cf")
+
+
+def test_zero_horizontal_resolution_raises(model_file_zero_hres) -> None:
+    with pytest.raises(ModelDataError, match="invalid horizontal_resolution"):
+        ModelManager(str(model_file_zero_hres), MODEL, "cf")
+
+
+def test_missing_product_variable_raises(model_file_no_clouds) -> None:
+    with pytest.raises(ModelDataError, match="cloud_fraction"):
+        ModelManager(str(model_file_no_clouds), MODEL, "cf")
 
 
 def test_calculate_wind_speed(model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    u = obj.getvar("uwind")
-    v = obj.getvar("vwind")
-    compare = np.sqrt(u.data**2 + v.data**2)  # type: ignore
-    x = obj._calculate_wind_speed()
-    testing.assert_array_almost_equal(x, compare)
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    u = model.getvar("uwind")
+    v = model.getvar("vwind")
+    expected = np.sqrt(u.data**2 + v.data**2)  # type: ignore
+    result = model._calculate_wind_speed()
+    testing.assert_array_almost_equal(result, expected)
diff --git a/cloudnetpy/model_evaluation/tests/unit/test_observation_products.py b/cloudnetpy/model_evaluation/tests/unit/test_observation_products.py
index 383e7855..c0f34777 100644
--- a/cloudnetpy/model_evaluation/tests/unit/test_observation_products.py
+++ b/cloudnetpy/model_evaluation/tests/unit/test_observation_products.py
@@ -23,46 +23,45 @@ def __init__(self) -> None:
 
 
 def test_get_date(obs_file) -> None:
-    obj = ObservationManager(PRODUCT, str(obs_file))
-    date = datetime(2019, 5, 23, 0, 0, 0, tzinfo=timezone.utc)
-    assert obj._get_date() == date
+    obs = ObservationManager(PRODUCT, str(obs_file))
+    expected_date = datetime(2019, 5, 23, 0, 0, 0, tzinfo=timezone.utc)
+    assert obs._get_date() == expected_date
 
 
 @pytest.mark.parametrize("key", ["iwc", "lwc", "cf"])
 def test_generate_product(key, obs_file) -> None:
-    obj = ObservationManager(key, str(obs_file))
-    obj._generate_product()
-    assert key in obj.data
+    obs = ObservationManager(key, str(obs_file))
+    obs._generate_product()
+    assert key in obs.data
 
 
 def test_add_height(obs_file) -> None:
-    obj = ObservationManager(PRODUCT, str(obs_file))
-    obj._generate_product()
-    assert "height" in obj.data
+    obs = ObservationManager(PRODUCT, str(obs_file))
+    obs._generate_product()
+    assert "height" in obs.data
 
 
 def test_generate_cf(obs_file) -> None:
-    obj = ObservationManager("cf", str(obs_file))
-    x = obj._generate_cf()
-    compare = ma.array(
+    obs = ObservationManager("cf", str(obs_file))
+    result = obs._generate_cf()
+    expected = ma.array(
         [
             [0, 1, 0, 0],
-            [0, 0, 0, 0],
-            [0, 0, 0, 1],
             [0, 0, 0, 1],
             [1, 0, 0, 0],
+            [0, 0, 0, 1],
+            [0, 0, 0, 0],
             [0, 0, 0, 0],
         ],
     )
-    compare[~obj._rain_index(), :] = ma.masked
-    testing.assert_array_almost_equal(compare, x)
+    testing.assert_array_almost_equal(expected, result)
 
 
 def test_basic_cloud_mask(obs_file) -> None:
     cat = CategorizeBits(str(obs_file))
-    obj = ObservationManager("cf", str(obs_file))
-    x = obj._classify_basic_mask(cat.category_bits)
-    compare = np.array(
+    obs = ObservationManager("cf", str(obs_file))
+    result = obs._classify_basic_mask(cat.category_bits)
+    expected = np.array(
         [
             [0, 1, 2, 0],
             [2, 0, 0, 1],
@@ -72,15 +71,15 @@ def test_basic_cloud_mask(obs_file) -> None:
             [7, 2, 0, 7],
         ],
     )
-    testing.assert_array_almost_equal(x, compare)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_mask_cloud_bits(obs_file) -> None:
     cat = CategorizeBits(str(obs_file))
-    obj = ObservationManager("cf", str(obs_file))
-    mask = obj._classify_basic_mask(cat.category_bits)
-    compare = obj._mask_cloud_bits(mask)
-    x = np.array(
+    obs = ObservationManager("cf", str(obs_file))
+    mask = obs._classify_basic_mask(cat.category_bits)
+    expected = obs._mask_cloud_bits(mask)
+    result = np.array(
         [
             [0, 1, 0, 0],
             [0, 0, 0, 1],
@@ -90,69 +89,30 @@ def test_mask_cloud_bits(obs_file) -> None:
             [0, 0, 0, 0],
         ],
     )
-    testing.assert_array_almost_equal(x, compare)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_basic_cloud_mask_all_values(obs_file) -> None:
     cat = CatBits()
-    obj = ObservationManager("cf", str(obs_file))
-    x = obj._classify_basic_mask(cat.category_bits)  # type: ignore
-    compare = np.array([[8, 7, 6, 1, 3, 1], [0, 1, 7, 5, 2, 4]])
-    testing.assert_array_almost_equal(x, compare)
+    obs = ObservationManager("cf", str(obs_file))
+    result = obs._classify_basic_mask(cat.category_bits)  # type: ignore
+    expected = np.array([[8, 7, 6, 1, 3, 1], [0, 1, 7, 5, 2, 4]])
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_mask_cloud_bits_all_values(obs_file) -> None:
     cat = CatBits()
-    obj = ObservationManager("cf", str(obs_file))
-    mask = obj._classify_basic_mask(cat.category_bits)  # type: ignore
-    x = obj._mask_cloud_bits(mask)
-    compare = np.array([[0, 0, 0, 1, 1, 1], [0, 1, 0, 1, 0, 1]])
-    testing.assert_array_almost_equal(x, compare)
-
-
-def test_check_rainrate(obs_file) -> None:
-    obj = ObservationManager("cf", str(obs_file))
-    x = obj._check_rainrate()
-    assert x is True
-
-
-def test_get_rainrate_threshold(obs_file) -> None:
-    obj = ObservationManager("cf", str(obs_file))
-    x = obj._get_rainrate_threshold()
-    assert x == 8
-
-
-def test_rain_index(obs_file) -> None:
-    obj = ObservationManager("cf", str(obs_file))
-    x = obj._rain_index()
-    compare = np.array([0, 0, 0, 1, 0, 1], dtype=bool)
-    testing.assert_array_almost_equal(x, compare)
-
-
-@pytest.mark.parametrize("key", ["iwc", "iwc_att", "iwc_rain"])
-def test_generate_iwc_masks(key, obs_file) -> None:
-    obj = ObservationManager(PRODUCT, str(obs_file))
-    obj._generate_iwc_masks()
-    assert key in obj.data
-
-
-def test_get_rain_iwc(obs_file) -> None:
-    obj = ObservationManager("iwc", str(obs_file))
-    iwc_status = obj.getvar("iwc_retrieval_status")
-    x = np.zeros(iwc_status.shape)
-    x[iwc_status == 5] = 1
-    x = np.any(x, axis=1)
-    obj._get_rain_iwc(iwc_status[:])
-    compare = obj.data["iwc_rain"]
-    testing.assert_array_almost_equal(x, compare[:])
-
-
-def test_mask_iwc_att(obs_file) -> None:
-    obj = ObservationManager("iwc", str(obs_file))
-    iwc = obj.getvar("iwc")
-    iwc_status = obj.getvar("iwc_retrieval_status")
-    x = ma.copy(iwc)
-    obj._mask_iwc_att(iwc, iwc_status)
-    compare = obj.data["iwc_att"]
-    x[iwc_status > 3] = ma.masked
-    testing.assert_array_almost_equal(x, compare[:])
+    obs = ObservationManager("cf", str(obs_file))
+    mask = obs._classify_basic_mask(cat.category_bits)  # type: ignore
+    result = obs._mask_cloud_bits(mask)
+    expected = np.array([[0, 0, 0, 1, 1, 1], [0, 1, 0, 1, 0, 1]])
+    testing.assert_array_almost_equal(result, expected)
+
+
+def test_mask_iwc(obs_file) -> None:
+    obs = ObservationManager("iwc", str(obs_file))
+    iwc_status = obs.getvar("iwc_retrieval_status")
+    expected = ma.copy(obs.getvar("iwc"))
+    expected[~np.isin(iwc_status, (1, 3))] = ma.masked
+    obs._mask_iwc()
+    testing.assert_array_almost_equal(expected, obs.data["iwc"][:])
diff --git a/cloudnetpy/model_evaluation/tests/unit/test_plot_tools.py b/cloudnetpy/model_evaluation/tests/unit/test_plot_tools.py
index 4807dea9..670a5a82 100644
--- a/cloudnetpy/model_evaluation/tests/unit/test_plot_tools.py
+++ b/cloudnetpy/model_evaluation/tests/unit/test_plot_tools.py
@@ -5,177 +5,136 @@
 
 
 def test_parse_wanted_names(regrid_file) -> None:
-    """nc_file: str, name: str, model: str,
-    vars: list | None = None,
-    advance: bool = False
-    """
-    compare = ["ecmwf_cf", "cf_ecmwf"]
-    x, _ = plt.parse_wanted_names(regrid_file, "cf", "ecmwf")
-    assert x == compare
+    expected = ["model_cf", "cf_V"]
+    result, _ = plt.parse_wanted_names(regrid_file, "cf")
+    assert result == expected
 
 
 def test_parse_wanted_names_adv(regrid_file) -> None:
-    compare = ["ecmwf_cf", "cf_adv_ecmwf"]
-    _, x_adv = plt.parse_wanted_names(regrid_file, "cf", "ecmwf")
-    assert x_adv == compare
+    expected = ["model_cf", "cf_V_adv"]
+    _, result_adv = plt.parse_wanted_names(regrid_file, "cf")
+    assert result_adv == expected
 
 
 def test_parse_wanted_names_advance_False(regrid_file) -> None:
-    compare = ["ecmwf_cf", "cf_ecmwf"]
-    x, _ = plt.parse_wanted_names(regrid_file, "cf", "ecmwf", advance=False)
-    assert x == compare
+    expected = ["model_cf", "cf_V"]
+    result, _ = plt.parse_wanted_names(regrid_file, "cf", advance=False)
+    assert result == expected
 
 
 def test_parse_wanted_names_advance_True(regrid_file) -> None:
-    compare = ["ecmwf_cf", "ecmwf_cf_cirrus", "ecmwf_cf_snow", "cf_ecmwf"]
-    x, _ = plt.parse_wanted_names(regrid_file, "cf", "ecmwf", advance=True)
-    assert x == compare
+    expected = ["model_cf", "model_cf_cirrus", "model_cf_snow", "cf_V"]
+    result, _ = plt.parse_wanted_names(regrid_file, "cf", advance=True)
+    assert result == expected
 
 
 def test_parse_wanted_names_adv_advance_True(regrid_file) -> None:
-    compare = ["ecmwf_cf", "ecmwf_cf_cirrus", "ecmwf_cf_snow", "cf_adv_ecmwf"]
-    _, x_adv = plt.parse_wanted_names(regrid_file, "cf", "ecmwf", advance=True)
-    assert x_adv == compare
+    expected = ["model_cf", "model_cf_cirrus", "model_cf_snow", "cf_V_adv"]
+    _, result_adv = plt.parse_wanted_names(regrid_file, "cf", advance=True)
+    assert result_adv == expected
 
 
 def test_parse_wanted_names_adv_advance_False(regrid_file) -> None:
-    compare = ["ecmwf_cf", "cf_adv_ecmwf"]
-    _, x_adv = plt.parse_wanted_names(regrid_file, "cf", "ecmwf", advance=False)
-    assert x_adv == compare
+    expected = ["model_cf", "cf_V_adv"]
+    _, result_adv = plt.parse_wanted_names(regrid_file, "cf", advance=False)
+    assert result_adv == expected
 
 
 def test_parse_wanted_names_fixed_list(regrid_file) -> None:
-    compare = ["ecmwf_cf", "ecmwf_cf_cirrus"]
-    x, _ = plt.parse_wanted_names(regrid_file, "cf", "ecmwf", variables=compare)
-    assert x == compare
+    expected = ["model_cf", "model_cf_cirrus"]
+    result, _ = plt.parse_wanted_names(regrid_file, "cf", variables=expected)
+    assert result == expected
 
 
 def test_parse_wanted_names_adv_fixed_list(regrid_file) -> None:
-    compare = ["ecmwf_cf_cirrus", "cf_adv_ecmwf"]
-    x_adv, _ = plt.parse_wanted_names(regrid_file, "cf", "ecmwf", variables=compare)
-    assert x_adv == compare
+    expected = ["model_cf_cirrus", "cf_V_adv"]
+    result_adv, _ = plt.parse_wanted_names(regrid_file, "cf", variables=expected)
+    assert result_adv == expected
 
 
 def test_sort_model2first_element() -> None:
-    a = ["ec_i", "cf_ec_i", "cf_ec_ii", "ec_ii"]
-    e = "ec"
-    compare = ["ec_i", "ec_ii", "cf_ec_i", "cf_ec_ii"]
-    x = plt.sort_model2first_element(a, e)
-    assert x == compare
-
-
-def test_sort_cycles_vars() -> None:
-    a = ["era5_cf_1-12", "era5_cf_7-18", "cf_era5_1-12", "cf_era5_7-18"]
-    e = "era5"
-    compare = [
-        ["era5_cf_1-12", "cf_era5_1-12"],
-        ["era5_cf_7-18", "cf_era5_7-18"],
-    ]
-    x, _ = plt.sort_cycles(a, e)
-    assert x == compare
-
-
-def test_sort_cycles_simo() -> None:
-    a = ["era5_cf_1-12", "era5_cf_7-18", "cf_era5_1-12", "cf_era5_7-18"]
-    e = "era5"
-    compare = ["1-12", "7-18"]
-    _, y = plt.sort_cycles(a, e)
-    assert y == compare
-
-
-def test_sort_cycles_vars_missing() -> None:
-    a = ["icon_cf_12-23", "icon_cf_36-47", "cf_icon_12-23", "cf_icon_36-47"]
-    e = "icon"
-    compare = [
-        ["icon_cf_12-23", "cf_icon_12-23"],
-        ["icon_cf_36-47", "cf_icon_36-47"],
-    ]
-    x, _ = plt.sort_cycles(a, e)
-    assert x == compare
-
-
-def test_sort_cycles_cycles_missing() -> None:
-    a = ["icon_cf_12-23", "icon_cf_36-47", "cf_icon_12-23", "cf_icon_36-47"]
-    e = "icon"
-    compare = ["12-23", "36-47"]
-    _, y = plt.sort_cycles(a, e)
-    assert y == compare
+    a = ["model_cf", "cf_V", "cf_A", "model_cf_cirrus"]
+    expected = ["model_cf", "model_cf_cirrus", "cf_V", "cf_A"]
+    result = plt.sort_model2first_element(a)
+    assert result == expected
 
 
 def test_read_data_characters(regrid_file) -> None:
-    t = np.array([[2, 2], [6, 6], [10, 10]])
-    h = np.array([[0.01, 0.014], [0.008, 0.014], [0.009, 0.015]])
+    time = np.array([[2, 2], [6, 6], [10, 10]])
+    height = np.array([[0.01, 0.014], [0.008, 0.014], [0.009, 0.015]])
     data = np.array([[0, 2], [3, 6], [5, 8]])
-    x, y, z = plt.read_data_characters(regrid_file, "ecmwf_cf", "ecmwf")
-    compare = [data, t, h]
-    test = [x, y, z]
+    result_data, result_time, result_height = plt.read_data_characters(
+        regrid_file, "model_cf"
+    )
+    expected = [data, time, height]
+    results = [result_data, result_time, result_height]
     for i in range(3):
-        testing.assert_array_almost_equal(compare[i], test[i])
+        testing.assert_array_almost_equal(expected[i], results[i])
 
 
 def test_mask_small_values_lwc() -> None:
     name = "lwc_lol"
     data = ma.array([[0, 1], [3, 6], [5, 8]])
     data = plt.mask_small_values(data, name)
-    compare = ma.array([[0, 1], [3, 6], [5, 8]])
-    compare[0, 0] = ma.masked
-    testing.assert_array_almost_equal(data.mask, compare.mask)
+    expected = ma.array([[0, 1], [3, 6], [5, 8]])
+    expected[0, 0] = ma.masked
+    testing.assert_array_almost_equal(data.mask, expected.mask)
 
 
 def test_mask_small_values_lwc_mask() -> None:
     name = "lwc_lol"
     data = ma.array([[0, 0.000001], [3, 6], [5, 8]])
     data = plt.mask_small_values(data, name)
-    compare = ma.array([[0, -0.000001], [3, 6], [5, 8]])
-    compare[0, 0] = ma.masked
-    compare[0, 1] = ma.masked
-    testing.assert_array_almost_equal(data.mask, compare.mask)
+    expected = ma.array([[0, -0.000001], [3, 6], [5, 8]])
+    expected[0, 0] = ma.masked
+    expected[0, 1] = ma.masked
+    testing.assert_array_almost_equal(data.mask, expected.mask)
 
 
 def test_mask_small_values_iwc() -> None:
     name = "iwc_lol"
     data = ma.array([[0, 1], [3, 6], [5, 8]])
     data = plt.mask_small_values(data, name)
-    compare = ma.array([[0, 1], [3, 6], [5, 8]])
-    compare[0, 0] = ma.masked
-    testing.assert_array_almost_equal(data.mask, compare.mask)
+    expected = ma.array([[0, 1], [3, 6], [5, 8]])
+    expected[0, 0] = ma.masked
+    testing.assert_array_almost_equal(data.mask, expected.mask)
 
 
 def test_mask_small_values_iwc_mask() -> None:
     name = "iwc_lol"
     data = ma.array([[0, 0.00000001], [3, 6], [5, 8]])
     data = plt.mask_small_values(data, name)
-    compare = ma.array([[0, -0.000001], [3, 6], [5, 8]])
-    compare[0, 0] = ma.masked
-    compare[0, 1] = ma.masked
-    testing.assert_array_almost_equal(data.mask, compare.mask)
+    expected = ma.array([[0, -0.000001], [3, 6], [5, 8]])
+    expected[0, 0] = ma.masked
+    expected[0, 1] = ma.masked
+    testing.assert_array_almost_equal(data.mask, expected.mask)
 
 
 def test_mask_small_values() -> None:
     name = "cf_lol"
     data = ma.array([[0, 1], [3, 6], [5, 8]])
     data = plt.mask_small_values(data, name)
-    compare = ma.array([[0, 1], [3, 6], [5, 8]])
-    compare[0, 0] = ma.masked
-    testing.assert_array_almost_equal(data.mask, compare.mask)
+    expected = ma.array([[0, 1], [3, 6], [5, 8]])
+    expected[0, 0] = ma.masked
+    testing.assert_array_almost_equal(data.mask, expected.mask)
 
 
 def test_mask_small_values_mask() -> None:
     name = "cf_lol"
     data = ma.array([[0, -0.000001], [3, 6], [5, 8]])
     data_mask = plt.mask_small_values(data, name)
-    compare = ma.array([[0, -0.000001], [3, 6], [5, 8]])
-    compare[0, 0] = ma.masked
-    compare[0, 1] = ma.masked
-    testing.assert_array_equal(data_mask.mask, compare.mask)
+    expected = ma.array([[0, -0.000001], [3, 6], [5, 8]])
+    expected[0, 0] = ma.masked
+    expected[0, 1] = ma.masked
+    testing.assert_array_equal(data_mask.mask, expected.mask)
 
 
 def test_reshape_1d2nd() -> None:
-    oned = np.array([1, 2, 3, 4])
-    twod = np.array([[0, 0], [0, 0], [0, 0], [0, 0]])
-    compare = np.array([[1, 1], [2, 2], [3, 3], [4, 4]])
-    x = plt.reshape_1d2nd(oned, twod)
-    testing.assert_array_almost_equal(x, compare)
+    one_d = np.array([1, 2, 3, 4])
+    two_d = np.array([[0, 0], [0, 0], [0, 0], [0, 0]])
+    expected = np.array([[1, 1], [2, 2], [3, 3], [4, 4]])
+    result = plt.reshape_1d2nd(one_d, two_d)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_create_segment_values() -> None:
@@ -183,89 +142,89 @@ def test_create_segment_values() -> None:
     model = ma.array(np.ones(model_mask.shape), mask=model_mask)
     obs_mask = np.array([[0, 0, 0, 0], [0, 0, 1, 0], [1, 0, 1, 1]], dtype=bool)
     obs = ma.array(np.ones(obs_mask.shape), mask=obs_mask)
-    x = plt.create_segment_values(model, obs)
-    compare = np.array([[2, 3, 3, 2], [2, 2, 0, 2], [0, 2, 1, 1]])
-    testing.assert_array_almost_equal(x, compare)
+    result = plt.create_segment_values(model, obs)
+    expected = np.array([[2, 3, 3, 2], [2, 2, 0, 2], [0, 2, 1, 1]])
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_rolling_mean() -> None:
     data = np.ma.array([1, 2, 7, 4, 2, 3, 8, 5])
-    x = plt.rolling_mean(data, 2)
-    compare = np.array([1.5, 4.5, 5.5, 3, 2.5, 5.5, 6.5, 5])
-    testing.assert_array_almost_equal(x, compare)
+    result = plt.rolling_mean(data, 2)
+    expected = np.array([1.5, 4.5, 5.5, 3, 2.5, 5.5, 6.5, 5])
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_rolling_mean_nan() -> None:
     data = np.ma.array([1, 2, np.nan, 4, 2, np.nan, 8, 5])
-    x = plt.rolling_mean(data, 2)
-    compare = np.array([1.5, 2, 4, 3, 2, 8, 6.5, 5])
-    testing.assert_array_almost_equal(x, compare)
+    result = plt.rolling_mean(data, 2)
+    expected = np.array([1.5, 2, 4, 3, 2, 8, 6.5, 5])
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_rolling_mean_mask() -> None:
     data = np.ma.array([1, 2, 7, 4, 2, 3, 8, 5])
     data.mask = np.array([0, 0, 1, 0, 1, 0, 0, 1])
-    x = plt.rolling_mean(data, 2)
-    compare = np.array([1.5, 2, 4, 4, 3, 5.5, 8, np.nan])
-    testing.assert_array_almost_equal(x, compare)
+    result = plt.rolling_mean(data, 2)
+    expected = np.array([1.5, 2, 4, 4, 3, 5.5, 8, np.nan])
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_rolling_mean_all_mask() -> None:
     data = np.ma.array([1, 2, 7, 4, 2, 3, 8, 5])
     data.mask = np.array([0, 1, 1, 1, 1, 0, 0, 1])
-    x = plt.rolling_mean(data, 2)
-    compare = np.array([1, np.nan, np.nan, np.nan, 3, 5.5, 8, np.nan])
-    testing.assert_array_almost_equal(x, compare)
+    result = plt.rolling_mean(data, 2)
+    expected = np.array([1, np.nan, np.nan, np.nan, 3, 5.5, 8, np.nan])
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_change2one_dim_axes_maskY() -> None:
-    x = np.ma.array(
+    axis_x = np.ma.array(
         [[1, 1, 1, 1, 1], [2, 2, 2, 2, 2], [3, 3, 3, 3, 3], [4, 4, 4, 4, 4]],
     )
-    y = np.ma.array(
+    axis_y = np.ma.array(
         [[1, 2, 3, 4, 5], [1, 2, 3, 4, 5], [1, 2, 3, 4, 5], [1, 2, 3, 4, 5]],
     )
-    y[1] = np.ma.masked
+    axis_y[1] = np.ma.masked
     data = np.ma.array(
         [[1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1]],
     )
-    x, y, data = plt.change2one_dim_axes(x, y, data)
-    compare_x = np.array([1, 2, 3, 4])
-    compare_y = np.array([1, 2, 3, 4, 5])
-    testing.assert_array_almost_equal(x, compare_x)
-    testing.assert_array_almost_equal(y, compare_y)
+    result_x, result_y, result_data = plt.change2one_dim_axes(axis_x, axis_y, data)
+    expected_x = np.array([1, 2, 3, 4])
+    expected_y = np.array([1, 2, 3, 4, 5])
+    testing.assert_array_almost_equal(result_x, expected_x)
+    testing.assert_array_almost_equal(result_y, expected_y)
 
 
 def test_change2one_dim_axes_maskX() -> None:
-    x = np.ma.array(
+    axis_x = np.ma.array(
         [[1, 1, 1, 1, 1], [2, 2, 2, 2, 2], [3, 3, 3, 3, 3], [4, 4, 4, 4, 4]],
     )
-    y = np.ma.array(
+    axis_y = np.ma.array(
         [[1, 2, 3, 4, 5], [1, 2, 3, 4, 5], [1, 2, 3, 4, 5], [1, 2, 3, 4, 5]],
     )
-    x[1] = np.ma.masked
+    axis_x[1] = np.ma.masked
     data = np.ma.array(
         [[1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1]],
     )
-    x, y, data = plt.change2one_dim_axes(x, y, data)
-    compare_x = np.array([1, 2, 3, 4])
-    compare_y = np.array([1, 2, 3, 4, 5])
-    testing.assert_array_almost_equal(x, compare_x)
-    testing.assert_array_almost_equal(y, compare_y)
+    result_x, result_y, result_data = plt.change2one_dim_axes(axis_x, axis_y, data)
+    expected_x = np.array([1, 2, 3, 4])
+    expected_y = np.array([1, 2, 3, 4, 5])
+    testing.assert_array_almost_equal(result_x, expected_x)
+    testing.assert_array_almost_equal(result_y, expected_y)
 
 
 def test_change2one_dim_axes() -> None:
-    x = np.ma.array(
+    axis_x = np.ma.array(
         [[1, 1, 1, 1, 1], [2, 2, 2, 2, 2], [3, 3, 3, 3, 3], [4, 4, 4, 4, 4]],
     )
-    compare_x = np.copy(x)
-    y = np.ma.array(
+    expected_x = np.copy(axis_x)
+    axis_y = np.ma.array(
         [[1, 2, 3, 4, 5], [1, 2, 3, 4, 5], [1, 2, 3, 4, 5], [1, 2, 3, 4, 5]],
     )
-    compare_y = np.copy(y)
+    expected_y = np.copy(axis_y)
     data = np.ma.array(
         [[1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1]],
     )
-    x, y, data = plt.change2one_dim_axes(x, y, data)
-    testing.assert_array_almost_equal(x, compare_x)
-    testing.assert_array_almost_equal(y, compare_y)
+    result_x, result_y, result_data = plt.change2one_dim_axes(axis_x, axis_y, data)
+    testing.assert_array_almost_equal(result_x, expected_x)
+    testing.assert_array_almost_equal(result_y, expected_y)
diff --git a/cloudnetpy/model_evaluation/tests/unit/test_plotting.py b/cloudnetpy/model_evaluation/tests/unit/test_plotting.py
index f81bf0a7..6d947df6 100644
--- a/cloudnetpy/model_evaluation/tests/unit/test_plotting.py
+++ b/cloudnetpy/model_evaluation/tests/unit/test_plotting.py
@@ -33,40 +33,6 @@ def test_get_cf_title_cycle(key) -> None:
     assert x == value
 
 
-@pytest.mark.parametrize(
-    "key, value",
-    [
-        ("iwc", "Product"),
-        ("iwc_att", "Product with good attenuation"),
-        ("iwc_rain", "Product with rain"),
-        ("iwc_adv", "Product"),
-        ("iwc_att_adv", "Product with good attenuation"),
-        ("iwc_rain_adv", "Product with rain"),
-    ],
-)
-def test_get_iwc_title(key, value) -> None:
-    field_name = key + "_" + MODEL
-    x = pl._get_iwc_title(field_name, VARIABLE_INFO)
-    assert x == value
-
-
-@pytest.mark.parametrize(
-    "key, value",
-    [
-        ("iwc", "Product"),
-        ("iwc_att", "Product with good attenuation"),
-        ("iwc_rain", "Product with rain"),
-        ("iwc_adv", "Product"),
-        ("iwc_att_adv", "Product with good attenuation"),
-        ("iwc_rain_adv", "Product with rain"),
-    ],
-)
-def test_get_iwc_title_cycle(key, value) -> None:
-    field_name = key + "_" + MODEL + "_001"
-    x = pl._get_iwc_title(field_name, VARIABLE_INFO)
-    assert x == value
-
-
 def test_get_product_title() -> None:
     value = "Product"
     x = pl._get_product_title(VARIABLE_INFO)
@@ -99,20 +65,6 @@ def test_get_cf_title_stat(key) -> None:
     assert x == value
 
 
-@pytest.mark.parametrize(
-    "key, value",
-    [
-        ("iwc", "Product"),
-        ("iwc_att", "Product with good attenuation"),
-        ("iwc_rain", "Product with rain"),
-    ],
-)
-def test_get_iwc_title_stat(key, value) -> None:
-    field_name = key + "_" + MODEL
-    x = pl._get_iwc_title_stat(field_name, VARIABLE_INFO)
-    assert x == value
-
-
 @pytest.mark.parametrize("key", ["lwc"])
 def test_get_product_title_stat(key) -> None:
     x = pl._get_product_title_stat(VARIABLE_INFO)
diff --git a/cloudnetpy/model_evaluation/tests/unit/test_statistical_methods.py b/cloudnetpy/model_evaluation/tests/unit/test_statistical_methods.py
index dac6ec4e..8f1eef6b 100644
--- a/cloudnetpy/model_evaluation/tests/unit/test_statistical_methods.py
+++ b/cloudnetpy/model_evaluation/tests/unit/test_statistical_methods.py
@@ -11,9 +11,9 @@
 def test_relative_error() -> None:
     model = ma.array([[1, 2, 2, 3], [2, 4, 10, 1]])
     observation = ma.array([[3, 2, 5, 4], [4, 6, 8, 4]])
-    x, _ = sts.relative_error(model, observation)
-    compare = ma.array([[-66.67, 0.0, -60.0, -25.0], [-50.0, -33.33, 25.0, -75.0]])
-    testing.assert_array_almost_equal(x, compare)
+    result_model, _ = sts.relative_error(model, observation)
+    expected = ma.array([[-66.67, 0.0, -60.0, -25.0], [-50.0, -33.33, 25.0, -75.0]])
+    testing.assert_array_almost_equal(result_model, expected)
 
 
 def test_relative_error_mask() -> None:
@@ -21,48 +21,48 @@ def test_relative_error_mask() -> None:
     model.mask = np.array([[0, 0, 0, 1], [0, 1, 0, 1]])
     observation = ma.array([[1, 2, 5, 4], [4, 6, 8, 1]])
     observation.mask = np.array([[1, 0, 0, 1], [0, 0, 0, 1]])
-    x, _ = sts.relative_error(model, observation)
-    compare = ma.array(
+    result_model, _ = sts.relative_error(model, observation)
+    expected = ma.array(
         [
             ma.array([-99, 0.0, -60.0, -99], mask=[True, False, False, True]),
             ma.array([-50.0, -83.33, 25.0, -99], mask=[False, False, False, True]),
         ],
     )
-    testing.assert_array_almost_equal(x, compare)
+    testing.assert_array_almost_equal(result_model, expected)
 
 
 def test_relative_error_nan() -> None:
     model = ma.array([[1, 2, 2, 3], [2, 4, 10, 1]])
     observation = ma.array([[3, 2, 5, np.nan], [4, np.nan, 8, 4]])
-    x, _ = sts.relative_error(model, observation)
-    compare = ma.array(
+    result_model, _ = sts.relative_error(model, observation)
+    expected = ma.array(
         [
             ma.array([-66.67, 0.0, -60.0, -99], mask=[False, False, False, True]),
             ma.array([-50.0, -99, 25.0, -75.0], mask=[False, True, False, False]),
         ],
     )
-    testing.assert_array_almost_equal(x, compare)
+    testing.assert_array_almost_equal(result_model, expected)
 
 
 def test_absolute_error() -> None:
     model = ma.array([[0.1, 0.2, 0.2, 0.3], [0.2, 0.4, 1.0, 0.0]])
     observation = ma.array([[0.2, 0.2, 0.1, 0.4], [0.4, 0.6, 0.8, 0.2]])
-    x, _ = sts.absolute_error(model, observation)
-    compare = ma.array([[10.0, 0.0, -10.0, 10.0], [20.0, 20.0, -20.0, 20.0]])
-    testing.assert_array_almost_equal(x, compare)
+    result_model, _ = sts.absolute_error(model, observation)
+    expected = ma.array([[10.0, 0.0, -10.0, 10.0], [20.0, 20.0, -20.0, 20.0]])
+    testing.assert_array_almost_equal(result_model, expected)
 
 
 def test_absolute_error_nan() -> None:
     model = ma.array([[0.1, 0.2, 0.2, 0.3], [0.2, 0.4, 1.0, 0.0]])
     observation = ma.array([[0.2, np.nan, 0.1, 0.4], [np.nan, 0.6, 0.8, 0.2]])
-    x, _ = sts.absolute_error(model, observation)
-    compare = ma.array(
+    result_model, _ = sts.absolute_error(model, observation)
+    expected = ma.array(
         [
             ma.array([10.0, -99, -10.0, 10.0], mask=[False, True, False, False]),
             ma.array([-99, 20.0, -20.0, 20.0], mask=[True, False, False, False]),
         ],
     )
-    testing.assert_array_almost_equal(x, compare)
+    testing.assert_array_almost_equal(result_model, expected)
 
 
 def test_absolute_error_mask() -> None:
@@ -70,46 +70,46 @@ def test_absolute_error_mask() -> None:
     model.mask = np.array([[0, 0, 0, 1], [0, 1, 0, 1]])
     observation = ma.array([[0.2, 0.2, 0.1, 0.4], [0.4, 0.6, 0.8, 0.2]])
     observation.mask = np.array([[0, 0, 0, 0], [0, 1, 0, 0]])
-    x, _ = sts.absolute_error(model, observation)
-    compare = ma.array(
+    result_model, _ = sts.absolute_error(model, observation)
+    expected = ma.array(
         [
             ma.array([10.0, 0.0, -10.0, -99], mask=[False, False, False, True]),
             ma.array([20.0, -99, -20.0, -99], mask=[False, True, False, True]),
         ],
     )
-    testing.assert_array_almost_equal(x, compare)
+    testing.assert_array_almost_equal(result_model, expected)
 
 
 def test_combine_masked_indices() -> None:
     model = ma.array([[1, 2, 2, 3], [2, 4, 10, 1]])
     observation = ma.array([[3, 2, 5, 4], [4, 6, 8, 4]])
-    x, y = sts.combine_masked_indices(model, observation)
-    compare_m = ma.array([[4, 2, 2, 3], [2, 4, 10, 6]])
-    compare_o = ma.array([[3, 2, 5, 4], [4, 6, 8, 4]])
-    testing.assert_array_almost_equal(x, compare_m)
-    testing.assert_array_almost_equal(y, compare_o)
+    result_model, result_obs = sts.combine_masked_indices(model, observation)
+    expected_model = ma.array([[4, 2, 2, 3], [2, 4, 10, 6]])
+    expected_obs = ma.array([[3, 2, 5, 4], [4, 6, 8, 4]])
+    testing.assert_array_almost_equal(result_model, expected_model)
+    testing.assert_array_almost_equal(result_obs, expected_obs)
 
 
 def test_combine_masked_indices_min() -> None:
     model = ma.array([[1, 2, 2, 3], [2, 4, 10, 1]])
     observation = ma.array([[3, 2, 5, 4], [4, 6, 8, 4]])
-    x, y = sts.combine_masked_indices(model, observation)
+    result_model, result_obs = sts.combine_masked_indices(model, observation)
 
-    compare_m = ma.array(
+    expected_model = ma.array(
         [
             ma.array([-99, 2, 2, 3], mask=[True, False, False, False]),
             ma.array([2, 4, 10, -99], mask=[True, False, False, True]),
         ],
     )
-    compare_o = ma.array(
+    expected_obs = ma.array(
         [
             ma.array([-99, 2, 5, 4], mask=[True, False, False, False]),
             ma.array([4, 6, 8, -99], mask=[True, False, False, True]),
         ],
     )
 
-    testing.assert_array_almost_equal(x, compare_m)
-    testing.assert_array_almost_equal(y, compare_o)
+    testing.assert_array_almost_equal(result_model, expected_model)
+    testing.assert_array_almost_equal(result_obs, expected_obs)
 
 
 def test_combine_masked_indices_mask() -> None:
@@ -117,42 +117,42 @@ def test_combine_masked_indices_mask() -> None:
     model.mask = ma.array([[1, 0, 0, 0], [0, 0, 0, 1]])
     observation = ma.array([[3, 2, 1, 4], [4, 6, 8, 4]])
     observation.mask = ma.array([[0, 1, 1, 0], [1, 0, 0, 0]])
-    x, y = sts.combine_masked_indices(model, observation)
-    model = ma.array(
+    result_model, result_obs = sts.combine_masked_indices(model, observation)
+    expected_model = ma.array(
         [
             ma.array([-99, -99, -99, 3], mask=[True, True, True, False]),
             ma.array([-99, 4, 10, -99], mask=[True, False, False, True]),
         ],
     )
-    observation = ma.array(
+    expected_obs = ma.array(
         [
             ma.array([-99, -99, -99, 4], mask=[True, True, True, False]),
             ma.array([-99, 6, 8, -99], mask=[True, False, False, True]),
         ],
     )
-    testing.assert_array_almost_equal(x, model)
-    testing.assert_array_almost_equal(y, observation)
+    testing.assert_array_almost_equal(result_model, expected_model)
+    testing.assert_array_almost_equal(result_obs, expected_obs)
 
 
 def test_combine_masked_indices_nan() -> None:
     model = ma.array([[1, 2, 2, 3], [2, 4, 10, 1]])
     observation = ma.array([[np.nan, 2, 5, 4], [4, 6, np.nan, 4]])
-    x, y = sts.combine_masked_indices(model, observation)
-    model = ma.array(
+    result_model, result_obs = sts.combine_masked_indices(model, observation)
+    expected_model = ma.array(
         [
             ma.array([-99, 2, 2, 3], mask=[True, False, False, False]),
             ma.array([2, 4, -99, 1], mask=[False, False, True, False]),
         ],
     )
-    testing.assert_array_almost_equal(x, model)
-    testing.assert_array_almost_equal(y, observation)
+    testing.assert_array_almost_equal(result_model, expected_model)
+    testing.assert_array_almost_equal(result_obs, observation)
 
 
 def test_calc_common_area_sum() -> None:
     model = ma.array([[1, 2, 2, 3], [2, 4, 10, 1]])
     observation = ma.array([[3, 2, 1, 4], [4, 6, 8, 4]])
-    x, _ = sts.calc_common_area_sum(model, observation)
-    testing.assert_almost_equal(x, 100)
+    result, _ = sts.calc_common_area_sum(model, observation)
+    testing.assert_almost_equal(result, 100)
 
 
 def test_calc_common_area_sum_min() -> None:
@@ -160,8 +160,8 @@ def test_calc_common_area_sum_min() -> None:
     model.mask = ma.array([[1, 0, 0, 0], [0, 0, 0, 1]])
     observation = ma.array([[3, 2, 1, 4], [4, 6, 8, 4]])
     observation.mask = ma.array([[0, 1, 1, 0], [1, 0, 0, 0]])
-    x, _ = sts.calc_common_area_sum(model, observation)
-    testing.assert_almost_equal(x, 60.0)
+    result, _ = sts.calc_common_area_sum(model, observation)
+    testing.assert_almost_equal(result, 60.0)
 
 
 def test_calc_common_area_sum_nan() -> None:
@@ -169,8 +169,8 @@ def test_calc_common_area_sum_nan() -> None:
     model.mask = ma.array([[1, 0, 0, 0], [0, 0, 0, 1]])
     observation = ma.array([[3, 2, 1, 4], [4, 6, np.nan, 4]])
     observation.mask = ma.array([[0, 0, 1, 0], [1, 0, 0, 0]])
-    x, _ = sts.calc_common_area_sum(model, observation)
-    testing.assert_almost_equal(x, 25.0)
+    result, _ = sts.calc_common_area_sum(model, observation)
+    testing.assert_almost_equal(result, 25.0)
 
 
 def test_calc_common_area_sum_mask() -> None:
@@ -178,18 +178,18 @@ def test_calc_common_area_sum_mask() -> None:
     model.mask = ma.array([[1, 0, 0, 0], [0, 0, 0, 0]])
     observation = ma.array([[3, 2, 1, 4], [4, 6, 8, 4]])
     observation.mask = ma.array([[0, 0, 1, 0], [1, 0, 0, 0]])
-    x, _ = sts.calc_common_area_sum(model, observation)
-    testing.assert_almost_equal(x, 50.0)
+    result, _ = sts.calc_common_area_sum(model, observation)
+    testing.assert_almost_equal(result, 50.0)
 
 
 def test_histogram() -> None:
     model = ma.array([[1, 2, 2, 3], [2, 4, 10, 1]])
     observation = ma.array([[3, 2, 1, 4], [4, 6, 8, 4]])
-    compare_x = np.array([1, 2, 2, 3, 2, 4, 8, 1])
-    compare_y = np.array([3, 2, 1, 4, 4, 6, 8, 4])
-    x, y = sts.histogram(PRODUCT_iwc, model, observation)
-    testing.assert_array_almost_equal(x, compare_x)
-    testing.assert_array_almost_equal(y, compare_y)
+    expected_x = np.array([1, 2, 2, 3, 2, 4, 8, 1])
+    expected_y = np.array([3, 2, 1, 4, 4, 6, 8, 4])
+    result_x, result_y = sts.histogram(PRODUCT_iwc, model, observation)
+    testing.assert_array_almost_equal(result_x, expected_x)
+    testing.assert_array_almost_equal(result_y, expected_y)
 
 
 def test_histogram_mask() -> None:
@@ -197,41 +197,41 @@ def test_histogram_mask() -> None:
     model.mask = ma.array([[1, 0, 0, 0], [0, 0, 0, 1]])
     observation = ma.array([[3, 2, 1, 4], [4, 6, 8, 4]])
     observation.mask = ma.array([[0, 1, 1, 0], [1, 0, 0, 0]])
-    compare_x = ma.array([2, 2, 3, 2, 4, 8])
-    compare_y = ma.array([3, 4, 6, 8, 4])
-    x, y = sts.histogram(PRODUCT_iwc, model, observation)
-    testing.assert_array_almost_equal(x, compare_x)
-    testing.assert_array_almost_equal(y, compare_y)
+    expected_x = ma.array([2, 2, 3, 2, 4, 8])
+    expected_y = ma.array([3, 4, 6, 8, 4])
+    result_x, result_y = sts.histogram(PRODUCT_iwc, model, observation)
+    testing.assert_array_almost_equal(result_x, expected_x)
+    testing.assert_array_almost_equal(result_y, expected_y)
 
 
 def test_histogram_nan() -> None:
     model = ma.array([[1, 2, 2, 3], [2, 4, 10, 1]])
     observation = ma.array([[3, np.nan, 1, 4], [np.nan, 6, np.nan, 4]])
-    compare_x = np.array([1, 2, 2, 3, 2, 4, 6, 1])
-    compare_y = np.array([3, 1, 4, 6, 4])
-    x, y = sts.histogram(PRODUCT_iwc, model, observation)
-    testing.assert_array_almost_equal(x, compare_x)
-    testing.assert_array_almost_equal(y, compare_y)
+    expected_x = np.array([1, 2, 2, 3, 2, 4, 6, 1])
+    expected_y = np.array([3, 1, 4, 6, 4])
+    result_x, result_y = sts.histogram(PRODUCT_iwc, model, observation)
+    testing.assert_array_almost_equal(result_x, expected_x)
+    testing.assert_array_almost_equal(result_y, expected_y)
 
 
 def test_vertical_profile() -> None:
     model = ma.array([[0, 2, 2, 3], [2, 4, 10, 1], [4, 6, 6, 8]])
     observation = ma.array([[3, 1, 1, 4], [4, 3, 8, 0], [5, 5, 3, 2]])
-    x, y = sts.vertical_profile(model, observation)
-    model = ma.array([2, 4, 6, 4])
-    observation = ma.array([4, 3, 4, 2])
-    testing.assert_array_almost_equal(x, model)
-    testing.assert_array_almost_equal(y, observation)
+    result_model, result_obs = sts.vertical_profile(model, observation)
+    expected_model = ma.array([2, 4, 6, 4])
+    expected_obs = ma.array([4, 3, 4, 2])
+    testing.assert_array_almost_equal(result_model, expected_model)
+    testing.assert_array_almost_equal(result_obs, expected_obs)
 
 
 def test_vertical_profile_nan() -> None:
     model = ma.array([[0, 2, 2, 3], [2, 4, 10, 1], [4, 6, 6, 8]])
     observation = ma.array([[3, 1, 1, np.nan], [np.nan, 3, 8, 0], [5, 5, 3, 2]])
-    x, y = sts.vertical_profile(model, observation)
-    model = ma.array([2, 4, 6, 4])
-    observation = ma.array([4, 3, 4, 1])
-    testing.assert_array_almost_equal(x, model)
-    testing.assert_array_almost_equal(y, observation)
+    result_model, result_obs = sts.vertical_profile(model, observation)
+    expected_model = ma.array([2, 4, 6, 4])
+    expected_obs = ma.array([4, 3, 4, 1])
+    testing.assert_array_almost_equal(result_model, expected_model)
+    testing.assert_array_almost_equal(result_obs, expected_obs)
 
 
 def test_vertical_profile_mask() -> None:
@@ -239,11 +239,11 @@ def test_vertical_profile_mask() -> None:
     model.mask = ma.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
     observation = ma.array([[3, 1, 1, 4], [4, 3, 8, 0], [5, 5, 3, 2]])
     observation.mask = ma.array([[0, 1, 0, 0], [0, 1, 0, 1], [0, 1, 0, 0]])
-    x, y = sts.vertical_profile(model, observation)
-    model = ma.array([3, 4, 6, 5.5])
-    observation = ma.array([4, -99, 4, 3], mask=[False, True, False, False])
-    testing.assert_array_almost_equal(x, model)
-    testing.assert_array_almost_equal(y, observation)
+    result_model, result_obs = sts.vertical_profile(model, observation)
+    expected_model = ma.array([3, 4, 6, 5.5])
+    expected_obs = ma.array([4, -99, 4, 3], mask=[False, True, False, False])
+    testing.assert_array_almost_equal(result_model, expected_model)
+    testing.assert_array_almost_equal(result_obs, expected_obs)
 
 
 @pytest.mark.parametrize(
@@ -254,5 +254,5 @@ def test_vertical_profile_mask() -> None:
     ],
 )
 def test_day_stat_title(method, title) -> None:
-    x = sts.day_stat_title(method, PRODUCT_cf)
-    assert x == title
+    result = sts.day_stat_title(method, PRODUCT_cf)
+    assert result == title
diff --git a/cloudnetpy/model_evaluation/tests/unit/test_tools.py b/cloudnetpy/model_evaluation/tests/unit/test_tools.py
index a24f13aa..4040c76c 100644
--- a/cloudnetpy/model_evaluation/tests/unit/test_tools.py
+++ b/cloudnetpy/model_evaluation/tests/unit/test_tools.py
@@ -8,119 +8,123 @@
 from cloudnetpy.model_evaluation.products.model_products import ModelManager
 
 MODEL = "ecmwf"
-OUTPUT_FILE = "/"
 PRODUCT = "iwc"
 
 
-def test_model_file_list() -> None:
-    name = "ec"
-    models = ["00_ec_1", "00_ec_2", "00_ec_3"]
-    tools.check_model_file_list(name, models)
-
-
-def test_model_file_list_fail() -> None:
-    name = "ec"
-    models = ["00_ec_1", "ac_1", "00_ec_2", "00_ec_3"]
-    with pytest.raises(AttributeError):
-        tools.check_model_file_list(name, models)
-
-
 def test_time2datetime() -> None:
-    time_list = np.array(range(10))
-    d = datetime(2020, 4, 7, 0, 0, 0, tzinfo=timezone.utc)
-    x = tools.time2datetime(time_list, d)
-    compare = [
-        datetime(2020, 4, 7, 0, 0, 0, tzinfo=timezone.utc) + timedelta(hours=1 * x)
-        for x in range(10)
+    time = np.array(range(10))
+    date = datetime(2020, 4, 7, 0, 0, 0, tzinfo=timezone.utc)
+    result = tools.time2datetime(time, date)
+    expected = [
+        datetime(2020, 4, 7, 0, 0, 0, tzinfo=timezone.utc) + timedelta(hours=1 * i)
+        for i in range(10)
     ]
-    assert all(a == b for a, b in zip(x, compare))
+    assert all(a == b for a, b in zip(result, expected))
 
 
 def test_rebin_edges() -> None:
     data = np.array([1, 3, 6, 10, 15, 21, 28])
-    compare = np.array([-1, 2, 4.5, 8, 12.5, 18, 24.5, 35])
-    x = tools.rebin_edges(data)
-    testing.assert_array_almost_equal(x, compare)
+    expected = np.array([-1, 2, 4.5, 8, 12.5, 18, 24.5, 35])
+    result = tools.rebin_edges(data)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_calculate_advection_time_hour(model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    h = obj.resolution_h
-    v = np.array([[1, 2, 3], [2, 3, 4], [3, 4, 5]])
-    s = 1
-    compare = h * 1000 / v / 60**2
-    compare[compare > 1 / s] = 1 / s
-    compare = np.asarray([[timedelta(hours=float(t)) for t in tt] for tt in compare])
-    x = tools.calculate_advection_time(int(h), ma.masked_array(v), s)
-    assert x.all() == compare.all()
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    resolution = model.resolution_h
+    wind = np.array([[1, 2, 3], [2, 3, 4], [3, 4, 5]])
+    sampling = 1
+    expected = resolution * 1000 / wind / 60**2
+    expected[expected > 1 / sampling] = 1 / sampling
+    expected = np.asarray([[timedelta(hours=float(t)) for t in tt] for tt in expected])
+    result = tools.calculate_advection_time(resolution, ma.masked_array(wind), sampling)
+    assert result.all() == expected.all()
 
 
 def test_calculate_advection_time_10min(model_file) -> None:
-    obj = ModelManager(str(model_file), MODEL, OUTPUT_FILE, PRODUCT)
-    h = obj.resolution_h
-    v = np.array([[1, 2, 3], [2, 3, 4], [3, 4, 5]])
-    s = 6
-    compare = h * 1000 / v / 60**2
-    compare[compare > 1 / s] = 1 / s
-    compare = np.asarray([[timedelta(hours=float(t)) for t in tt] for tt in compare])
-    x = tools.calculate_advection_time(int(h), ma.masked_array(v), s)
-    assert x.all() == compare.all()
+    model = ModelManager(str(model_file), MODEL, PRODUCT)
+    resolution = model.resolution_h
+    wind = np.array([[1, 2, 3], [2, 3, 4], [3, 4, 5]])
+    sampling = 6
+    expected = resolution * 1000 / wind / 60**2
+    expected[expected > 1 / sampling] = 1 / sampling
+    expected = np.asarray([[timedelta(hours=float(t)) for t in tt] for tt in expected])
+    result = tools.calculate_advection_time(resolution, ma.masked_array(wind), sampling)
+    assert result.all() == expected.all()
+
+
+def test_calculate_advection_time_fractional_resolution() -> None:
+    # A sub-kilometre / fractional resolution must not be truncated to int.
+    resolution = 0.5
+    wind = ma.masked_array([[2.0]])
+    sampling = 6
+    result = tools.calculate_advection_time(resolution, wind, sampling)
+    expected = timedelta(hours=resolution * 1000 / 2.0 / 60**2)
+    assert result[0, 0] == expected
 
 
 def test_get_1d_indices() -> None:
-    w = (1, 5)
-    d = ma.array([0, 1, 2, 3, 4, 5, 6, 7])
-    compare = ma.array([0, 1, 1, 1, 1, 0, 0, 0])
-    x = tools.get_1d_indices(w, d)
-    testing.assert_array_almost_equal(x, compare)
+    window = (1, 5)
+    data = ma.array([0, 1, 2, 3, 4, 5, 6, 7])
+    expected = ma.array([0, 1, 1, 1, 1, 0, 0, 0])
+    result = tools.get_1d_indices(window, data)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_get_1d_indices_mask() -> None:
-    w = (1, 5)
-    d = ma.array([0, 1, 2, 3, 4, 5, 6, 7])
-    m = np.array([0, 0, 1, 0, 0, 1, 0, 1], dtype=bool)
-    x = tools.get_1d_indices(w, d, m)
-    d[m] = ma.masked
-    compare = ma.array(
+    window = (1, 5)
+    data = ma.array([0, 1, 2, 3, 4, 5, 6, 7])
+    mask = np.array([0, 0, 1, 0, 0, 1, 0, 1], dtype=bool)
+    result = tools.get_1d_indices(window, data, mask)
+    data[mask] = ma.masked
+    expected = ma.array(
         [0, 1, -99, 1, 1, -99, 0, -99],
         mask=[False, False, True, False, False, True, False, True],
     )
-    testing.assert_array_almost_equal(x, compare)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_get_adv_indices() -> None:
-    mt = 3
-    at = 4
-    d = ma.array([0, 1, 2, 3, 4, 5, 6, 7])
-    compare: ma.MaskedArray = ma.array([0, 1, 1, 1, 1, 0, 0, 0], dtype=bool)
-    x = tools.get_adv_indices(mt, at, d)
-    testing.assert_array_almost_equal(x, compare)
+    model_t = 3
+    adv_t = 4
+    data = ma.array([0, 1, 2, 3, 4, 5, 6, 7])
+    expected: ma.MaskedArray = ma.array([0, 1, 1, 1, 1, 0, 0, 0], dtype=bool)
+    result = tools.get_adv_indices(model_t, adv_t, data)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_get_adv_indices_mask() -> None:
-    mt = 3
-    at = 4
-    d: ma.MaskedArray = ma.array([0, 1, 2, 3, 4, 5, 6, 7])
-    m: ma.MaskedArray = ma.array([0, 0, 1, 0, 0, 1, 0, 1], dtype=bool)
-    x = tools.get_adv_indices(mt, at, d, m)
-    d[m] = ma.masked
-    compare = ma.array(
+    model_t = 3
+    adv_t = 4
+    data: ma.MaskedArray = ma.array([0, 1, 2, 3, 4, 5, 6, 7])
+    mask: ma.MaskedArray = ma.array([0, 0, 1, 0, 0, 1, 0, 1], dtype=bool)
+    result = tools.get_adv_indices(model_t, adv_t, data, mask)
+    data[mask] = ma.masked
+    expected = ma.array(
         [0, 1, -99, 1, 1, -99, 0, 0],
         mask=[False, False, True, False, False, True, False, False],
     )
-    testing.assert_array_almost_equal(x, compare)
+    testing.assert_array_almost_equal(result, expected)
 
 
 def test_obs_windows_size() -> None:
-    i = np.array([0, 0, 1, 1, 1, 1, 0], dtype=bool)
-    j = np.array([0, 1, 1, 1, 0, 0, 0], dtype=bool)
-    x = tools.get_obs_window_size(i, j)
-    assert x is not None
-    testing.assert_almost_equal(x, (4, 3))
+    ind_x = np.array([0, 0, 1, 1, 1, 1, 0], dtype=bool)
+    ind_y = np.array([0, 1, 1, 1, 0, 0, 0], dtype=bool)
+    result = tools.get_obs_window_size(ind_x, ind_y)
+    assert result is not None
+    testing.assert_almost_equal(result, (4, 3))
 
 
 def test_obs_windows_size_none() -> None:
-    i = np.array([0, 0, 1, 1, 1, 1, 0], dtype=bool)
-    j = np.array([0, 0, 0, 0, 0, 0, 0], dtype=bool)
-    x = tools.get_obs_window_size(i, j)
-    assert x is None
+    ind_x = np.array([0, 0, 1, 1, 1, 1, 0], dtype=bool)
+    ind_y = np.array([0, 0, 0, 0, 0, 0, 0], dtype=bool)
+    result = tools.get_obs_window_size(ind_x, ind_y)
+    assert result is None
+
+
+def test_obs_windows_size_first_index() -> None:
+    ind_x = np.array([1, 0, 0, 0], dtype=bool)
+    ind_y = np.array([1, 0, 0, 0], dtype=bool)
+    result = tools.get_obs_window_size(ind_x, ind_y)
+    assert result is not None
+    testing.assert_almost_equal(result, (1, 1))
diff --git a/cloudnetpy/plotting/plotting.py b/cloudnetpy/plotting/plotting.py
index d7b9fc5f..15caaa97 100644
--- a/cloudnetpy/plotting/plotting.py
+++ b/cloudnetpy/plotting/plotting.py
@@ -153,14 +153,7 @@ def initialize_figure(self) -> tuple[Figure, list[Axes]]:
         return fig, axes_list
 
     def add_subtitle(self, fig: Figure) -> None:
-        fig.suptitle(
-            self._get_subtitle_text(),
-            fontsize=13,
-            y=0.885,
-            x=0.07,
-            horizontalalignment="left",
-            verticalalignment="bottom",
-        )
+        add_subtitle(fig, self.datetime().date(), self.file.location)
 
     def datetime(self) -> datetime:
         return datetime(
@@ -170,13 +163,6 @@ def datetime(self) -> datetime:
             tzinfo=timezone.utc,
         )
 
-    def _get_subtitle_text(self) -> str:
-        measurement_date = date(
-            int(self.file.year), int(self.file.month), int(self.file.day)
-        )
-        site_name = self.file.location.replace("-", " ")
-        return f"{site_name}, {measurement_date.strftime('%d %b %Y').lstrip('0')}"
-
     def _get_valid_variables_and_indices(
         self, requested_variables: list[str]
     ) -> tuple[list[netCDF4.Variable], list[int | None]]:
@@ -291,12 +277,9 @@ def set_xax(self, figure_data: FigureData) -> None:
         if self.file_type in ("cpr-validation", "cpr-tc-validation"):
             return
         resolution = 4
-        x_tick_labels = [
-            f"{int(i):02d}:00"
-            if (24 >= i >= 0 if self.options.edge_tick_labels else 24 > i > 0)
-            else ""
-            for i in np.arange(0, 24.01, resolution)
-        ]
+        x_tick_labels = get_time_tick_labels(
+            resolution, edge_tick_labels=self.options.edge_tick_labels
+        )
         self.ax.set_xticks(np.arange(0, 25, resolution, dtype=int))
         self.ax.set_xticklabels(x_tick_labels, fontsize=12)
         if self.options.minor_ticks:
@@ -1207,6 +1190,29 @@ def generate_figure(
             plt.close(fig)
 
 
+def add_subtitle(fig: Figure, case_date: date, site_name: str) -> None:
+    """Adds a date/site subtitle to a figure."""
+    text = f"{site_name}, {case_date.strftime('%d %b %Y').lstrip('0')}"
+    fig.suptitle(
+        text,
+        fontsize=13,
+        y=0.885,
+        x=0.07,
+        horizontalalignment="left",
+        verticalalignment="bottom",
+    )
+
+
+def get_time_tick_labels(
+    resolution: int = 4, *, edge_tick_labels: bool = False
+) -> list[str]:
+    """Returns hourly tick labels for a 0-24 h time axis."""
+    return [
+        f"{int(i):02d}:00" if (24 >= i >= 0 if edge_tick_labels else 24 > i > 0) else ""
+        for i in np.arange(0, 24.01, resolution)
+    ]
+
+
 def lin2log(*args: npt.ArrayLike) -> list[ma.MaskedArray]:
     return [ma.log10(x) for x in args]
 
diff --git a/cloudnetpy/products/iwc.py b/cloudnetpy/products/iwc.py
index 373f6b08..b9cca5d9 100644
--- a/cloudnetpy/products/iwc.py
+++ b/cloudnetpy/products/iwc.py
@@ -85,7 +85,7 @@ def _calc_random_error() -> npt.NDArray:
             return self.getvar("Z_error") * scaled_temperature * 10
 
         def _calc_error_in_uncorrected_ice() -> float:
-            spec_liq_atten = 1.0 if self.wl_band == 0 else 4.5
+            spec_liq_atten = 1.0 if self.wl_band == "Ka" else 4.5
             liq_atten_scaled = spec_liq_atten * self.coefficients.Z
             return lwp_prior * G_TO_KG * liq_atten_scaled * 2 * 10
 
diff --git a/cloudnetpy/products/product_tools.py b/cloudnetpy/products/product_tools.py
index c6a6acbf..2d76ade4 100644
--- a/cloudnetpy/products/product_tools.py
+++ b/cloudnetpy/products/product_tools.py
@@ -25,6 +25,44 @@ class IceCoefficients(NamedTuple):
     c: float
 
 
+def get_ice_coefficients(product: str, wl_band: str) -> IceCoefficients:
+    """Returns coefficients for ice retrieval.
+
+    References:
+        Hogan et.al. 2006, https://doi.org/10.1175/JAM2340.1
+    """
+    msg = f"Unsupported band: {wl_band}"
+    if product == "ier":
+        if wl_band == "Ka":
+            return IceCoefficients(0.878, -0.000205, -0.0015, 0.0016, -1.52)
+        if wl_band == "W":
+            return IceCoefficients(0.669, -0.000296, -0.00193, -0.000, -1.502)
+        raise ValueError(msg)
+    if wl_band == "Ka":
+        return IceCoefficients(0.878, 0.000242, -0.0186, 0.0699, -1.63)
+    if wl_band == "W":
+        return IceCoefficients(0.669, 0.000580, -0.00706, 0.0923, -0.992)
+    raise ValueError(msg)
+
+
+def z_to_iwc(
+    coefficients: IceCoefficients,
+    z: npt.NDArray,
+    temperature: npt.NDArray,
+) -> npt.NDArray:
+    """Converts reflectivity to ice water content with the Z-T method.
+
+    References:
+        Hogan et.al. 2006, https://doi.org/10.1175/JAM2340.1
+    """
+    return 10 ** (
+        coefficients.ZT * z * temperature
+        + coefficients.T * temperature
+        + coefficients.Z * z
+        + coefficients.c
+    )
+
+
 @dataclass
 class CategoryBits:
     droplet: NDArray[np.bool_]
@@ -214,23 +252,7 @@ def append_status(self, ice_classification: IceClassification) -> None:
         self.append_data(retrieval_status, f"{self.product}_retrieval_status")
 
     def _get_coefficients(self) -> IceCoefficients:
-        """Returns coefficients for ice effective radius retrieval.
-
-        References:
-            Hogan et.al. 2006, https://doi.org/10.1175/JAM2340.1
-        """
-        msg = f"Unsupported band: {self.wl_band}"
-        if self.product == "ier":
-            if self.wl_band == "Ka":
-                return IceCoefficients(0.878, -0.000205, -0.0015, 0.0016, -1.52)
-            if self.wl_band == "W":
-                return IceCoefficients(0.669, -0.000296, -0.00193, -0.000, -1.502)
-            raise ValueError(msg)
-        if self.wl_band == "Ka":
-            return IceCoefficients(0.878, 0.000242, -0.0186, 0.0699, -1.63)
-        if self.wl_band == "W":
-            return IceCoefficients(0.669, 0.000580, -0.00706, 0.0923, -0.992)
-        raise ValueError(msg)
+        return get_ice_coefficients(self.product, self.wl_band)
 
     def _convert_z(self, z_variable: str = "Z") -> npt.NDArray:
         """Calculates temperature weighted z, i.e. ice effective radius [m]."""
@@ -249,16 +271,7 @@ def _convert_z(self, z_variable: str = "Z") -> npt.NDArray:
         scale = (
             g_to_kg if self.product == "iwc" else 3 / (2 * constants.RHO_ICE) * m_to_mu
         )
-        return (
-            10
-            ** (
-                self.coefficients.ZT * z_scaled * temperature
-                + self.coefficients.T * temperature
-                + self.coefficients.Z * z_scaled
-                + self.coefficients.c
-            )
-            * scale
-        )
+        return z_to_iwc(self.coefficients, z_scaled, temperature) * scale
 
     def _get_z_factor(self) -> float:
         """Returns empirical scaling factor for radar echo."""