diff --git a/R/CST_SplitDim.R b/R/CST_SplitDim.R
index bac378be1202675fd1fb08c080023e717d0d4184..4ab9e5a4dcf769bfb17848cd2d9220cc2cc0f24e 100644
--- a/R/CST_SplitDim.R
+++ b/R/CST_SplitDim.R
@@ -124,13 +124,16 @@ CST_SplitDim <- function(data, split_dim = 'time', indices = NULL,
       }
     }
   }
+  old_dim_name <- names(dim(data$data))
   # Call the function
   res <- SplitDim(data = data$data, split_dim = split_dim, 
                   indices = indices, freq = freq, 
                   new_dim_name = new_dim_name, 
                   dates = data$attrs$Dates, 
                   return_indices = return_indices)
-  if (inherits(res, 'list')) {
+
+  
+    if (inherits(res, 'list')) {
     data$data <- res$data
     # Split dim on Dates
     if (!is.null(res$dates)) {
@@ -139,10 +142,16 @@ CST_SplitDim <- function(data, split_dim = 'time', indices = NULL,
   } else {
     data$data <- res
   }
+  
   data$dims <- dim(data$data)
-
+  new_dim_name <- names(data$dims)[!names(data$dims) %in% (old_dim_name)]
+  
   # Coordinates
-  # TO DO: Subset splitted coordinate and add the new dimension coordinate.
+  data$coords[[split_dim]] <- 1:dim(data$data)[[split_dim]]
+  attr(data$coords[[split_dim]], 'indices') <- TRUE
+  data$coords[[new_dim_name]] <- 1:dim(data$data)[[new_dim_name]] 
+  attr(data$coords[[new_dim_name]], 'indices') <- TRUE
+  
   if (return_indices) {
     return(list(data = data, indices = res$indices))
   } else {
@@ -190,7 +199,7 @@ CST_SplitDim <- function(data, split_dim = 'time', indices = NULL,
 #'@importFrom ClimProjDiags Subset
 #'@export
 SplitDim <- function(data, split_dim = 'time', indices, freq = 'monthly',
-                     new_dim_name = NULL, dates = NULL, 
+                     new_dim_name = NULL, dates = NULL,
                      return_indices = FALSE) {
   # check data
   if (is.null(data)) {
@@ -325,7 +334,7 @@ SplitDim <- function(data, split_dim = 'time', indices, freq = 'monthly',
     }
     dates_exist <- TRUE
   }
-
+  
   # Return objects
   if (all(dates_exist, return_indices)) {
     return(list(data = data, dates = dates, indices = indices))
@@ -334,7 +343,7 @@ SplitDim <- function(data, split_dim = 'time', indices, freq = 'monthly',
   } else if (all(!dates_exist, return_indices)) {
     return(list(data = data, indices = indices))
   } else {
-    return(data)
+    return(data = data)
   }
 }
 
diff --git a/tests/testthat/test-CST_SplitDim.R b/tests/testthat/test-CST_SplitDim.R
index f8ad88ee09329e6106b630359643641d79b65f63..f461de52fdc1adca34cd995118d396163afffd14 100644
--- a/tests/testthat/test-CST_SplitDim.R
+++ b/tests/testthat/test-CST_SplitDim.R
@@ -8,9 +8,13 @@ class(data1) <- 's2dv_cube'
 
 indices1 <- c(rep(1,5), rep(2,5), rep (3, 5), rep(4, 5))
 output1 <- matrix(data1$data, nrow = 5, ncol = 4)
+output1 <- matrix(data1$data, nrow = 5, ncol = 4)
 names(dim(output1)) <- c('time', 'monthly')
 output1 <- list(data = output1) 
 output1$dims <- dim(output1$data)
+output1$coords <- list(time = 1:5, monthly = 1:4)
+attr(output1$coords$time, 'indices') <- TRUE
+attr(output1$coords$monthly, 'indices') <- TRUE
 class(output1) <- 's2dv_cube'
 
 exp_cor <- 1 : 20
@@ -23,6 +27,9 @@ output2 <- matrix(data1$data, nrow = 5, ncol = 4)
 names(dim(output2)) <- c('time', 'index')
 output2 <- list(data = output2)
 output2$dims <- dim(output2$data)
+output2$coords <- list(time = 1:5, index = 1:4)
+attr(output2$coords$time, 'indices') <- TRUE
+attr(output2$coords$index, 'indices') <- TRUE
 class(output2) <- 's2dv_cube'
 
 time2 <- c(seq(ISOdate(1903, 1, 1), ISOdate(1903, 1, 4), "days"),
@@ -44,7 +51,9 @@ dim(output3) <- c(time = 8, monthly = 3)
 result3 <- data3
 result3$data <- output3
 result3$dims <- dim(result3$data)
-
+result3$coords <- list(time = 1:8, monthly = 1:3)
+attr(result3$coords$time, 'indices') <- TRUE
+attr(result3$coords$monthly, 'indices') <- TRUE
 # dat4
 data4 <- list(data = array(rnorm(10), dim = c(sdate = 2, lon = 5)))
 class(data4) <- 's2dv_cube'
@@ -107,6 +116,10 @@ test_that("3. Output checks: sample data", {
   result$dims <- c(dataset = 1, member = 15, sdate = 6, ftime = 1, lat = 22, 
                    lon = 53, monthly = 3)
   attributes(result$attrs$Dates)$end <- NULL
+  result$coords$ftime <- 1 
+  result$coords$monthly <- 1:3
+  attr(result$coords$ftime, 'indices') <- TRUE
+  attr(result$coords$monthly, 'indices') <- TRUE
   expect_equal(
     CST_SplitDim(data = lonlat_temp$exp, split_dim = 'ftime', ftime_dim = 'ftime'),
     result