How to put 2-D fitted SPI6 back to 3-D array, add lon and lat dimensions and write it as nc file

Question

I have Netcdf data file called "P_monthly.nc", which contains monthly CHIRPS rainfall data (1981-2022). The file is light (just 49 MB). I have reshaped the array to matrix to calculate SPI6 using SPEIpackage. I then calculated SPI6 by splitting data into chunks to parallelize using snowfall package. I got some hints from some discussions to write the my code, like here and here.

I want to put the fitted SPI6, which contain a 2-D matrix, back to a 3-D array and add other dimensions like lon and lat, and write it as .ncfile.

Can anyone help please?

I used the code below:

rm(list=ls(all=TRUE))

# Load required packages
require(utils)
require(RNetCDF)
require(rasterVis)
require(rgdal)
library(ncdf4)
library(metR)
require(SPEI)
library(snowfall)

###############################################################################
# Upload and preprocess data
prec <- nc_open("data/chirps/P_monthly.nc")
lon <- ncvar_get(prec,"longitude")
lat <- ncvar_get(prec,"latitude")
nlon <- length(lon)
nlat <- length(lat)
ppt <- ncvar_get(prec, "Precip")
dim(ppt) # 138 165 504 (1981-2022)

#Reshape data to matrix
precip <- sapply(1:dim(ppt)[3], function(x) c(ppt[,,x]))
precip <- t(precip)
dim(precip) # nice, 504 x 22770

# Initialize a cluster using `snowfall` package.
cores <- 6
print(cores)
sfInit(parallel=TRUE, cpus=cores)
sfLibrary(SPEI)

# Create a matrix to store the results
spi6 <- matrix(NA, nrow = nrow(precip), ncol = ncol(precip))

# Calculate SPI-6 by splitting data into chunks
system.time(
  for (i in 0:(cores - 1)) {
    chunk <- ((i * ncol(precip) %/% cores) + 1) : ((i + 1) * ncol(precip) %/% cores)
    spi6[, chunk] <- sfApply(precip[, chunk], 2, function(x) spi(x, 6, distribution = "Gamma", na.rm = TRUE)$fitted)
  }
)

sfStop()

class(spi6)
dim(spi6) #504 22770

My P_monthly.nc file can be found here. It is light file that can be downloaded (about 49 MB)