Extracting max bottom depth in NetCDF files

Question

I have global 4 dimenions (lon, lat, time and depth) with two variables NetCDF files for environmental marine variables. I am looking to extract data from the bottom-most depth for each cell that does not contain an NA, and calculate a mean over all years.

You'll find attached my current code but that for level = 1 (surface depth of -1m), but I don't know how to do for the max bottom depth of each cell...

Do you know how I could do that in R (not on CDO) ?

library(raster)
library(ncdf4)


setwd("E:/Modélisation - Eric Goberville/PhD/Environmental data/Raw data/Western_MedSea/Merged")

#CMCC files

rm(list=ls())
graphics.off()

# Load .nc file and access informaiton
Info.Data <- ncdf4::nc_open("med-cmcc-tem-rean-m_01_2000_01_2021.nc")
Data_date <- ncvar_get(Info.Data, "time") # to have an idea of time span
head(Data_date)
#Data_date <- as.POSIXct("1900-01-01 00:00:00", tz = "UTC") + as.difftime(Data_date, units = "mins")#to have the time in a standard date format
print(Data_date)
Data_depth <- ncvar_get(Info.Data, "depth") # How many depth... use for level

# Extract data
names(Info.Data[["var"]])
Data_raster <- brick("med-cmcc-tem-rean-m_01_2000_01_2021.nc", varname = "bottomT", level = 1)
# Data_raster

# Calculate annual mean from monthly means
# A quick way, as we know time span
#Matrix creation  creating sequences of indices by starting at 1 and incrementing by 12,
ID <- cbind(seq(1, dim(Data_date), by = 12), ID <- seq(12, dim(Data_date), by = 12))

Data_raster_annual <- stack() #Empty rasterstack to store the annual means
for (i in 1:nrow(ID)) {
  Data_raster_annual <- stack(Data_raster_annual,
                              calc(Data_raster[[ID[i,1]:ID[i,2]]], mean)) #mean for each cell across the layers corresponding to the current 12-month period
}
# Calculate climatology (from annual means)
Raster_climatology <- calc(Data_raster_annual, mean)
#Raster_climatology

# Mapping of your environmental parameters 
graphics.off()
my.colors = colorRampPalette(c("#5E85B8","#EDF0C0","#C13127")) 
plot(Raster_climatology, main = "Bottom temperature",col=my.colors(1000),axes=FALSE, box=FALSE)