Index: trunk/LMDZ.MARS/util/display_netcdf.py
===================================================================
--- trunk/LMDZ.MARS/util/display_netcdf.py	(revision 3798)
+++ trunk/LMDZ.MARS/util/display_netcdf.py	(revision 3808)
@@ -11,6 +11,6 @@
   - 2D latitude/longitude map
   - 2D (Time × another dimension)
-  - Variables with dimension “physical_points” as 2D map if lat/lon present,
-    or generic 2D plot if the remaining axes are spatial
+  - Variables with dimension “physical_points” displayed on a 2D map if lat/lon are present
+  - Optionally average over latitude and plot longitude vs. time heatmap
   - Scalar output (ndim == 0 after slicing)
 
@@ -18,19 +18,20 @@
   1) Command-line mode:
        python display_netcdf.py /path/to/your_file.nc --variable VAR_NAME \
-           [--time-index 0] [--cmap viridis] [--output out.png] \
-           [--extra-indices '{"dim1": idx1, "dim2": idx2}']
+           [--time-index 0] [--alt-index 0] [--cmap viridis] [--avg-lat] \
+           [--output out.png] [--extra-indices '{"nslope": 1}']
 
     --variable     : Name of the variable to visualize.
-    --time-index   : Index along the Time dimension (ignored for purely 1D time series).
-    --alt-index    : Index along the altitude dimension, if present.
-    --cmap         : Matplotlib colormap for contourf (default: "jet").
+    --time-index   : Index along the Time dimension (0-based, ignored for purely 1D time series).
+    --alt-index    : Index along the altitude dimension (0-based), if present.
+    --cmap         : Matplotlib colormap (default: "jet").
+    --avg-lat      : Average over latitude and plot longitude vs. time heatmap.
     --output       : If provided, save the figure to this filename instead of displaying.
-    --extra-indices: JSON string to fix indices for any dimensions other than Time, lat, lon, or altitude.
-                     Example: '{"nslope": 0, "physical_points": 2}'
-                     Omitting a dimension means it remains unfixed (useful to plot a 1D profile).
+    --extra-indices: JSON string to fix indices for any other dimensions.
+                     For any dimension whose name contains "slope", use 1-based numbering here.
+                     Example: '{"nslope": 1, "physical_points": 3}'
 
   2) Interactive mode:
        python display_netcdf.py
-       (The script will prompt for the NetCDF file, the variable, etc.)
+       (The script will prompt for everything, including averaging option.)
 """
 
@@ -43,4 +44,5 @@
 import numpy as np
 import matplotlib.pyplot as plt
+import matplotlib.tri as mtri
 from netCDF4 import Dataset
 
@@ -105,561 +107,492 @@
 
 
-def plot_variable(dataset, varname, time_index=None, alt_index=None, colormap="jet",
-                  output_path=None, extra_indices=None):
-    """
-    Extracts the requested slice from the variable and plots it according to the data shape:
-
-    - Pure 1D time series → time-series plot
-    - After slicing:
-        • If data_slice.ndim == 0 → print the scalar value
-        • If data_slice.ndim == 1:
-            • If the remaining dimension is “subsurface_layers” (or another known coordinate) → vertical profile
-            • Else → simple plot vs. index
-        • If data_slice.ndim == 2:
-            • If lat/lon exist → contourf map
-            • Else → imshow generic 2D plot
-    - If data_slice.ndim is neither 0, 1, nor 2 → error message
-
-    Parameters
-    ----------
-    dataset       : netCDF4.Dataset object (already open)
-    varname       : name of the variable to plot
-    time_index    : int or None (if variable has a time dimension, ignored for pure time series)
-    alt_index     : int or None (if variable has an altitude dimension)
-    colormap      : string colormap name (passed to plt.contourf)
-    output_path   : string filepath to save figure, or None to display interactively
-    extra_indices : dict { dimension_name (str) : index (int) } for slicing all
-                     dimensions except Time/lat/lon/alt. If a dimension is not
-                     included, it remains “slice(None)” (useful for 1D plots).
+def plot_variable(dataset, varname, time_index=None, alt_index=None,
+                  colormap="jet", output_path=None, extra_indices=None,
+                  avg_lat=False):
+    """
+    Core plotting logic: reads the variable, handles masks,
+    determines dimensionality, and creates the appropriate plot:
+      - 1D time series
+      - 1D profiles or physical_points maps
+      - 2D lat×lon or generic 2D
+      - Time×lon heatmap if avg_lat=True
+      - Scalar printing
     """
     var = dataset.variables[varname]
-    dims = var.dimensions  # tuple of dimension names
-
-    # Read the full data (could be a masked array)
+    dims = var.dimensions
+
+    # Read full data
     try:
         data_full = var[:]
     except Exception as e:
-        print(f"Error: Cannot read data for variable '{varname}': {e}")
-        return
-
-    # Convert masked array to NaN
+        print(f"Error: Cannot read data for '{varname}': {e}")
+        return
     if hasattr(data_full, "mask"):
         data_full = np.where(data_full.mask, np.nan, data_full.data)
 
-    # ------------------------------------------------------------------------
-    # 1) Pure 1D time series (dims == ('Time',) or equivalent)
+    # Pure 1D time series
     if len(dims) == 1 and find_dim_index(dims, TIME_DIMS) is not None:
-        # Plot the time series directly
-        time_varname = find_coord_var(dataset, TIME_DIMS)
-        if time_varname:
-            time_vals = dataset.variables[time_varname][:]
-            if hasattr(time_vals, "mask"):
-                time_vals = np.where(time_vals.mask, np.nan, time_vals.data)
-        else:
-            time_vals = np.arange(data_full.shape[0])
-
+        time_var = find_coord_var(dataset, TIME_DIMS)
+        tvals = (dataset.variables[time_var][:] if time_var
+                 else np.arange(data_full.shape[0]))
+        if hasattr(tvals, "mask"):
+            tvals = np.where(tvals.mask, np.nan, tvals.data)
         plt.figure()
-        plt.plot(time_vals, data_full, marker='o')
-        xlabel = time_varname if time_varname else "Time Index"
-        plt.xlabel(xlabel)
-        ylabel = varname
-        if hasattr(var, "units"):
-            ylabel += f" ({var.units})"
-        plt.ylabel(ylabel)
-        plt.title(f"{varname} vs {xlabel}")
-
+        plt.plot(tvals, data_full, marker="o")
+        plt.xlabel(time_var or "Time Index")
+        plt.ylabel(varname + (f" ({var.units})" if hasattr(var, "units") else ""))
+        plt.title(f"{varname} vs {time_var or 'Index'}")
         if output_path:
-            try:
-                plt.savefig(output_path, bbox_inches="tight")
-                print(f"Figure saved to '{output_path}'")
-            except Exception as e:
-                print(f"Error saving figure: {e}")
+            plt.savefig(output_path, bbox_inches="tight")
+            print(f"Saved to {output_path}")
         else:
             plt.show()
-        plt.close()
-        return
-    # ------------------------------------------------------------------------
-
-    # Identify special dimension indices
+        return
+
+    # Identify dims
     t_idx = find_dim_index(dims, TIME_DIMS)
-    a_idx = find_dim_index(dims, ALT_DIMS)
     lat_idx = find_dim_index(dims, LAT_DIMS)
     lon_idx = find_dim_index(dims, LON_DIMS)
-
-    # Build the slicer list
+    a_idx = find_dim_index(dims, ALT_DIMS)
+
+    # Average over latitude & plot time × lon heatmap
+    if avg_lat and t_idx is not None and lat_idx is not None and lon_idx is not None:
+        # mean over lat axis
+        data_avg = np.nanmean(data_full, axis=lat_idx)
+        # data_avg shape: (time, lon, ...)
+        # we assume no other unfixed dims
+        # get coordinates
+        time_var = find_coord_var(dataset, TIME_DIMS)
+        lon_var = find_coord_var(dataset, LON_DIMS)
+        tvals = dataset.variables[time_var][:]
+        lons = dataset.variables[lon_var][:]
+        if hasattr(tvals, "mask"):
+            tvals = np.where(tvals.mask, np.nan, tvals.data)
+        if hasattr(lons, "mask"):
+            lons = np.where(lons.mask, np.nan, lons.data)
+        plt.figure(figsize=(10, 6))
+        plt.pcolormesh(lons, tvals, data_avg, shading="auto", cmap=colormap)
+        plt.xlabel(f"Longitude ({getattr(dataset.variables[lon_var], 'units', 'deg')})")
+        plt.ylabel(time_var)
+        cbar = plt.colorbar()
+        cbar.set_label(varname + (f" ({var.units})" if hasattr(var, "units") else ""))
+        plt.title(f"{varname} averaged over latitude")
+        if output_path:
+            plt.savefig(output_path, bbox_inches="tight")
+            print(f"Saved to {output_path}")
+        else:
+            plt.show()
+        return
+
+    # Build slicer for other cases
     slicer = [slice(None)] * len(dims)
-
-    # Apply slicing on Time and altitude if specified
     if t_idx is not None:
         if time_index is None:
-            print("Error: Variable has a time dimension; please supply a time index.")
+            print("Error: please supply a time index.")
             return
         slicer[t_idx] = time_index
     if a_idx is not None:
         if alt_index is None:
-            print("Error: Variable has an altitude dimension; please supply an altitude index.")
+            print("Error: please supply an altitude index.")
             return
         slicer[a_idx] = alt_index
 
-    # Apply slicing on all “extra” dimensions (except Time/lat/lon/alt)
     if extra_indices is None:
         extra_indices = {}
-    for dim_name, idx_val in extra_indices.items():
-        if dim_name in dims:
-            dim_index = dims.index(dim_name)
-            slicer[dim_index] = idx_val
-
-    # Extract the sliced data
+    for dn, idx_val in extra_indices.items():
+        if dn in dims:
+            slicer[dims.index(dn)] = idx_val
+
+    # Extract slice
     try:
-        data_slice = data_full[tuple(slicer)]
+        dslice = data_full[tuple(slicer)]
     except Exception as e:
-        print(f"Error: Could not slice variable '{varname}': {e}")
-        return
-
-    # CASE: After slicing, if data_slice.ndim == 0 → scalar
-    if np.ndim(data_slice) == 0:
-        try:
-            scalar_val = float(data_slice)
-        except Exception:
-            scalar_val = data_slice
-        print(f"Scalar result for '{varname}': {scalar_val}")
-        return
-
-    # CASE: After slicing, if data_slice.ndim == 1 (vertical profile or simple vector)
-    if data_slice.ndim == 1:
-        # Identify the remaining dimension
-        rem_dim = None
-        for di, dname in enumerate(dims):
-            if slicer[di] == slice(None):
-                rem_dim = (di, dname)
-                break
-
-        if rem_dim is not None:
-            di, dname = rem_dim
-            coord_var = None
-
-            # If it's "subsurface_layers", look for coordinate "soildepth"
+        print(f"Error slicing '{varname}': {e}")
+        return
+
+    # Scalar
+    if np.ndim(dslice) == 0:
+        print(f"Scalar '{varname}': {float(dslice)}")
+        return
+
+    # 1D: vector, profile, or physical_points
+    if dslice.ndim == 1:
+        rem = [(i, name) for i, name in enumerate(dims) if slicer[i] == slice(None)]
+        if rem:
+            di, dname = rem[0]
+            # physical_points → interpolated map
+            if dname.lower() == "physical_points":
+                latv = find_coord_var(dataset, LAT_DIMS)
+                lonv = find_coord_var(dataset, LON_DIMS)
+                if latv and lonv:
+                    lats = dataset.variables[latv][:]
+                    lons = dataset.variables[lonv][:]
+                    if hasattr(lats, "mask"):
+                        lats = np.where(lats.mask, np.nan, lats.data)
+                    if hasattr(lons, "mask"):
+                        lons = np.where(lons.mask, np.nan, lons.data)
+                    triang = mtri.Triangulation(lons, lats)
+                    plt.figure(figsize=(8, 6))
+                    cf = plt.tricontourf(triang, dslice, cmap=colormap)
+                    cbar = plt.colorbar(cf)
+                    cbar.set_label(varname + (f" ({var.units})" if hasattr(var, "units") else ""))
+                    plt.xlabel(f"Longitude ({getattr(dataset.variables[lonv], 'units', 'deg')})")
+                    plt.ylabel(f"Latitude ({getattr(dataset.variables[latv], 'units', 'deg')})")
+                    plt.title(f"{varname} (interpolated map over physical_points)")
+                    if output_path:
+                        plt.savefig(output_path, bbox_inches="tight")
+                        print(f"Saved to {output_path}")
+                    else:
+                        plt.show()
+                    return
+            # vertical profile?
+            coord = None
             if dname.lower() == "subsurface_layers" and "soildepth" in dataset.variables:
-                coord_var = "soildepth"
-            # If there is a variable with the same name, use it
+                coord = "soildepth"
             elif dname in dataset.variables:
-                coord_var = dname
-
-            if coord_var:
-                coord_vals = dataset.variables[coord_var][:]
-                if hasattr(coord_vals, "mask"):
-                    coord_vals = np.where(coord_vals.mask, np.nan, coord_vals.data)
-                x = data_slice
-                y = coord_vals
-
+                coord = dname
+            if coord:
+                coords = dataset.variables[coord][:]
+                if hasattr(coords, "mask"):
+                    coords = np.where(coords.mask, np.nan, coords.data)
                 plt.figure()
-                plt.plot(x, y, marker='o')
-                # Invert Y-axis if it's a depth coordinate
+                plt.plot(dslice, coords, marker="o")
                 if dname.lower() == "subsurface_layers":
                     plt.gca().invert_yaxis()
-
-                xlabel = varname
-                if hasattr(var, "units"):
-                    xlabel += f" ({var.units})"
-                plt.xlabel(xlabel)
-
-                ylabel = coord_var
-                if hasattr(dataset.variables[coord_var], "units"):
-                    ylabel += f" ({dataset.variables[coord_var].units})"
-                plt.ylabel(ylabel)
-
-                plt.title(f"{varname} vs {coord_var}")
-
+                plt.xlabel(varname + (f" ({var.units})" if hasattr(var, "units") else ""))
+                plt.ylabel(coord + (f" ({dataset.variables[coord].units})" if hasattr(dataset.variables[coord], "units") else ""))
+                plt.title(f"{varname} vs {coord}")
                 if output_path:
-                    try:
-                        plt.savefig(output_path, bbox_inches="tight")
-                        print(f"Figure saved to '{output_path}'")
-                    except Exception as e:
-                        print(f"Error saving figure: {e}")
+                    plt.savefig(output_path, bbox_inches="tight")
+                    print(f"Saved to {output_path}")
                 else:
                     plt.show()
-                plt.close()
                 return
-            else:
-                # No known coordinate found → simple plot vs index
-                plt.figure()
-                plt.plot(data_slice, marker='o')
-                plt.xlabel("Index")
-                ylabel = varname
-                if hasattr(var, "units"):
-                    ylabel += f" ({var.units})"
-                plt.ylabel(ylabel)
-                plt.title(f"{varname} (1D)")
-
-                if output_path:
-                    try:
-                        plt.savefig(output_path, bbox_inches="tight")
-                        print(f"Figure saved to '{output_path}'")
-                    except Exception as e:
-                        print(f"Error saving figure: {e}")
-                else:
-                    plt.show()
-                plt.close()
-                return
-
-        else:
-            # Unable to identify the remaining dimension → error
-            print(f"Error: After slicing, data for '{varname}' is 1D but remaining dimension is unknown.")
-            return
-
-    # CASE: After slicing, if data_slice.ndim == 2
-    if data_slice.ndim == 2:
-        # If lat and lon exist in the original dims, re-find their indices
+        # generic 1D
+        plt.figure()
+        plt.plot(dslice, marker="o")
+        plt.xlabel("Index")
+        plt.ylabel(varname + (f" ({var.units})" if hasattr(var, "units") else ""))
+        plt.title(f"{varname} (1D)")
+        if output_path:
+            plt.savefig(output_path, bbox_inches="tight")
+            print(f"Saved to {output_path}")
+        else:
+            plt.show()
+        return
+
+    # 2D: map or generic
+    if dslice.ndim == 2:
         lat_idx2 = find_dim_index(dims, LAT_DIMS)
         lon_idx2 = find_dim_index(dims, LON_DIMS)
-
         if lat_idx2 is not None and lon_idx2 is not None:
-            # We have a 2D variable on a lat×lon grid
-            lat_varname = find_coord_var(dataset, LAT_DIMS)
-            lon_varname = find_coord_var(dataset, LON_DIMS)
-            if lat_varname is None or lon_varname is None:
-                print("Error: Could not locate latitude/longitude variables in the dataset.")
-                return
-
-            lat_var = dataset.variables[lat_varname][:]
-            lon_var = dataset.variables[lon_varname][:]
-            if hasattr(lat_var, "mask"):
-                lat_var = np.where(lat_var.mask, np.nan, lat_var.data)
-            if hasattr(lon_var, "mask"):
-                lon_var = np.where(lon_var.mask, np.nan, lon_var.data)
-
-            # Build 2D coordinate arrays
-            if lat_var.ndim == 1 and lon_var.ndim == 1:
-                lon2d, lat2d = np.meshgrid(lon_var, lat_var)
-            elif lat_var.ndim == 2 and lon_var.ndim == 2:
-                lat2d, lon2d = lat_var, lon_var
+            latv = find_coord_var(dataset, LAT_DIMS)
+            lonv = find_coord_var(dataset, LON_DIMS)
+            lats = dataset.variables[latv][:]
+            lons = dataset.variables[lonv][:]
+            if hasattr(lats, "mask"):
+                lats = np.where(lats.mask, np.nan, lats.data)
+            if hasattr(lons, "mask"):
+                lons = np.where(lons.mask, np.nan, lons.data)
+            if lats.ndim == 1 and lons.ndim == 1:
+                lon2d, lat2d = np.meshgrid(lons, lats)
             else:
-                print("Error: Latitude and longitude must both be either 1D or 2D.")
-                return
-
+                lat2d, lon2d = lats, lons
             plt.figure(figsize=(10, 6))
-            cf = plt.contourf(lon2d, lat2d, data_slice, cmap=colormap)
+            cf = plt.contourf(lon2d, lat2d, dslice, cmap=colormap)
             cbar = plt.colorbar(cf)
-            if hasattr(var, "units"):
-                cbar.set_label(f"{varname} ({var.units})")
-            else:
-                cbar.set_label(varname)
-
-            lon_label = f"Longitude ({getattr(dataset.variables[lon_varname], 'units', 'degrees')})"
-            lat_label = f"Latitude ({getattr(dataset.variables[lat_varname], 'units', 'degrees')})"
-            plt.xlabel(lon_label)
-            plt.ylabel(lat_label)
+            cbar.set_label(varname + (f" ({var.units})" if hasattr(var, "units") else ""))
+            plt.xlabel(f"Longitude ({getattr(dataset.variables[lonv], 'units', 'deg')})")
+            plt.ylabel(f"Latitude ({getattr(dataset.variables[latv], 'units', 'deg')})")
             plt.title(f"{varname} (lat × lon)")
-
             if output_path:
-                try:
-                    plt.savefig(output_path, bbox_inches="tight")
-                    print(f"Figure saved to '{output_path}'")
-                except Exception as e:
-                    print(f"Error saving figure: {e}")
+                plt.savefig(output_path, bbox_inches="tight")
+                print(f"Saved to {output_path}")
             else:
                 plt.show()
-            plt.close()
             return
-
-        else:
-            # No lat/lon → two non-geographical dimensions; plot with imshow
-            plt.figure(figsize=(8, 6))
-            plt.imshow(data_slice, aspect='auto')
-            cb_label = varname
-            if hasattr(var, "units"):
-                cb_label += f" ({var.units})"
-            plt.colorbar(label=cb_label)
-            plt.xlabel("Dimension 2 Index")
-            plt.ylabel("Dimension 1 Index")
-            plt.title(f"{varname} (2D without lat/lon)")
-
-            if output_path:
-                try:
-                    plt.savefig(output_path, bbox_inches="tight")
-                    print(f"Figure saved to '{output_path}'")
-                except Exception as e:
-                    print(f"Error saving figure: {e}")
-            else:
-                plt.show()
-            plt.close()
-            return
-
-    # CASE: data_slice.ndim is neither 0, 1, nor 2
-    print(f"Error: After slicing, data for '{varname}' has ndim={data_slice.ndim}, which is not supported.")
-    return
+        # generic 2D
+        plt.figure(figsize=(8, 6))
+        plt.imshow(dslice, aspect="auto")
+        plt.colorbar(label=varname + (f" ({var.units})" if hasattr(var, "units") else ""))
+        plt.xlabel("Dim 2 index")
+        plt.ylabel("Dim 1 index")
+        plt.title(f"{varname} (2D)")
+        if output_path:
+            plt.savefig(output_path, bbox_inches="tight")
+            print(f"Saved to {output_path}")
+        else:
+            plt.show()
+        return
+
+    print(f"Error: ndim={dslice.ndim} not supported.")
 
 
 def visualize_variable_interactive(nc_path=None):
     """
-    Interactive mode: prompts for the NetCDF file if not provided, then for the variable,
-    then for Time/altitude indices (skipped entirely if variable is purely 1D over Time),
-    and for each other dimension offers to fix an index or to plot along that dimension (by typing 'f').
-
-    If a dimension has length 1, the index 0 is chosen automatically.
-    """
-    # Determine file path
+    Interactive mode: prompts for file, variable, displays dims,
+    handles special case of pure time series, then guides user
+    through any needed index selections.
+    """
+    # File selection
     if nc_path:
-        file_input = nc_path
+        path = nc_path
     else:
         readline.set_completer(complete_filename)
         readline.parse_and_bind("tab: complete")
-        file_input = input("Enter the path to the NetCDF file: ").strip()
-
-    if not file_input:
-        print("No file specified. Exiting.")
-        return
-    if not os.path.isfile(file_input):
-        print(f"Error: '{file_input}' not found.")
-        return
-
-    try:
-        ds = Dataset(file_input, mode="r")
-    except Exception as e:
-        print(f"Error: Unable to open '{file_input}': {e}")
-        return
-
-    varnames = list(ds.variables.keys())
-    if not varnames:
-        print("Error: No variables found in the dataset.")
+        path = input("Enter path to NetCDF file: ").strip()
+    if not os.path.isfile(path):
+        print(f"Error: '{path}' not found."); return
+    ds = Dataset(path, "r")
+
+    # Variable selection with autocomplete
+    vars_ = list(ds.variables.keys())
+    if not vars_:
+        print("No variables found."); ds.close(); return
+    if len(vars_) == 1:
+        var = vars_[0]; print(f"Selected '{var}'")
+    else:
+        print("Available variables:")
+        for v in vars_:
+            print(f"  - {v}")
+        readline.set_completer(make_varname_completer(vars_))
+        readline.parse_and_bind("tab: complete")
+        var = input("Variable name: ").strip()
+        if var not in ds.variables:
+            print("Unknown variable."); ds.close(); return
+
+    # DISPLAY DIMENSIONS AND SIZES
+    dims  = ds.variables[var].dimensions
+    shape = ds.variables[var].shape
+    print(f"\nVariable '{var}' has {len(dims)} dimensions:")
+    for name, size in zip(dims, shape):
+        print(f"  - {name}: size {size}")
+    print()
+
+    # Identify dimension indices
+    t_idx   = find_dim_index(dims, TIME_DIMS)
+    lat_idx = find_dim_index(dims, LAT_DIMS)
+    lon_idx = find_dim_index(dims, LON_DIMS)
+    a_idx   = find_dim_index(dims, ALT_DIMS)
+
+    # SPECIAL CASE: time-only series (all others singleton) → plot directly
+    if (
+        t_idx is not None and shape[t_idx] > 1 and
+        all(shape[i] == 1 for i in range(len(dims)) if i != t_idx)
+    ):
+        print("Detected single-point spatial dims; plotting time series…")
+        # récupérer les valeurs
+        var_obj = ds.variables[var]
+        data = var_obj[:].squeeze()   # shape (time,)
+        # temps
+        time_var = find_coord_var(ds, TIME_DIMS)
+        if time_var:
+            tvals = ds.variables[time_var][:]
+        else:
+            tvals = np.arange(data.shape[0])
+        # masque éventuel
+        if hasattr(data, "mask"):
+            data = np.where(data.mask, np.nan, data.data)
+        if hasattr(tvals, "mask"):
+            tvals = np.where(tvals.mask, np.nan, tvals.data)
+        # tracé
+        plt.figure()
+        plt.plot(tvals, data, marker="o")
+        plt.xlabel(time_var or "Time Index")
+        plt.ylabel(var + (f" ({var_obj.units})" if hasattr(var_obj, "units") else ""))
+        plt.title(f"{var} vs {time_var or 'Index'}")
+        plt.show()
         ds.close()
         return
 
-    # Auto-select if only one variable
-    if len(varnames) == 1:
-        var_input = varnames[0]
-        print(f"Automatically selected the only variable: '{var_input}'")
-    else:
-        print("\nAvailable variables:")
-        for name in varnames:
-            print(f"  - {name}")
-        print()
-        readline.set_completer(make_varname_completer(varnames))
-        var_input = input("Enter the name of the variable to visualize: ").strip()
-        if var_input not in ds.variables:
-            print(f"Error: Variable '{var_input}' not found. Exiting.")
-            ds.close()
-            return
-
-    dims = ds.variables[var_input].dimensions  # tuple of dimension names
-
-    # If the variable is purely 1D over Time, plot immediately without asking for time index
-    if len(dims) == 1 and find_dim_index(dims, TIME_DIMS) is not None:
-        plot_variable(
-            dataset=ds,
-            varname=var_input,
-            time_index=None,
-            alt_index=None,
-            colormap="jet",
-            output_path=None,
-            extra_indices=None
-        )
+    # Ask average over latitude only if Time, lat AND lon each >1
+    avg_lat = False
+    if (
+        t_idx   is not None and shape[t_idx]  > 1 and
+        lat_idx is not None and shape[lat_idx] > 1 and
+        lon_idx is not None and shape[lon_idx] > 1
+    ):
+        u = input("Average over latitude & plot lon vs time? [y/n]: ").strip().lower()
+        avg_lat = (u == "y")
+
+    # Time index prompt
+    ti = None
+    if t_idx is not None:
+        L = shape[t_idx]
+        if L > 1:
+            while True:
+                u = input(f"Enter time index [0..{L-1}]: ").strip()
+                try:
+                    ti = int(u)
+                    if 0 <= ti < L:
+                        break
+                except:
+                    pass
+                print("Invalid.")
+        else:
+            ti = 0; print("Only one time; using 0.")
+
+    # Altitude index prompt
+    ai = None
+    if a_idx is not None:
+        L = shape[a_idx]
+        if L > 1:
+            while True:
+                u = input(f"Enter altitude index [0..{L-1}]: ").strip()
+                try:
+                    ai = int(u)
+                    if 0 <= ai < L:
+                        break
+                except:
+                    pass
+                print("Invalid.")
+        else:
+            ai = 0; print("Only one altitude; using 0.")
+
+    # Other dims
+    extra = {}
+    for idx, dname in enumerate(dims):
+        if idx in (t_idx, a_idx):
+            continue
+        if dname.lower() in LAT_DIMS + LON_DIMS and shape[idx] == 1:
+            extra[dname] = 0
+            continue
+        L = shape[idx]
+        if L == 1:
+            extra[dname] = 0
+            continue
+        if "slope" in dname.lower():
+            prompt = f"Enter slope number [1..{L}] for '{dname}': "
+        else:
+            prompt = f"Enter index [0..{L-1}] or 'f' to plot '{dname}': "
+        while True:
+            u = input(prompt).strip().lower()
+            if u == "f" and "slope" not in dname.lower():
+                break
+            try:
+                iv = int(u)
+                if "slope" in dname.lower():
+                    if 1 <= iv <= L:
+                        extra[dname] = iv - 1
+                        break
+                else:
+                    if 0 <= iv < L:
+                        extra[dname] = iv
+                        break
+            except:
+                pass
+            print("Invalid.")
+
+    plot_variable(ds, var, time_index=ti, alt_index=ai,
+                  colormap="jet", output_path=None,
+                  extra_indices=extra, avg_lat=avg_lat)
+    ds.close()
+
+
+def visualize_variable_cli(nc_file, varname, time_index, alt_index,
+                           colormap, output_path, extra_json, avg_lat):
+    """
+    Command-line mode: visualize directly, parsing the --extra-indices argument (JSON string).
+    """
+    if not os.path.isfile(nc_file):
+        print(f"Error: '{nc_file}' not found."); return
+    ds = Dataset(nc_file, "r")
+    if varname not in ds.variables:
+        print(f"Variable '{varname}' not in file."); ds.close(); return
+
+    # DISPLAY DIMENSIONS AND SIZES
+    dims  = ds.variables[varname].dimensions
+    shape = ds.variables[varname].shape
+    print(f"\nVariable '{varname}' has {len(dims)} dimensions:")
+    for name, size in zip(dims, shape):
+        print(f"  - {name}: size {size}")
+    print()
+
+    # SPECIAL CASE: time-only → plot directly
+    t_idx = find_dim_index(dims, TIME_DIMS)
+    if (
+        t_idx is not None and shape[t_idx] > 1 and
+        all(shape[i] == 1 for i in range(len(dims)) if i != t_idx)
+    ):
+        print("Detected single-point spatial dims; plotting time series…")
+        # même logique que ci‑dessus
+        var_obj = ds.variables[varname]
+        data = var_obj[:].squeeze()
+        time_var = find_coord_var(ds, TIME_DIMS)
+        if time_var:
+            tvals = ds.variables[time_var][:]
+        else:
+            tvals = np.arange(data.shape[0])
+        if hasattr(data, "mask"):
+            data = np.where(data.mask, np.nan, data.data)
+        if hasattr(tvals, "mask"):
+            tvals = np.where(tvals.mask, np.nan, tvals.data)
+        plt.figure()
+        plt.plot(tvals, data, marker="o")
+        plt.xlabel(time_var or "Time Index")
+        plt.ylabel(varname + (f" ({var_obj.units})" if hasattr(var_obj, "units") else ""))
+        plt.title(f"{varname} vs {time_var or 'Index'}")
+        if output_path:
+            plt.savefig(output_path, bbox_inches="tight")
+            print(f"Saved to {output_path}")
+        else:
+            plt.show()
         ds.close()
         return
 
-    # Otherwise, proceed to prompt for time/altitude and other dimensions
-    time_idx = None
-    alt_idx = None
-
-    # Prompt for time index if applicable
-    t_idx = find_dim_index(dims, TIME_DIMS)
-    if t_idx is not None:
-        time_len = ds.variables[var_input].shape[t_idx]
-        if time_len > 1:
-            while True:
-                try:
-                    user_t = input(f"Enter time index [0..{time_len - 1}]: ").strip()
-                    if user_t == "":
-                        print("No time index entered. Exiting.")
-                        ds.close()
-                        return
-                    ti = int(user_t)
-                    if 0 <= ti < time_len:
-                        time_idx = ti
-                        break
-                except ValueError:
-                    pass
-                print(f"Invalid index. Enter an integer between 0 and {time_len - 1}.")
-        else:
-            time_idx = 0
-            print("Only one time step available; using index 0.")
-
-    # Prompt for altitude index if applicable
-    a_idx = find_dim_index(dims, ALT_DIMS)
-    if a_idx is not None:
-        alt_len = ds.variables[var_input].shape[a_idx]
-        if alt_len > 1:
-            while True:
-                try:
-                    user_a = input(f"Enter altitude index [0..{alt_len - 1}]: ").strip()
-                    if user_a == "":
-                        print("No altitude index entered. Exiting.")
-                        ds.close()
-                        return
-                    ai = int(user_a)
-                    if 0 <= ai < alt_len:
-                        alt_idx = ai
-                        break
-                except ValueError:
-                    pass
-                print(f"Invalid index. Enter an integer between 0 and {alt_len - 1}.")
-        else:
-            alt_idx = 0
-            print("Only one altitude level available; using index 0.")
-
-    # Identify other dimensions (excluding Time/lat/lon/alt)
-    other_dims = []
-    for idx_dim, dim_name in enumerate(dims):
-        if idx_dim == t_idx or idx_dim == a_idx:
-            continue
-        if dim_name.lower() in (d.lower() for d in LAT_DIMS + LON_DIMS):
-            continue
-        other_dims.append((idx_dim, dim_name))
-
-    # For each other dimension, ask user to fix an index or type 'f' to plot along that dimension
-    extra_indices = {}
-    for idx_dim, dim_name in other_dims:
-        dim_len = ds.variables[var_input].shape[idx_dim]
-        if dim_len == 1:
-            extra_indices[dim_name] = 0
-            print(f"Dimension '{dim_name}' has length 1; using index 0.")
-        else:
-            while True:
-                prompt = (
-                    f"Enter index for '{dim_name}' [0..{dim_len - 1}] "
-                    f"or 'f' to plot along '{dim_name}': "
-                )
-                user_i = input(prompt).strip().lower()
-                if user_i == 'f':
-                    # Leave this dimension unfixed → no key in extra_indices
-                    break
-                if user_i == "":
-                    print("No index entered. Exiting.")
-                    ds.close()
-                    return
-                try:
-                    idx_val = int(user_i)
-                    if 0 <= idx_val < dim_len:
-                        extra_indices[dim_name] = idx_val
-                        break
-                except ValueError:
-                    pass
-                print(f"Invalid index. Enter an integer between 0 and {dim_len - 1}, or 'f'.")
-
-    # Finally, call plot_variable with collected indices
-    plot_variable(
-        dataset=ds,
-        varname=var_input,
-        time_index=time_idx,
-        alt_index=alt_idx,
-        colormap="jet",
-        output_path=None,
-        extra_indices=extra_indices
-    )
-    ds.close()
-
-
-def visualize_variable_cli(nc_path, varname, time_index, alt_index, colormap, output_path, extra_json):
-    """
-    Command-line mode: visualize directly, parsing the --extra-indices argument (JSON string).
-    """
-    if not os.path.isfile(nc_path):
-        print(f"Error: '{nc_path}' not found.")
-        return
-
-    try:
-        ds = Dataset(nc_path, mode="r")
-    except Exception as e:
-        print(f"Error: Unable to open '{nc_path}': {e}")
-        return
-
-    if varname not in ds.variables:
-        print(f"Error: Variable '{varname}' not found in '{nc_path}'.")
-        ds.close()
-        return
-
-    # Parse extra_indices if provided
-    extra_indices = {}
+    # Si --avg-lat mais lat/lon/Time non compatibles → désactive
+    lat_idx = find_dim_index(dims, LAT_DIMS)
+    lon_idx = find_dim_index(dims, LON_DIMS)
+    if avg_lat and not (
+        t_idx   is not None and shape[t_idx]  > 1 and
+        lat_idx is not None and shape[lat_idx] > 1 and
+        lon_idx is not None and shape[lon_idx] > 1
+    ):
+        print("Note: disabling --avg-lat (requires Time, lat & lon each >1).")
+        avg_lat = False
+
+    # Parse extra indices JSON
+    extra = {}
     if extra_json:
         try:
             parsed = json.loads(extra_json)
-            if isinstance(parsed, dict):
-                for k, v in parsed.items():
-                    if isinstance(k, str) and isinstance(v, int):
-                        extra_indices[k] = v
-            else:
-                print("Warning: --extra-indices is not a JSON object. Ignored.")
-        except json.JSONDecodeError:
-            print("Warning: --extra-indices is not valid JSON. Ignored.")
-
-    plot_variable(
-        dataset=ds,
-        varname=varname,
-        time_index=time_index,
-        alt_index=alt_index,
-        colormap=colormap,
-        output_path=output_path,
-        extra_indices=extra_indices
-    )
+            for k, v in parsed.items():
+                if isinstance(v, int):
+                    if "slope" in k.lower():
+                        extra[k] = v - 1
+                    else:
+                        extra[k] = v
+        except:
+            print("Warning: bad extra-indices.")
+
+    plot_variable(ds, varname, time_index, alt_index,
+                  colormap, output_path, extra, avg_lat)
     ds.close()
 
 
 def main():
-    parser = argparse.ArgumentParser(
-        description="Visualize a 1D/2D slice of a NetCDF variable on various dimension types."
-    )
-    parser.add_argument(
-        "nc_file",
-        nargs="?",
-        help="Path to the NetCDF file (interactive if omitted)."
-    )
-    parser.add_argument(
-        "--variable", "-v",
-        help="Name of the variable to visualize."
-    )
-    parser.add_argument(
-        "--time-index", "-t",
-        type=int,
-        help="Index on the Time dimension, if applicable (ignored for pure 1D time series)."
-    )
-    parser.add_argument(
-        "--alt-index", "-a",
-        type=int,
-        help="Index on the altitude dimension, if applicable."
-    )
-    parser.add_argument(
-        "--cmap", "-c",
-        default="jet",
-        help="Matplotlib colormap (default: 'jet')."
-    )
-    parser.add_argument(
-        "--output", "-o",
-        help="If provided, save the figure to this file instead of displaying it."
-    )
-    parser.add_argument(
-        "--extra-indices", "-e",
-        help="JSON string to fix indices of dimensions outside Time/lat/lon/alt. "
-             "Example: '{\"nslope\":0, \"physical_points\":2}'."
-    )
-
+    parser = argparse.ArgumentParser()
+    parser.add_argument("nc_file", nargs="?", help="NetCDF file (omit for interactive)")
+    parser.add_argument("-v", "--variable", help="Variable name")
+    parser.add_argument("-t", "--time-index", type=int, help="Time index (0-based)")
+    parser.add_argument("-a", "--alt-index", type=int, help="Altitude index (0-based)")
+    parser.add_argument("-c", "--cmap", default="jet", help="Colormap")
+    parser.add_argument("--avg-lat", action="store_true",
+                        help="Average over latitude (time × lon heatmap)")
+    parser.add_argument("-o", "--output", help="Save figure path")
+    parser.add_argument("-e", "--extra-indices", help="JSON for other dims")
     args = parser.parse_args()
 
-    # If both nc_file and variable are provided → CLI mode
     if args.nc_file and args.variable:
         visualize_variable_cli(
-            nc_path=args.nc_file,
-            varname=args.variable,
-            time_index=args.time_index,
-            alt_index=args.alt_index,
-            colormap=args.cmap,
-            output_path=args.output,
-            extra_json=args.extra_indices
+            args.nc_file, args.variable,
+            args.time_index, args.alt_index,
+            args.cmap, args.output,
+            args.extra_indices, args.avg_lat
         )
     else:
-        # Otherwise → fully interactive mode
-        visualize_variable_interactive(nc_path=args.nc_file)
+        visualize_variable_interactive(args.nc_file)
 
 
 if __name__ == "__main__":
     main()
-
