How can I add the time dimension to a series of .tif files?

I have been trying to arrange a series of .tif files into a stars object, without success. I have three climatic variables: temperature (temp), precipitation (prec), and maximum temperature (tmax) downloaded from WorldClim using the following function

# devtools required
if(!require(devtools)) install.packages("devtools")
library(devtools)
devtools::install_github("HelgeJentsch/ClimDatDownloadR")
library(ClimDatDownloadR)

# shapefile and bounding box
aoi_abruzzo <- st_read("abruzzo.shp") %>% .$geometry 

# Bounding Box 
abruzzo_bb <- st_bbox(aoi_abruzzo)

WorldClim.HistClim.download(
  # save.location = "./",
  parameter = c("prec", "temp", "tmax"),
  
  # 4 months, from jan to april
  month.var = c(1:4),

  # here, 10 arc-minutes are chosen as resolution 
  resolution = c("10m"),
  version.var = c("2.1"),
  clipping = TRUE,
  clip.shapefile = aoi_abruzzo,
  clip.extent = abruzzo_bb,
  # buffer = 0,
  convert.files.to.asc = FALSE,
  stacking.data = TRUE,
  keep.raw.zip = FALSE,
  delete.raw.data = FALSE,
  save.bib.file = TRUE
)

# Import raster
# String containing the names of raster files
rastlist <- list.files(path ="my/path/prec/WorldClim_2.1_prec_30s/clipped_2024-10-09_16-35-06", pattern = "wc2.1", full.names = TRUE)

# Using the list of names, all the files are imported into a single raster package
stack_prec <- stack(rastlist)

I obtain three folders: prec, temp, and tmax. Inside each of these, there are 4 .tif files that correspond to the monthly average of the variable for the months of January, February, March, and April. By stacking each variable, I get an object like this:

> stack_prec
class      : RasterStack 
dimensions : 145, 212, 30740, 4  (nrow, ncol, ncell, nlayers)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 13.01667, 14.78333, 41.68333, 42.89167  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : wc2.1_30s_prec_01, wc2.1_30s_prec_02, wc2.1_30s_prec_03, wc2.1_30s_prec_04 
min values :                37,                34,                35,                32 
max values :               102,                99,                84,               106

Each layer represents a different month. I would like to convert it to a stars object to include the time dimension. However, I am unable to match the time dimension to each layer. This is what I have done, following this response this:

# stack into stars
stars_prec <- st_as_stars(stack_prec)

> stars_prec
stars object with 3 dimensions and 1 attribute
attribute(s):
                   Min. 1st Qu. Median     Mean 3rd Qu. Max.  NA's
wc2.1_30s_prec_01    32      48     53 55.19539      60  106 24480
dimension(s):
     from  to offset     delta                       refsys                                  values x/y
x       1 212  13.02  0.008333 +proj=longlat +datum=WGS8...                                    NULL [x]
y       1 145  42.89 -0.008333 +proj=longlat +datum=WGS8...                                    NULL [y]
band    1   4     NA        NA                          NA wc2.1_30s_prec_01,...,wc2.1_30s_prec_04 

Already at this point, I can see that the problem is that it only considers the first month of precipitation (wc2.1_30s_prec_01) as an 'attribute' and does not include all 4 of them. To add the temporal dimension, I tried the following steps

# add temporal dimension by repeating 4 times 
l = vector("list", 4)
l[] = list(stars_prec)
r = do.call("c", l) 
r = st_redimension(r)                                                                                                            
r <- st_set_dimensions(.x = r, which = "attributes", values = as.Date(c("1980-01-01", "1980-02-01", "1980-03-01","1980-04-01")), names = "time") 

And obtaining this result

> r
stars object with 4 dimensions and 1 attribute
attribute(s):
                                Min. 1st Qu. Median     Mean 3rd Qu. Max.  NA's
wc2.1_30s_tavg_01_clipped.t...  -7.2     2.3    5.1 4.930197     7.4 13.6 97920
dimension(s):
        from  to offset     delta refsys point                                                            values x/y
x          1 212  13.02  0.008333 WGS 84 FALSE                                                              NULL [x]
y          1 145  42.89 -0.008333 WGS 84 FALSE                                                              NULL [y]
time       1   4     NA        NA   Date    NA                                         1980-01-01,...,1980-04-01    
new_dim    1   4     NA        NA     NA    NA wc2.1_30s_tavg_01_clipped.tif,...,wc2.1_30s_tavg_01_clipped.tif.3 

I set the temporal dimension, but the attribute still only includes the precipitation values for the first month.

If I apply the split function on the "new_dim" dimension, I get this:

> split(r,3)
stars object with 3 dimensions and 4 attributes
attribute(s):
                   Min. 1st Qu. Median     Mean 3rd Qu. Max.  NA's
wc2.1_30s_prec_01    37      44     48 49.30158      51  102 24480
wc2.1_30s_prec_02    34      49     52 53.25581      56   99 24480
wc2.1_30s_prec_03    35      50     53 53.35853      56   84 24480
wc2.1_30s_prec_04    32      61     66 64.86564      70  106 24480
dimension(s):
     from  to offset     delta                       refsys                    values x/y
x       1 212  13.02  0.008333 +proj=longlat +datum=WGS8...                      NULL [x]
y       1 145  42.89 -0.008333 +proj=longlat +datum=WGS8...                      NULL [y]
time    1   4     NA        NA                         Date 1980-01-01,...,1980-04-01  

Each month becomes a separate attribute. Ultimately, I would like to achieve something like this with all the variables (I apologize for the somewhat rough way of describing the result):

and be able to select them simply by taking the date.

I tried to incorporate the months of the same variable into a stack and associate a date with each layer, but it didn't work. The time dimension seems to exist on its own and is not correctly associated with each layer

Edit: shapefile and bounding box as example

# download shapefile
install.packages("devtools")
devtools::install_github("ropensci/rnaturalearthhires")
library(rnaturalearth)

map1 <- ne_countries(type = "countries", country = "Malta",
                     scale = "medium", returnclass = "sf")
malta <- map1$geometry
malta_bb <- st_bbox(malta)

Example data

library(geodata)
tmin = worldclim_country("Lux", "tmin", path=tempdir())
tmax = worldclim_country("Lux", "tmax", path=tempdir())
prec = worldclim_country("Lux", "prec", path=tempdir())

Set time (not entirely clear what type of time you have in mind, so here I use "raw" values of 1 to 12)

time(tmin) <- time(tmax) <- time(prec) <- 1:12

With "terra" you might combine these so:

s <- sds(tmin, tmax, prec)
s
#class       : SpatRasterDataset 
#subdatasets : 3 
#dimensions  : 180, 180 (nrow, ncol)
#nlyr        : 12, 12, 12 
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : 5.5, 7, 49, 50.5  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#source(s)   : LUX_wc2.1_30s_tmin.tif, LUX_wc2.1_30s_tmax.tif, LUX_wc2.1_30s_prec.tif 
#names       : LUX_wc2.1_30s_tmin, LUX_wc2.1_30s_tmax, LUX_wc2.1_30s_prec 

You could then get the second month for all variables like this

s[,2]
#or more generally
s[, which(time(s)[[1]] == 2)]

To make a stars object you can do

x <- list(tmin, tmax, prec)
do.call("c", lapply(x, stars::st_as_stars))

#stars object with 3 dimensions and 3 attributes
#attribute(s):
#                        Min. 1st Qu. Median      Mean 3rd Qu.  Max.
#LUX_wc2.1_30s_tmin.tif  -3.2     0.4    4.4  4.917869     9.4  14.1
#LUX_wc2.1_30s_tmax.tif   0.9     5.9   12.5 12.636610    18.9  25.1
#LUX_wc2.1_30s_prec.tif  46.0    67.0   75.0 76.737418    84.0 143.0
#dimension(s):
#     from  to offset     delta refsys point x/y
#x       1 180    5.5  0.008333 WGS 84 FALSE [x]
#y       1 180   50.5 -0.008333 WGS 84 FALSE [y]
#time    1  12      1         1     NA    NA