Work with MOD07 - Retrieved Moisture Profile layer - data in R

I downloaded a complete month of daily images of MODIS MOD07_L2 product (Atmospheric profile). I’ve worked with *.hdf files for really long (using both R and tools such as HEG), nevertheless, I am now stuck with MOD07 product. I need to work with Dew-point Temperature at different pressure levels, but I don’t understand how the data is actually saved in the *.hdf layer (which is Retrieved_Moisture_Profile). I tried using HEG software with all possible options (for output file types), but none of them worked. For every iteration I tried in HEG I got the "Could not find TimeofDay attribute, CalendarDate, Sensorname, pointingangle (and many more) attribute in input hdf file". So I believe HEG is not capable to work with this type of data. The result for every image was a black square with just one pixel value.

I opened the *.hdf file directly in QGIS (version 3.10.4-A Coruña) and, more specifically, I opened the Retrieved_Moisture_Profile layer. When I studied the histogram of the latter image it is quite clear that there are 20 different Dew-point temperature profiles (each of one of them is for an specific pressure level), as you can see in the following figure.

I was thinking that maybe I can just save the *.hdf to *.nc, but I tried with rhdf5 and the latter package is depricated for R version 3.6.3. I also tried opening the *.hdf file with gdalUtils, which actually works, but I couldn’t find how to ‘divide’ each layer of the raster or even save them separately.

So, what I need is the following (in R).

1. I need to open the MOD07 *.hdf file
2. I need to open the Retrieved_Moisture_Profile layer within the *.hdf file
3. I need to save the dew-point temperature profile for each pressure level i.e. for each ‘band’ (from 1 to 20), as a GeoTIFF file.

FYI: I’m running R version 3.6.3 in Windows 10 64 bit machine.
Thanks in advance for your help!
Sincerely, Jorge.

library(raster)
library(rgdal)
library(gdalUtils)

# hdf files directory
main_dir <- file.path("D:/path_to_your_hdf_files...")
# List all *.hdf files within main_dir
hdf_files <- list.files(path = main_dir, pattern = "*.hdf", full.names = T)

# Get the info of the first hdf within the hdf_files list
hdf_info <- GDALinfo(hdf_files[1], returnScaleOffset = FALSE)
# Obtain the metadata of the latter hdf file
meta <- attr(hdf_info, "mdata")
# Name the extent name attributes to search
crd_str <- paste0(c("WEST", "EAST", "SOUTH", "NORTH"), "BOUNDINGCOORDINATE")
# Function to search for the attributes
crd_id <- sapply(crd_str, function(i) grep(i, meta))
# Extract the coordinates of each bounding coordinate
crd <- meta[crd_id]
crd <- as.numeric(sapply(strsplit(crd, "="), "[[", 2))
# Define the coordinates as a raster::extent to later add to the raster
ext <- extent(crd)

# Open the first hdf file
one <- getSds(hdf_files[1])
# Open band 16 (which is Retrieved Moisture Profile) as a raster brick
b16 <- brick(readGDAL(one$SDS4gdal[16], as.is = TRUE))
# Set crd as the extent of b16
extent(b16) <- crd
# Set EPSG:4326 as the projection of b16
projection(b16) <- "+init=epsg:4326"