Parallel computing

General
1

Dear all, I am trying to develop a parallel computing code in R, I know how to do grid by grid, but its my first time running codes in parallel. can you please check my code, if this looks ok.

library(ncdf4)
library(trend)
library(parallel)
library(doParallel)
library(foreach)

calc_trend <- function(infile, outfile, spinup = 6, ncores = 40) {
# Read only coordinates in master
  nc  <- nc_open(infile)
  lon <- ncvar_get(nc, "lon")
  lat <- ncvar_get(nc, "lat")
  ny  <- length(lat)
  nt  <- dim(ncvar_get(nc, "tavg", start=c(1,1,1), count=c(1,1,-1)))[1]  # probe nt cheaply
  nc_close(nc)

  # Re-read nt properly
  nc  <- nc_open(infile)
  nt  <- dim(nc$var[["tavg"]]$varsize)[3]
  nc_close(nc)

  nx <- length(lon)

  cat("Dataset:", nx, "x", ny, "x", nt, "\n")
  cat("Cores  :", ncores, "\n")

  cl <- makeCluster(ncores)
  registerDoParallel(cl)

  results <- foreach(
    i         = 1:nx,
    .packages = c("ncdf4", "trend", "modifiedmk"),
    .combine  = "rbind",
    .inorder  = TRUE
  ) %dopar% {

    slope_row <- rep(NA_real_, ny)
    tau_row   <- rep(NA_real_, ny)
    pval_row  <- rep(NA_real_, ny)

    tryCatch({
      nc_w     <- nc_open(infile)
      Rr_strip <- ncvar_get(nc_w, "tavg",
                            start = c(i, 1, 1),
                            count = c(1, ny, -1))
      nc_close(nc_w)

      for (j in 1:ny) {
        Rr <- as.numeric(Rr_strip[j, ])
        Rr <- Rr[-(1:spinup)]
        Rr <- Rr[is.finite(Rr)]

        if (length(Rr) > 10 && sd(Rr) > 0) {
          tryCatch({
            slp          <- sens.slope(Rr)
            slope_row[j] <- slp$estimates * 12
            mk_res        <- mmkh(Rr)
            tau_row[j]   <- mk_res["Tau"][[1]]
            pval_row[j]  <- mk_res["new P-value"][[1]]
          }, error = function(e) {})
        }
      }
    }, error = function(e) {
      cat("Strip error at i=", i, ":", conditionMessage(e), "\n")
    })

    c(slope_row, tau_row, pval_row)
  }

  stopCluster(cl)

  # Rebuild arrays
  slope  <- matrix(results[, 1:ny],              nrow = nx)
  mk_tau <- matrix(results[, (ny+1):(2*ny)],     nrow = nx)
  pvalue <- matrix(results[, (2*ny+1):(3*ny)],   nrow = nx)

  # Write NetCDF
  lon_dim <- ncdim_def("lon", "degrees_east",  lon)
  lat_dim <- ncdim_def("lat", "degrees_north", lat)

  var1 <- ncvar_def("sen_slope", "tavg/year", list(lon_dim, lat_dim), NA, compression=5)
  var2 <- ncvar_def("mk_tau",    "1",          list(lon_dim, lat_dim), NA, compression=5)
  var3 <- ncvar_def("pvalue",    "1",          list(lon_dim, lat_dim), NA, compression=5)

  ncnew <- nc_create(outfile, list(var1, var2, var3))
  ncvar_put(ncnew, var1, slope)
  ncvar_put(ncnew, var2, mk_tau)
  ncvar_put(ncnew, var3, pvalue)
  nc_close(ncnew)

  cat("Finished:", outfile, "\n")
}

calc__trend(
  "Global_temperature.nc",
  "Global_temperature_trend.nc",
  spinup = 6,
  ncores = 10
)