Funções do aplicativo Raster com várias pilhas e índice compartilhado

Tenho duas pilhas raster, ambas representando uma sequência de meses. Stack1 é uma sequência repetida em vários anos (por exemplo, Jan, Fev, Março, ..., Jan, Fev, Março, ..., Jan, Fev, Março, ...) e Stack2 é uma única sequência de valores de referência para esses mesmos meses (por exemplo, Jan, Fev, Março, ...). Quero executar uma operação combinando as duas pilhas usando seu índice comum - por exemplo, Stack1/Stack2, mas específico para o índice (então Stack1Jan1/Stack2Jan, Stack1Feb1/Stack2Feb, ... Stack1Jan2/Stack2Jan).

Tenho quase certeza de que isso é possível usando um dos métodos baseados em apply (app, tapp, rapp, etc.), mas apesar de gastar muito tempo na documentação, não consigo descobrir qual. Sei que usando tapp, posso aplicar uma função a diferentes camadas do raster usando um índice, como no exemplo abaixo, mas não tenho certeza de como fazer esse mesmo índice aplicar ao envolver uma segunda pilha.

#sample data
r <- rast(ncols=9, nrows=9)
values(r) <- 1:ncell(r)
s <- c(r, r, r, r, r, r)
s <- s * 1:6
Stack1 <- c(s, s, s)
Stack2 <- rast(ncols=9, nrows=9, nlyrs=6)
values(Stack2) <- rep(1, ncell(Stack2)*6)
Stack2 <- Stack2 * 1:6

#tapp lets you apply an index for layers within a single stack
idx <- rep(1:6, 3)
s1_mean <- tapp(Stack1, idx, fun=mean)

# what I'd like would be something like:
myfun <- function(x) {x/Stack2[[idx]]}
s1_out <- tapp(Stack1, index=idx, fun=myfun)

Eu sei que poderia fazer isso subdividindo cada raster em seus meses e então fazendo o cálculo dessa forma, mas tenho certeza de que há uma maneira mais eficiente!