library(ncdf4)

#----------------
# Make dimensions
#----------------
xlon <- as.numeric(lon)
xlat <- as.numeric(lat)
  
nx <- length(xlon)
ny <- length(xlat)

xdim <- ncdim_def( 'lon', 'degreesE', xlon )
ydim <- ncdim_def( 'lat', 'degreesN', xlat )
 
#---------
# Define the variables 
#---------
missing_val <- -9999     # missing value

pblw <- ncvar_def( 'below', '%', list(xdim,ydim), missing_val )
pabv <- ncvar_def( 'above', '%', list(xdim,ydim), missing_val )
pnor <- ncvar_def( 'normal','%', list(xdim,ydim), missing_val )
hrps <- ncvar_def( 'hrpss', '%', list(xdim,ydim), missing_val )
trps <- ncvar_def( 'trpss', '%', list(xdim,ydim), missing_val )
thss <- ncvar_def( 'hhss', '%', list(xdim,ydim), missing_val )

pnsel <- ncvar_def(
  "n_selected",
  "count",
  list(xdim, ydim),
  missing_val,
  longname = "number of predictors selected in final training model",
  prec = "integer"
)

output_fname <- paste0(dirdat, fileName)

ncid <- nc_create(
  output_fname,
  list(pblw, pabv, pnor, hrps, trps, thss, pnsel)
)


#-------------------------------
# Write data 
#-------------------------------

ncvar_put( ncid, pblw, below, start=c(1,1), count=c(nx,ny))
ncvar_put( ncid, pabv, above, start=c(1,1), count=c(nx,ny))
ncvar_put( ncid, pnor, normal, start=c(1,1), count=c(nx,ny))
ncvar_put( ncid, hrps, hrpss  , start=c(1,1), count=c(nx,ny))
ncvar_put( ncid, trps, trpss  , start=c(1,1), count=c(nx,ny))
ncvar_put( ncid, thss, hhss  , start=c(1,1), count=c(nx,ny))


n_selected_out <- n_selected
n_selected_out[is.na(n_selected_out)] <- missing_val
storage.mode(n_selected_out) <- "integer"

ncvar_put(ncid, pnsel, n_selected_out, start = c(1,1), count = c(nx,ny))

nc_close( ncid )


