library(RNetCDF)

outfile <- paste0(workingDir,"/skill_forecast_output_", SEASON, "_lead", iLead, FTYPEF,".nc")

#--------------------------------------------------
# Basic sizes
#--------------------------------------------------
nlon    <- length(lon)
nlat    <- length(lat)
ntest   <- dim(best_probs_by_model)[3]
ncat    <- dim(best_probs_by_model)[4]
nmodels <- dim(best_probs_by_model)[5]
nrun    <- dim(test_probs_all)[5]

# Better to store actual years if available
test_vals <- if (exists("test_years")) as.integer(test_years) else as.integer(seq_len(ntest))
run_vals  <- if (exists("run_ids"))    as.integer(run_ids)    else as.integer(seq_len(nrun))
cat_vals  <- 1:ncat
cat_names <- c("below", "normal", "above")
model_vals <- seq_len(nmodels)
if (!exists("model_names")) model_names <- paste0("model", model_vals)

fill_dbl <- -9999.0
fill_int <- -9999L

#--------------------------------------------------
# Create file
#--------------------------------------------------
if (file.exists(outfile)) file.remove(outfile)
nc <- create.nc(outfile)

#--------------------------------------------------
# Define dimensions
#--------------------------------------------------
dim.def.nc(nc, "lon",   nlon)
dim.def.nc(nc, "lat",   nlat)
dim.def.nc(nc, "test",  ntest)
dim.def.nc(nc, "cat",   ncat)
dim.def.nc(nc, "run",   nrun)
dim.def.nc(nc, "model", nmodels)

dim.def.nc(nc, "nchar_model", max(nchar(model_names)))
dim.def.nc(nc, "nchar_cat",   max(nchar(cat_names)))

#--------------------------------------------------
# Define coordinate variables
#--------------------------------------------------
var.def.nc(nc, "lon",   "NC_DOUBLE", "lon")
var.def.nc(nc, "lat",   "NC_DOUBLE", "lat")
var.def.nc(nc, "test",  "NC_INT",    "test")
var.def.nc(nc, "cat",   "NC_INT",    "cat")
var.def.nc(nc, "run",   "NC_INT",    "run")
var.def.nc(nc, "model", "NC_INT",    "model")

var.def.nc(nc, "model_name", "NC_CHAR", c("nchar_model", "model"))
var.def.nc(nc, "cat_name",   "NC_CHAR", c("nchar_cat", "cat"))

#--------------------------------------------------
# Define data variables
#--------------------------------------------------
var.def.nc(nc, "best_run_by_model",   "NC_INT",    c("lon", "lat", "model"))
var.def.nc(nc, "best_hhss_by_model",  "NC_DOUBLE", c("lon", "lat", "model"))
var.def.nc(nc, "best_probs_by_model", "NC_DOUBLE", c("lon", "lat", "test", "cat", "model"))
var.def.nc(nc, "test_probs_all",      "NC_DOUBLE", c("lon", "lat", "test", "cat", "run", "model"))
var.def.nc(nc, "ens_trpss",           "NC_DOUBLE", c("lon", "lat"))

# Optional, since you also have them
var.def.nc(nc, "hhss_all",            "NC_DOUBLE", c("lon", "lat", "run", "model"))
var.def.nc(nc, "hsharpness_all",      "NC_DOUBLE", c("lon", "lat", "run", "model"))

#--------------------------------------------------
# Attributes
#--------------------------------------------------
att.put.nc(nc, "NC_GLOBAL", "title",   "NC_CHAR", "Forecast skill outputs")
att.put.nc(nc, "NC_GLOBAL", "season",  "NC_CHAR", as.character(SEASON))
att.put.nc(nc, "NC_GLOBAL", "lead",    "NC_INT",  as.integer(iLead))
att.put.nc(nc, "NC_GLOBAL", "history", "NC_CHAR",
           paste("Created", format(Sys.time(), "%Y-%m-%d %H:%M:%S")))

att.put.nc(nc, "lon", "units", "NC_CHAR", "degrees_east")
att.put.nc(nc, "lat", "units", "NC_CHAR", "degrees_north")

att.put.nc(nc, "best_run_by_model",   "missing_value", "NC_INT",    fill_int)
att.put.nc(nc, "best_hhss_by_model",  "missing_value", "NC_DOUBLE", fill_dbl)
att.put.nc(nc, "best_probs_by_model", "missing_value", "NC_DOUBLE", fill_dbl)
att.put.nc(nc, "test_probs_all",      "missing_value", "NC_DOUBLE", fill_dbl)
att.put.nc(nc, "ens_trpss",           "missing_value", "NC_DOUBLE", fill_dbl)
att.put.nc(nc, "hhss_all",            "missing_value", "NC_DOUBLE", fill_dbl)
att.put.nc(nc, "hsharpness_all",      "missing_value", "NC_DOUBLE", fill_dbl)

att.put.nc(nc, "best_hhss_by_model", "long_name", "NC_CHAR", "Best hindcast HSS by model")
att.put.nc(nc, "best_probs_by_model","long_name", "NC_CHAR", "Forecast probabilities from best run")
att.put.nc(nc, "test_probs_all",     "long_name", "NC_CHAR", "Forecast probabilities for all runs")
att.put.nc(nc, "best_run_by_model",  "long_name", "NC_CHAR", "Best run index by model")
att.put.nc(nc, "ens_trpss",          "long_name", "NC_CHAR", "Ensemble tercile RPSS")
att.put.nc(nc, "hhss_all",           "long_name", "NC_CHAR", "Hindcast HSS for all runs")
att.put.nc(nc, "hsharpness_all",     "long_name", "NC_CHAR", "Hindcast sharpness for all runs")

#--------------------------------------------------
# Write coordinates
#--------------------------------------------------
var.put.nc(nc, "lon", lon)
var.put.nc(nc, "lat", lat)
var.put.nc(nc, "test", test_vals)
var.put.nc(nc, "cat", as.integer(cat_vals))
var.put.nc(nc, "run", run_vals)
var.put.nc(nc, "model", as.integer(model_vals))

#--------------------------------------------------
# Helper for char vars
#--------------------------------------------------
#write_char_var <- function(nc, varname, strings, strlen_dim) {
#  strlen <- dim.inq.nc(nc, strlen_dim)$length
#  mat <- matrix(" ", nrow = strlen, ncol = length(strings))
#  for (i in seq_along(strings)) {
#    s <- sprintf(paste0("%-", strlen, "s"), strings[i])
#    mat[, i] <- strsplit(substr(s, 1, strlen), "")[[1]]
#  }
#  var.put.nc(nc, varname, mat)
#}

write_char_var <- function(nc, varname, strings, strlen_dim) {
  strlen <- dim.inq.nc(nc, strlen_dim)$length
  vinfo  <- var.inq.nc(nc, varname)

  nd <- length(vinfo$dimids)
  if (nd != 2) stop(varname, " must be a 2-D char variable")

  dim1_name <- dim.inq.nc(nc, vinfo$dimids[1])$name
  dim2_name <- dim.inq.nc(nc, vinfo$dimids[2])$name

  if (dim1_name == strlen_dim) {
    mat <- matrix(" ", nrow = strlen, ncol = length(strings))
    for (i in seq_along(strings)) {
      s <- sprintf(paste0("%-", strlen, "s"), strings[i])
      mat[, i] <- strsplit(substr(s, 1, strlen), "", fixed = TRUE)[[1]]
    }
  } else if (dim2_name == strlen_dim) {
    mat <- matrix(" ", nrow = length(strings), ncol = strlen)
    for (i in seq_along(strings)) {
      s <- sprintf(paste0("%-", strlen, "s"), strings[i])
      mat[i, ] <- strsplit(substr(s, 1, strlen), "", fixed = TRUE)[[1]]
    }
  } else {
    stop("Could not find strlen_dim in variable dimensions")
  }

  var.put.nc(nc, varname, mat, start = rep(1, 2), count = dim(mat))
}

write_char_var(nc, "model_name", model_names, "nchar_model")
write_char_var(nc, "cat_name",   cat_names,   "nchar_cat")

#--------------------------------------------------
# Replace NA with fill values
#--------------------------------------------------
x_best_run <- best_run_by_model
x_best_run[is.na(x_best_run)] <- fill_int
storage.mode(x_best_run) <- "integer"

x_best_hhss <- best_hhss_by_model
x_best_hhss[is.na(x_best_hhss)] <- fill_dbl

x_best_probs <- best_probs_by_model
x_best_probs[is.na(x_best_probs)] <- fill_dbl

x_test_probs <- test_probs_all
x_test_probs[is.na(x_test_probs)] <- fill_dbl

x_ens_trpss <- ens_trpss
x_ens_trpss[is.na(x_ens_trpss)] <- fill_dbl

x_hhss_all <- hhss_all
x_hhss_all[is.na(x_hhss_all)] <- fill_dbl

x_hsharpness_all <- hsharpness_all
x_hsharpness_all[is.na(x_hsharpness_all)] <- fill_dbl

#--------------------------------------------------
# Write arrays
#--------------------------------------------------
var.put.nc(nc, "best_run_by_model",   x_best_run)
var.put.nc(nc, "best_hhss_by_model",  x_best_hhss)
var.put.nc(nc, "best_probs_by_model", x_best_probs)
var.put.nc(nc, "test_probs_all",      x_test_probs)
var.put.nc(nc, "ens_trpss",           x_ens_trpss)
var.put.nc(nc, "hhss_all",            x_hhss_all)
var.put.nc(nc, "hsharpness_all",      x_hsharpness_all)

close.nc(nc)

cat("Wrote file:", outfile, "\n")


