Set up FOCUS PELMO runs and run them in parallel

- This works on Linux using wine
- PELMO runs (including pelmo.inp files) are correctly generated
- The PLM files for FOCUS Pesticide_D in the test data archive are
  correctly reproduced
- The data files (including FOCUS groundwater scenario data) are now
  created and documented in R files

3 files changed, 446 insertions, 2 deletions

diff --git a/R/FOCUS_GW_scenarios_2012.R b/R/FOCUS_GW_scenarios_2012.R
index 49cf7dd..e26ea50 100644
--- a/R/FOCUS_GW_scenarios_2012.R
+++ b/R/FOCUS_GW_scenarios_2012.R
@@ -1,11 +1,69 @@
 #' A very small subset of the FOCUS Groundwater scenario defitions
 #' 
-#' Currently, only a small subset of the soil definitions are provided.
+#' Currently, only scenario names with acronyms and a small subset of the soil definitions are provided. The 
+#' soil definitions are from page 46ff. from FOCUS (2012).
 #' 
 #' @name FOCUS_GW_scenarios_2012
 #' @references FOCUS (2012) Generic guidance for Tier 1 FOCUS ground water assessments. Version 2.1. 
 #'  FOrum for the Co-ordination of pesticde fate models and their USe. 
 #'  http://focus.jrc.ec.europa.eu/gw/docs/Generic_guidance_FOCV2_1.pdf
+#' @export
 #' @examples
 #' FOCUS_GW_scenarios_2012
-NULL
+FOCUS_GW_scenarios_2012 = list()
+
+n_layers = c(7, 6, 6, 5, 5, 6, 4, 6, 6)
+acronyms = c("Cha", "Ham", "Jok", "Kre", "Oke", "Pia", "Por", "Sev", "Thi")
+FOCUS_GW_scenarios_2012$names = c("Ch\u00e2teaudun", "Hamburg", "Jokioinen",
+                                  "Kremsm\u00fcnster", "Okehampton",
+                                  "Piacenza", "Porto", "Sevilla", "Thiva")
+
+names(FOCUS_GW_scenarios_2012$names) = acronyms
+FOCUS_GW_scenarios_2012$soils <- data.frame(
+  location= rep(acronyms, times = n_layers),
+  horizon = c("Ap", "B1", "B2", "II C1", "II C1", "II C2", "M",
+              "Ap", "BvI", "BvII", "Bv/Cv", "Cv", "Cv",
+              "Ap", "Bs", "BC1", "BC2", "BC2", "Cg",
+              rep(NA, 5),
+              "A", "Bw1", "BC", "C", "C",
+              "Ap", "Ap", "Bw", "Bw", "2C", "2C",
+              rep(NA, 4),
+              rep(NA, 6),
+              "Ap1", "Ap2", "Bw", "Bw", "Ck1", "Ck1"),
+  number = unlist(sapply(n_layers, function(x) 1:x)),
+  pH_H2O = c(8.0, 8.1, 8.2, 8.5, 8.5, 8.5, 8.3,
+             6.4, 5.6, 5.6, 5.7, 5.5, 5.5,
+             6.2, 5.6, 5.4, 5.4, 5.4, 5.3,
+             7.7, 7.0, 7.1, 7.1, 7.1,
+             5.8, 6.3, 6.5, 6.6, 6.6,
+             7, 7, 6.3, 6.3, 6.4, 6.4,
+             4.9, 4.8, 4.8, 4.8,
+             7.3, 7.3, 7.8, 8.1, 8.1, 8.2,
+             7.7, 7.7, 7.8, 7.8, 7.8, 7.8),
+  perc_clay = c(30, 31, 25, 26, 26, 24, 31,
+                7.2, 6.7, 0.9, 0, 0, 0,
+                3.6, 1.8, 1.2, 1.7, 1.7, 1.9,
+                14, 25, 27, 27, 27,
+                18, 17, 14, 9, 9,
+                15, 15, 7, 7, 0, 0,
+                10, 8, 8, 8,
+                14, 13, 15, 16, 16, 22,
+                25.3, 25.3, 29.6, 31.9, 32.9, 32.9),
+  perc_oc = c(1.39, 0.93, 0.7, 0.3, 0.3, 0.27, 0.21,
+              1.5, 1, 0.2, 0, 0, 0,
+              4.06, 0.84, 0.36, 0.29, 0.29, 0.21,
+              3.6, 1.0, 0.5, 0.5, 0.5,
+              2.2, 0.7, 0.4, 0.1, 0.1,
+              1.26, 1.26, 0.47, 0.47, 0, 0,
+              1.42, 0.78, 0.78, 0.78,
+              0.93, 0.93, 0.70, 0.58, 0.58, 0.49,
+              0.74, 0.74, 0.57, 0.31, 0.18, 0.18),
+  rel_deg = c(1, 0.5, 0.5, 0.3, 0, 0, 0,
+              1, 0.5, 0.3, 0.3, 0.3, 0,
+              1, 0.5, 0.3, 0.3, 0, 0,
+              1, 0.5, 0.5, 0.3, 0,
+              1, 0.5, 0.3, 0.3, 0,
+              1, 0.5, 0.5, 0.3, 0.3, 0,
+              1, 0.5, 0.3, 0,
+              1, 1, 0.5, 0.3, 0, 0,
+              1, 0.5, 0.5, 0.3, 0.3, 0))
diff --git a/R/FOCUS_PELMO_data.R b/R/FOCUS_PELMO_data.R
new file mode 100644
index 0000000..28888e7
--- /dev/null
+++ b/R/FOCUS_PELMO_data.R
@@ -0,0 +1,90 @@
+#' FOCUS PELMO crop names
+#'
+#' A named character vector with the crop names used in the PELMO 5.5.3 GUI.
+#' For the names, three letter codes were constructed by generally taking the
+#' first three letters in lower case. Only when there is an expression in 
+#' parentheses, the first letter of this expression becomes the first letter 
+#' in the three letter code, i.e. 'Peas (animals)' has the code \code{ape}.
+#' 
+#' @docType data
+#' @export
+#' @examples
+#' print(FOCUS_PELMO_crops)
+FOCUS_PELMO_crops = c(
+  app = "Apples",
+  gra = "Grass and alfalfa",
+  pot = "Potatoes",
+  sug = "Sugar beets",
+  win = "Winter cereals",
+  fbe = "Beans (field)",
+  vbe = "Beans (vegetables)",
+  bus = "Bushberries",
+  cab = "Cabbage",
+  car = "Carrots",
+  cit = "Citrus",
+  cot = "Cotton",
+  lin = "Linseed",
+  mai = "Maize",
+  soi = "Oil seed rape (summer)",
+  woi = "Oil seed rape (winter)",
+  oni = "Onions",
+  ape = "Peas (animals)",
+  soy = "Soybeans",
+  spr = "Spring cereals",
+  str = "Strawberries",
+  sun = "Sunflower",
+  tob = "Tobacco",
+  tom = "Tomatoes",
+  vin = "Vines")
+
+#' FOCUS PELMO crop acronyms used for naming sceneario files 
+#'
+#' A named character vector with the crop acronyms used in FOCUS PELMO 5.5.3
+#' for naming the .sze files located in the FOCUS directory. The crop acronyms
+#' in the file names are sometimes capitalized, sometimes not. The scenario
+#' files used for Beans (field) and Beans (vegetable) are the same.
+#' 
+#' @docType data
+#' @export
+#' @examples
+#' print(FOCUS_PELMO_crop_sze_names)
+FOCUS_PELMO_crop_sze_names = c(
+  app = "apples",
+  gra = "grass",
+  pot = "potato",
+  sug = "sbeets",
+  win = "wcerea",
+  fbe = "beans",
+  vbe = "beans", # Same sze as for fbe is used, with 'irrigation' which has no effect on the run
+  bus = "bushb",
+  cab = "cabbag",
+  car = "carrot",
+  cit = "citrus",
+  cot = "cotton",
+  lin = "linse",
+  mai = "maize",
+  soi = "rapesu",
+  woi = "rapewi",
+  oni = "onions",
+  ape = "peas",
+  soy = "soyb",
+  spr = "scerea",
+  str = "strawb",
+  sun = "sunflo",
+  tob = "tobacc",
+  tom = "tomato",
+  vin = "vines")
+
+#' Location codes in FOCUS PELMO
+#'
+#' A named character vector of one letter location codes used in FOCUS PELMO,
+#' indexed by three letter acronyms.
+
+#' @docType data
+#' @export
+#' @examples
+#' print(FOCUS_PELMO_location_codes)
+FOCUS_PELMO_location_codes = c(
+  Cha = "C", Ham = "H", Jok = "J", 
+  Kre = "K", Oke = "N", Pia = "P", 
+  Por = "O", Sev = "S", Thi = "T")
diff --git a/R/PELMO_runs.R b/R/PELMO_runs.R
new file mode 100644
index 0000000..7430b40
--- /dev/null
+++ b/R/PELMO_runs.R
@@ -0,0 +1,296 @@
+#' Set up runs for FOCUS PELMO
+#'
+#' Per default, the runs are also executed with FOCUS PELMO, and the results are processed
+#' and returned.
+#'
+#' @param runs A list of lists. Each inner lists has an element named 'psm'
+#'   that holds the psm string, and elements named using three letter crop acronyms,
+#'   as used in \code{\link{FOCUS_PELMO_crops}},
+#'   that hold character vectors of three letter scenario acronyms
+#'   as used in \code{\link{FOCUS_GW_scenarios_2012}}.
+#' @param psm_dir The directory where the psm files are located
+#' @param version The FOCUS PELMO version
+#' @param PELMO_base Where the FOCUS PELMO installation is located
+#' @param execute Should PELMO be executed directly?
+#' @param cores The number of cores to execute PELMO runs in parallel
+#' @param evaluate Should the results be returned?
+#' @param overwrite Should an existing run directories be overwritten
+#' @export
+PELMO_runs <- function(runs, psm_dir = ".", version = "5.5.3", PELMO_base = "auto",
+                       execute = TRUE, cores = getOption("mc.cores", 2L),
+                       evaluate = TRUE, overwrite = FALSE)
+{
+  if (PELMO_base[1] == "auto") {
+    PELMO_base = file.path(system.file(package = "PELMO.installeR"),
+                                       paste0("FOCUSPELMO.", gsub("\\.", "", version)))
+  }
+  PELMO_exe = list.files(PELMO_base, "^pelmo.*.exe", full.names = TRUE)
+  if (length(PELMO_exe) == 0) {
+    stop("Could not find PELMO executable. Did you run PELMO.installeR::install_PELMO()?")
+  }
+
+  run_list <- create_run_list(runs, psm_dir = psm_dir, check_psm_files = TRUE)
+
+  setup_run <- function(x) {
+    psm <- x$psm
+    crop <- x$crop
+    scenario <- x$scenario
+
+    location_code <- FOCUS_PELMO_location_codes[scenario]
+
+    run_dir <- file.path(PELMO_base, "FOCUS", PELMO_path(psm, crop, scenario))
+
+    if (dir.exists(run_dir)) {
+      if (overwrite) {
+        unlink(run_dir, recursive = TRUE)
+      } else {
+        stop("Run directory for ", crop, " in ", scenario, "\n", run_dir, "\nalready exists")
+      }
+    }
+
+    # Create the run directory
+    dir.create(run_dir, recursive = TRUE)
+
+    # Copy the psm file
+    psm_file <- file.path(psm_dir, paste0(psm, ".psm"))
+    file.copy(psm_file, run_dir)
+
+    # Copy Haude factors
+    file.copy(file.path(PELMO_base, "FOCUS", "HAUDE.DAT"), file.path(run_dir, "Haude.dat"))
+
+    # Copy scenario file
+    sze_file_upper <- paste0(location_code, "_",
+                             toupper(FOCUS_PELMO_crop_sze_names[crop]), ".sze")
+    sze_file_lower <- paste0(location_code, "_",
+                             FOCUS_PELMO_crop_sze_names[crop], ".sze")
+    if (file.exists(file.path(PELMO_base, "FOCUS", sze_file_lower))) {
+      sze_file <- sze_file_lower
+    } else {
+      if (file.exists(file.path(PELMO_base, "FOCUS", sze_file_upper))) {
+        sze_file <- sze_file_upper
+      } else {
+        stop("Could not find szenario file for ", crop, " in ", scenario)
+      }
+    }
+    file.copy(file.path(PELMO_base, "FOCUS", sze_file), file.path(run_dir, sze_file_lower))
+
+    # Generate PELMO input file
+    input_file <- file(file.path(run_dir, "pelmo.inp"), encoding = "latin1", open = "w+")
+    on.exit(close(input_file))
+
+    add <- function(x) cat(paste0(x, "\r\n"), file = input_file, append = TRUE)
+
+    # How many years do we calculate (26, 46 or 66)?
+    psm <- readLines(psm_file, encoding = "latin1")
+    number_of_apps_lines <- grep("number of application location", psm)
+    absolute_apps_line <- grep(location_code,
+                               psm[number_of_apps_lines])
+    period <- 1 # application every year
+    if (length(absolute_apps_line) == 1) {
+      apps_root <- number_of_apps_lines[absolute_apps_line]
+    } else {
+      apps_root <- number_of_apps_lines[1]
+    }
+
+    number_of_apps <- as.integer(substr(psm[apps_root], 1, 3))
+    last_app_line <- psm[apps_root + number_of_apps]
+    last_app_year <- as.integer(gsub("^.{2,3}  ..  (..)   .*", "\\1",
+                                     last_app_line))
+    if (last_app_year > 26) period <- 2
+    if (last_app_year > 46) period <- 3
+
+    n_years <- switch(as.character(period),
+      "1" = 26,
+      "2" = 46,
+      "3" = 66)
+
+    # First line with number of climate years
+    add(paste0(formatC(n_years, width = 3, flag = "0"),
+               " 01 01 31 12 Version ", substr(version, 1, 1)))
+
+    # Second line with scenario file (why did we copy it, if this line refers to the FOCUS dir...
+    add(paste0("..\\..\\..\\", sze_file_lower))
+
+    # Third line with psm file
+    add(basename(psm_file))
+
+    # Lines with climate files
+    for (year in 1:n_years) {
+      add(paste0("..\\..\\..\\", location_code, "_", formatC(year, width = 2, flag = "0"), ".cli"))
+    }
+
+    # Output control section
+    add("00000015 00.")
+    add("    PRSN    TSER  000.0     1.0 00000001")
+    add("    TETD    TSER  000.0     1.0 00000001")
+    add("    INFL    TSER  100.0     1.0 00000001")
+    add("    RUNF    TSER  000.0     1.0 00000001")
+    add("    THET    TSER  000.0     1.0 00000001")
+    add("    THET    TSER  030.0     1.0 00000001")
+    add("    TEMP    TSER  000.0     1.0 00000001")
+    add("    TEMP    TSER  030.0     1.0 00000001")
+    add("    TPAP    TSER  000.0  1.0E05 00000001")
+    add("    TDKF    TSER  000.0  1.0E05 00000001")
+    add("    TUPF    TSER  000.0  1.0E05 00000001")
+    add("    TPST    TSER  005.0  1.0E06 00000001")
+    add("    PFLX    TSER  100.0  1.0E05 00000001")
+    add("    RFLX    TSER  000.0  1.0E05 00000001")
+    add("    LEAC    TSER  100.0  1.0E09 00000001")
+
+    # Copy pelmo executable
+    file.copy(PELMO_exe, run_dir)
+
+    # In addition to the files copied by the FOCUS PELMO GUI, we
+    # need the error file in the run directory, as we start
+    # the exe file from this directory
+    file.copy(file.path(PELMO_base, "lf90.eer"), run_dir)
+  }
+
+  lapply(run_list, setup_run)
+
+  if (execute) {
+    run_PELMO(runs, version = version, PELMO_base = PELMO_base)
+  }
+}
+
+#' Run PELMO
+#'
+#' @inheritParams PELMO_runs
+#' @importFrom parallel mclapply
+#' @export
+run_PELMO <- function(runs, psm_dir = ".", version = "5.5.3", PELMO_base = "auto",
+                      cores = getOption("mc.cores", 2L))
+{
+
+  if (PELMO_base[1] == "auto") {
+    PELMO_base = file.path(system.file(package = "PELMO.installeR"),
+                                       paste0("FOCUSPELMO.", gsub("\\.", "", version)))
+  }
+
+  pelmo_exe = list.files(PELMO_base, "^pelmo.*.exe")
+
+  # Create list of runs to traverse using mclapply
+  run_list <- create_run_list(runs, check_psm_files = FALSE)
+
+  # In order to run in parallel, we need to make a directory tree starting
+  # from the base of the PELMO installation, as pelmo401.exe seems to look
+  # for ..\..\..\summary.PLM and does not start if this is in use by another
+  # instance
+  run_exe <- function(x) {
+    psm <- x$psm
+    crop <- x$crop
+    scenario <- x$scenario
+    run_dir <- file.path(PELMO_base, "FOCUS", PELMO_path(psm, crop, scenario))
+    exe_dir <- file.path(PELMO_base, paste(sample(c(0:9, letters, LETTERS),
+                                                  size = 8, replace = TRUE),
+                                           collapse = ""))
+    dir.create(exe_dir)
+
+    # Fresh PELMO base directory for every run with random name
+    run_dir_exe <- file.path(exe_dir, "FOCUS", PELMO_path(psm, crop, scenario))
+    dir.create(run_dir_exe, recursive = TRUE)
+
+    # Copy the contents of the run directory
+    run_files <- list.files(run_dir, full.names = TRUE)
+    file.copy(run_files, run_dir_exe)
+
+    # Copy FOCUS files
+    location_code <- FOCUS_PELMO_location_codes[x$scenario]
+    location_files <- list.files(file.path(PELMO_base, "FOCUS"), paste0("^", location_code),
+                                 full.names = TRUE)
+    focus_dir_exe <- file.path(exe_dir, "FOCUS")
+    file.copy(location_files, focus_dir_exe)
+    file.copy(file.path(PELMO_base, "FOCUS", "FOCUSCHECK.DAT"), focus_dir_exe)
+
+    # We need to go the directory to simplify calling pelmo with wine
+    setwd(run_dir_exe)
+    psm_file <- file.path(psm_dir, paste0(psm, ".psm"))
+    system(paste("wine", pelmo_exe, psm_file), ignore.stdout = TRUE)
+
+    # Copy the results to the original run directory
+    plm_files <- list.files(run_dir_exe, ".plm$", full.names = TRUE)
+    PLM_files <- list.files(run_dir_exe, ".PLM$", full.names = TRUE)
+    file.copy(c(plm_files, PLM_files), run_dir)
+
+    # Clean up
+    unlink(exe_dir, recursive = TRUE)
+  }
+
+  mclapply(run_list, run_exe, mc.cores = cores)
+}
+
+#' Create a path of run directories as the PELMO GUI does
+#'
+#' @export
+#' @importFrom utils data
+#' @param psm The psm identifier
+#' @param crop The PELMO crop acronym
+#' @param scenario The scenario
+PELMO_path <- function(psm, crop, scenario) {
+  if (crop %in% names(FOCUS_PELMO_crops)) {
+    crop <- FOCUS_PELMO_crops[crop]
+  }
+  if (!crop %in% FOCUS_PELMO_crops) stop(crop, " is not in FOCUS_PELMO_crops.")
+  psm_dir <- paste0(psm, ".run")
+  crop_dir <- paste0(gsub(" ", "_-_", crop), ".run")
+
+  # Deal with 'irrigation'. It only affects the naming of the scenario directory,
+  # but it does not appear to affect the PELMO run. This can be seen when
+  # comparing runs set up for Beans (field) and Beans (vegetable) for the Porto
+  # scenario
+  irrigation_string <- ""
+
+  # 'Irrigation' is only possible in the GUI for five scenarios
+  if (scenario %in% c("Cha", "Pia", "Por", "Sev", "Thi")) {
+    crop_acronyms <- names(FOCUS_PELMO_crops)
+    names(crop_acronyms) <- FOCUS_PELMO_crops
+    crop_acronym <- crop_acronyms[[crop]]
+    # Some crops are not 'irrigated' according to the GUI
+    if (crop_acronym %in% c("win", "fbe", "woi", "ape", "spr")) {
+      irrigation_string <- "--_-_no_-_irrigation"
+    } else {
+      irrigation_string <- "--_-_irrigated"
+    }
+  }
+
+  scenario_dir <- paste0(
+    FOCUS_GW_scenarios_2012$names[scenario], "_-_(", FOCUS_PELMO_location_codes[scenario], ")",
+    irrigation_string, ".run")
+
+  return(file.path(psm_dir, crop_dir, scenario_dir))
+}
+
+#' Create a list of runs that we can traverse
+#'
+#' @inheritParams PELMO_runs
+#' @param check_psm_files Should we check if the psm file exists
+create_run_list <- function(runs, psm_dir = ".", check_psm_files = FALSE) {
+  i <- 0
+  run_list <- list()
+  for (run in runs) {
+    psm <- run$psm
+    if (check_psm_files) {
+      psm_file <- file.path(psm_dir, paste0(psm, ".psm"))
+      if (file.access(psm_file) == -1) {
+        stop(psm_file, " is not readable")
+      }
+    }
+    crops <- setdiff(names(run), "psm")
+    for (crop in crops) {
+      crop_acronyms <- names(FOCUS_PELMO_crops)
+      names(crop_acronyms) <- FOCUS_PELMO_crops
+      if (!crop %in% crop_acronyms) {
+        if (crop %in% FOCUS_PELMO_crops) {
+          crop <- crop_acronyms[crop]
+        } else {
+          stop("Invalid crop specification ", crop)
+        }
+      }
+      for (scenario in run[[crop]]) {
+        i <- i + 1
+        run_list[[i]] <- list(psm = psm, crop = crop, scenario = scenario)
+      }
+    }
+  }
+  return(run_list)
+}