diff --git a/DESCRIPTION b/DESCRIPTION
index 0c5151978a47f47d02a0e1dfe2ea051c28b6d17a..d2444dab773613047218c7002202b443924f499b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ClimProjDiags
 Title: Set of Tools to Compute Various Climate Indices
-Version: 0.3.0
+Version: 0.3.1
 Authors@R: c(
     person("BSC-CNS", role = c("aut", "cph")),
     person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut"), comment = c(ORCID = "0000-0001-8568-3071")),
diff --git a/NEWS.md b/NEWS.md
index dd95442f5ba4fc20e1e46608107618ec7440c725..d72e8eba87f3961e8e032cfe8eb563b16944b2b4 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,6 @@
+# 0.3.1 (Release date: 2023-03-23)
+- Subset(): Prioritize the dimension names from names(dim(x)) rather than attribute 'dimensions'; If the input data doesn't have dimension names, the output doesn't have either.
+
 # 0.3.0 (Release date: 2023-02-28)
 - SelBox() and ShiftLon() to accept non-numerical data input  
 - SelBox() uses the latitude and longitude dimension name instead of index  
diff --git a/R/Subset.R b/R/Subset.R
index 37cb327e1fe7eee61061c49f6db140cf2d11d11b..d9cb384e2daf32cfb6aca0294f101db81596eff0 100644
--- a/R/Subset.R
+++ b/R/Subset.R
@@ -48,9 +48,11 @@ Subset <- function(x, along, indices, drop = FALSE) {
   }
 
   # Take the input array dimension names
+  x_has_names_dim <- TRUE
   dim_names <- names(dim(x))
   if (is.null(dim_names)) {
     dim_names <- attr(x, 'dimensions')
+    x_has_names_dim <- FALSE
   }
 
   if (!is.character(dim_names)) {
@@ -63,6 +65,8 @@ Subset <- function(x, along, indices, drop = FALSE) {
         warning("Found attribute 'dimensions' containing different dimension names from ",
                 "dim(names(x)). Use the latter one only.")
       }
+    } else if (is.null(names(dim(x))) & !is.null(attr(x, 'dimensions'))) {
+      names(dim(x)) <- dim_names
     }
   }
   # Check along
@@ -138,6 +142,10 @@ Subset <- function(x, along, indices, drop = FALSE) {
       }
     }
   }
+  # if names(dim(x)) doesn't exist, remove the dimension names
+  if (!x_has_names_dim) {
+    names(dim(subset)) <- NULL
+  }
 
   # Amend the final dimensions and put dimnames and attributes
   metadata <- attributes(x)
@@ -157,12 +165,14 @@ Subset <- function(x, along, indices, drop = FALSE) {
       }
     }
   } else if (is.character(dim_names) & !identical(dim_names, character(0))) {
-    names(metadata[['dim']]) <- dim_names
+    if (x_has_names_dim) {
+      names(metadata[['dim']]) <- dim_names
+    }
     if ('dimensions' %in% names(attributes(x))) {
       metadata[['dimensions']] <- dim_names
     }
   }
-
   attributes(subset) <- metadata
-  subset
+
+  return(subset)
 }
diff --git a/tests/testthat/test-Subset.R b/tests/testthat/test-Subset.R
index 202112d4bf28c1046c292e1dc91c92fcbf89e465..8144f96262ee9d76f32684565a9f1644a9979ba0 100644
--- a/tests/testthat/test-Subset.R
+++ b/tests/testthat/test-Subset.R
@@ -253,7 +253,17 @@ check.attributes = F
 )
 expect_equal(
 attributes(Subset(dat3, 'dat', 1, drop = FALSE)),
-list(dim = c(dat = 1, lat = 2, lon = 10), dimensions = c('dat', 'lat', 'lon'))
+list(dim = c(1, 2, 10), dimensions = c('dat', 'lat', 'lon'))
+)
+
+# drop: dat, lat
+expect_equal(
+attributes(Subset(dat3, c('dat', 'lat'), list(1, 1), drop = FALSE)),
+list(dim = c(1, 1, 10), dimensions = c('dat', 'lat', 'lon'))
+)
+expect_equal(
+attributes(Subset(dat3, c('dat', 'lat'), list(1, 1), drop = TRUE)),
+list(dim = c(10), dimensions = c('lon'))
 )
 
 # drop: lat, lon