LakeEnsemblR: Run Ensemble of Lake Models

Documented in export_flow

#' Export LakeEnsemblR standardised flow files to model specific driver format
#'
#' Export in- anbd out-flow driver files for each model
#'
#' @param config_file filepath; to LakeEnsemblr yaml master config file
#' @param model vector; model to export driving data. Options include
#'    c("GOTM", "GLM", "Simstrat", "FLake", "MyLake")
#' @param folder filepath; to folder which contains the model folders generated by export_config()
#'
#' @examples
#' \dontrun{
#' export_flow(config_file, model = c("GOTM", "GLM", "Simstrat", "FLake", "MyLake"))
#' }
#' @importFrom gotmtools get_yaml_value calc_cc input_yaml
#' @importFrom glmtools read_nml set_nml write_nml
#' @importFrom configr read.config
#' @importFrom lubridate floor_date ceiling_date
#'
#' @export
export_flow <- function(config_file, model = c("GOTM", "GLM", "Simstrat", "FLake", "MyLake"),
                          folder = ".") {

  # It's advisable to set timezone to GMT in order to avoid errors when reading time
  original_tz  <-  Sys.getenv("TZ")
  Sys.setenv(TZ = "UTC")

  # Set working directory
  oldwd <- getwd()
  setwd(folder)

  # this way if the function exits for any reason, success or failure, these are reset:
  on.exit({
    setwd(oldwd)
    Sys.setenv(TZ = original_tz)
  })

  # check model input
  model <- check_models(model)

  ##-------------Read settings---------------
  # initial water level
  init_lvl <- get_yaml_value(config_file, "location", "init_depth")
  # surface elevation
  surf_lvl <- get_yaml_value(config_file, "location", "elevation")
  # bottom elevation
  bot_lvl <- surf_lvl - get_yaml_value(config_file, "location", "depth")

  # Get start & stop dates
  start_date <- get_yaml_value(config_file, "time", "start")
  stop_date <- get_yaml_value(config_file, "time", "stop")

  # Use inflows
  use_inflows <- get_yaml_value(config_file, "inflows", "use")
  
  # Use outflows
  tryCatch({get_yaml_value(config_file, "inflows", "mass-balance")
    warning(paste0("The 'mass-balance' argument is no longer used after ",
                   "version 1.1. If you would like to have outflows ",
                   "matching the inflows, please add the inflow file ",
                   "manually to the 'outflows' section. You can use the same ",
                   "or a different file as for inflows."))},
    error = function(e) { })
  
  # Fix water level - assumed to be FALSE if not present in config_file
  fix_wlvl <- tryCatch(get_yaml_value(config_file, "inflows", "fix_wlvl"),
                       error = function(e) {return(FALSE)})
  
  # read in yaml file
  yml_fl <- read.config(config_file)
  
  if("outflows" %in% names(yml_fl)) {
    use_outflows <- get_yaml_value(config_file, "outflows", "use")
  } else {
    use_outflows <- FALSE
  }
  

  if(use_outflows){
    # number of outflows
    num_outflows <- yml_fl$outflows$number_outflows
    # outflow depths
    lvl_outflows <- yml_fl$outflows$outflow_lvl
    # Get scaling parameter
    if(!is.null(yml_fl$scaling_factors$all$outflow)){
      scale_param_out <- get_yaml_value(file = config_file, "all", "outflow")
    }else{
      scale_param_out <- rep(1, num_outflows)
      warning(paste0("No outflow scaling found in section 'scaling_factors/all'",
                     " of the config file '", config_file,
                     "'. Inflow scaling was assumed to be 1.",
                     " If you provide model specific scaling they",
                     " will overwrite this factor."), call. = FALSE)
    }
    if(!is.null(yml_fl$outflows$scale_param)){
      warning(paste0("Outflow scaling found in section 'outflow/scale_param'",
                     " of the config file '", config_file,
                     "'. outflow scaling was moved to the ",
                     "'scaling_factors' section, please update",
                     " your yaml file."), call. = FALSE)
    }
  }

  if(use_inflows){
    # number of inflows
    num_inflows <- get_yaml_value(config_file, "inflows", "number_inflows")
    # inflow depths - always assumed to be at the surface
    lvl_inflows <- rep(-1, num_inflows)
    # Get scaling parameter
    if(!is.null(yml_fl$scaling_factors$all$inflow)) {
      scale_param_inf <- get_yaml_value(file = config_file, "all", "inflow")
    }else {
      scale_param_inf <- rep(1, num_inflows)
      warning(paste0("No inflow scaling found in section 'scaling_factors/all'",
                     " of the config file '", config_file,
                     "'. Inflow scaling was assumed to be 1.",
                     " If you provide model specific scaling they",
                     " will overwrite this factor."), call. = FALSE)
    }
    if(!is.null(yml_fl$inflows$scale_param)) {
      warning(paste0("Inflow scaling found in section 'inflow/scale_param'",
                     " of the config file '", config_file,
                   "'. inflow scaling was moved to the ",
                     "'scaling_factors' section, please update",
                     " your yaml file."), call. = FALSE)
    }
  }


  ##---------------FLake-------------

  if("FLake" %in% model){
    fla_fil <- file.path(folder, get_yaml_value(config_file, "config_files", "FLake"))

    if(!use_inflows){
      input_nml(fla_fil, label = "inflow", key = "Qfromfile",  ".false.")
    }else{
      input_nml(fla_fil, label = "inflow", key = "Qfromfile",  ".true.")
    }
  }

  ##---------------GLM-------------

  if("GLM" %in% model){
    glm_nml <- file.path(folder, get_yaml_value(config_file, "config_files", "GLM"))

    # Read in nml and input parameters
    nml <- read_nml(glm_nml)
    
    # 2022-02-03: template file in GLM3r, ggplot_overhaul branch does not have
    # entries flt_off_sw, crest_width, and crest_factor. And with outlet_type,
    # the model crashes. 
    if(is.null(nml[["outflow"]][["flt_off_sw"]])){
      nml[["outflow"]][["flt_off_sw"]] <- FALSE
    }
    if(is.null(nml[["outflow"]][["crest_width"]])){
      nml[["outflow"]][["crest_width"]] <- 100
    }
    if(is.null(nml[["outflow"]][["crest_factor"]])){
      nml[["outflow"]][["crest_factor"]] <- 0.61
    }
    if(!is.null(nml[["outflow"]][["outlet_type"]])){
      nml[["outflow"]][["outlet_type"]] <- NULL
    }
    
    # if no inflow or outflow is used  this list is keep otherwise it is changed
    inp_list <- list("num_inflows" = 0, "num_outlet" = 0)
    # set inflow
    if (use_inflows) {
      inp_list$num_inflows  <-  num_inflows
      inp_list <- c(inp_list, list("names_of_strms" = paste0("inflow_", 1:num_inflows),
                                   "strm_hf_angle" = rep(65, num_inflows),
                                   "strmbd_slope" = rep(2, num_inflows),
                                   "strmbd_drag" = rep(0.016, num_inflows),
                                   "inflow_factor" = rep(1, num_inflows),
                                   "inflow_fl" = paste0("inflow_", 1:num_inflows, ".csv")))
    }

    # set outflows
    if (use_outflows){
      outf_surf <- rep(FALSE, num_outflows)
      outf_surf[lvl_outflows == -1] <- TRUE
      #!! outflow elevations need to be in meters above sea level!!
      lvl_outflows_glm <- lvl_outflows + bot_lvl
      # outflow lvl for floating outflows is set to 0
      lvl_outflows_glm[lvl_outflows == -1] <- 0
      inp_list$num_outlet <- num_outflows
      inp_list <- c(inp_list, list("flt_off_sw" = outf_surf,
                                   "outl_elvs" = lvl_outflows_glm,
                                   "outflow_fl" = paste0("outflow_", 1:num_outflows, ".csv"),
                                   "outflow_factor" = rep(1, num_outflows)))
    }

    nml <- glmtools::set_nml(nml, arg_list = inp_list)
    write_nml(nml, glm_nml)

  }

  ##---------------GOTM-------------

  if("GOTM" %in% model) {
    got_yaml <- file.path(folder, get_yaml_value(config_file, "config_files", "GOTM"))
    yml_no_comment <- unname(sapply(readLines(got_yaml), function(x) strsplit(x, "#")[[1]][1]))
    # number of inflows in the yaml file so far
    num_inf_yaml <- length(grep("inflow\\_*\\d*:", yml_no_comment, value = TRUE))

    if(use_outflows) {
      # number of outflows in the yaml file so far
      num_outf_yaml <- length(grep("outflow\\_*\\d*:", yml_no_comment, value = TRUE))
    } else {
      num_outf_yaml <- 0
    }


    ## Switch off streams
    if(!use_inflows) {
      # switch of flexible water level
      input_yaml_multiple(got_yaml, key1 = "mimic_3d", key2 = "zeta", key3 = "method",
                          value = 0)
      input_yaml_multiple(got_yaml, key1 = "water_balance_method", value = 0)
      # remove all inflows but one
      if (num_inf_yaml > 1) {
        for (i in 2:(num_inf_yaml)) {
          rm_yaml_sec(got_yaml, paste0("inflow_", i))
          num_inf_yaml <- 1
        }
      }
      # streams_switch(file = got_yaml, method = "off")
      input_yaml_multiple(got_yaml, key1 = "streams", key2 = "inflow", key3 = "flow", key4 =
                            "method", value = 0)
      input_yaml_multiple(got_yaml, key1 = "streams", key2 = "inflow", key3 = "temp", key4 =
                            "method", value = 0)
      input_yaml_multiple(got_yaml, key1 = "streams", key2 = "inflow", key3 = "salt", key4 =
                            "method", value = 0)
    }
    if (!use_outflows) {
      # remove all outflows but one
      if (num_outf_yaml > 1) {
        for (i in 2:(num_outf_yaml)) {
          rm_yaml_sec(got_yaml, paste0("outflow_", i))
          num_outf_yaml <- 1
        }
      }
      input_yaml_multiple(got_yaml, key1 = "streams", key2 = "outflow", key3 = "flow", key4 =
                            "method", value = 0)
      input_yaml_multiple(got_yaml, key1 = "streams", key2 = "outflow", key3 = "temp", key4 =
                            "method", value = 0)
      input_yaml_multiple(got_yaml, key1 = "streams", key2 = "outflow", key3 = "salt", key4 =
                            "method", value = 0)
    }
    # set inflows
    if (use_inflows) {
      # switch on flexible water level
      input_yaml_multiple(got_yaml, key1 = "mimic_3d", key2 = "zeta", key3 = "method",
                          value = 3)
      input_yaml_multiple(got_yaml, key1 = "water_balance_method", value = 3)
      # remove additional inflows that are not needed
      if (num_inf_yaml > num_inflows) {
        for (i in 2:(num_inf_yaml)) {
          rm_yaml_sec(got_yaml, paste0("inflow_", i))
          num_inf_yaml <- 1
        }
      }
      # add additional inflows if necessary
      if(num_inflows > 1 & num_inflows != num_inf_yaml) {
        for (i in num_inflows:2) {
          doubl_yaml_sec(got_yaml, "inflow", paste0("_", i))
          num_inf_yaml <- num_inflows
        }
      }
        # set inflow settings for all inflows
      for (i in 1:num_inflows) {

        if(i == 1) {
          inf_sec <- "inflow"
        } else {
          inf_sec <- paste0("inflow_", i)
        }

        # streams_switch(file = got_yaml, method = "on")
        input_yaml_multiple(got_yaml, key1 = "streams", key2 = inf_sec, key3 = "flow", key4 =
                              "method", value = 2)
        input_yaml_multiple(got_yaml, key1 = "streams", key2 = inf_sec, key3 = "temp", key4 =
                              "method", value = 2)
        input_yaml_multiple(got_yaml, key1 = "streams", key2 = inf_sec, key3 = "salt", key4 =
                              "method", value = 2)
        input_yaml_multiple(got_yaml, key1 = "streams", key2 = inf_sec, key3 = "flow", key4 =
                              "file", value = paste0("inflow_file_", i, ".dat"))
        input_yaml_multiple(got_yaml, key1 = "streams", key2 = inf_sec, key3 = "temp", key4 =
                              "file", value = paste0("inflow_file_", i, ".dat"))
        input_yaml_multiple(got_yaml, key1 = "streams", key2 = inf_sec, key3 = "salt", key4 =
                              "file", value = paste0("inflow_file_", i, ".dat"))
      }
      
      # Fix wlvls if fix_wlvl is TRUE
      if(fix_wlvl){
        input_yaml_multiple(got_yaml, key1 = "mimic_3d", key2 = "zeta", key3 = "method",
                            value = 0)
        input_yaml_multiple(got_yaml, key1 = "water_balance_method", value = 0)
      }
      
    }

    # set outflows
    if (use_outflows) {
      if (!use_inflows & !fix_wlvl) {
        # switch on flexible water level
        input_yaml_multiple(got_yaml, key1 = "mimic_3d", key2 = "zeta", key3 = "method",
                            value = 3)
        input_yaml_multiple(got_yaml, key1 = "water_balance_method", value = 3)
      }
      # remove additional outflows that are not needed
      if (num_outf_yaml > num_outflows) {
        for (i in 2:(num_outf_yaml)) {
          rm_yaml_sec(got_yaml, paste0("outflow_", i))
          num_outf_yaml <- 1
        }
      }
      outf_surf <- rep(FALSE, num_outflows)
      outf_surf[lvl_outflows == -1] <- TRUE
      # outflow lvl in GOTM are meters below initial surface lvl
      lvl_outflows_gotm <- lvl_outflows - init_lvl
      # add additional outflows if necessary
      if(num_outflows > 1 & num_outflows != num_outf_yaml) {
        for (i in num_outflows:2) {
          doubl_yaml_sec(got_yaml, "outflow", paste0("_", i))
          num_outf_yaml <- num_outflows
        }
      }
        # set outflow settings for all outflows
        for (i in 1:num_outflows) {
          if(i == 1) {
            outf_sec <- "outflow"
          } else {
            outf_sec <- paste0("outflow_", i)
          }
          # streams_switch(file = got_yaml, method = "on")
          input_yaml_multiple(got_yaml, key1 = "streams", key2 = outf_sec, key3 = "method",
                              value = ifelse(outf_surf[i], 1, 3))
          input_yaml_multiple(got_yaml, key1 = "streams", key2 = outf_sec, key3 = "zl",
                              value = lvl_outflows_gotm[i] - 0.5)
          input_yaml_multiple(got_yaml, key1 = "streams", key2 = outf_sec, key3 = "zu",
                              value = lvl_outflows_gotm[i] + 0.5)
          input_yaml_multiple(got_yaml, key1 = "streams", key2 = outf_sec, key3 = "flow", key4 =
                                "method", value = 2)
          input_yaml_multiple(got_yaml, key1 = "streams", key2 = outf_sec, key3 = "temp", key4 =
                                "method", value = 0)
          input_yaml_multiple(got_yaml, key1 = "streams", key2 = outf_sec, key3 = "salt", key4 =
                                "method", value = 0)
          input_yaml_multiple(got_yaml, key1 = "streams", key2 = outf_sec, key3 = "flow", key4 =
                                "file", value = paste0("outflow_file_", i, ".dat"))
          input_yaml_multiple(got_yaml, key1 = "streams", key2 = outf_sec, key3 = "temp", key4 =
                                "file", value = paste0("outflow_file_", i, ".dat"))
          input_yaml_multiple(got_yaml, key1 = "streams", key2 = outf_sec, key3 = "salt", key4 =
                                "file", value = paste0("outflow_file_", i, ".dat"))
      }
    }
   }

  ##---------------Simstrat-------------

  if("Simstrat" %in% model){
    sim_par <- file.path(folder, get_yaml_value(config_file, "config_files", "Simstrat"))
    # Turn off inflow
    if(!use_inflows){
      ## Set InflowMode to 0
      input_json(sim_par, "ModelConfig", "InflowMode", 0)
      
      ## Set Qin, Qout, Tin, and Sin to 0
      # See format_flow_simstrat function in helpers_flow.R
      qin_file <- format_flow_simstrat(flow_file = NULL, levels = NULL,
                                       surf_flow = NULL,
                                       in_or_out = "inflow", type = "discharge",
                                       sim_par = sim_par)
      qout_file <- format_flow_simstrat(flow_file = NULL, levels = NULL,
                                        surf_flow = NULL,
                                        in_or_out = "outflow", type = "discharge",
                                        sim_par = sim_par)
      tin_file <- format_flow_simstrat(flow_file = NULL, levels = NULL,
                                       surf_flow = NULL,
                                       in_or_out = "inflow", type = "temp",
                                       sim_par = sim_par)
      sin_file <- format_flow_simstrat(flow_file = NULL, levels = NULL,
                                       surf_flow = NULL,
                                       in_or_out = "inflow", type = "salt",
                                       sim_par = sim_par)
      
      writeLines(qin_file, file.path(folder, "Simstrat/Qin.dat"))
      writeLines(qout_file, file.path(folder, "Simstrat/Qout.dat"))
      writeLines(tin_file, file.path(folder, "Simstrat/Tin.dat"))
      writeLines(sin_file, file.path(folder, "Simstrat/Sin.dat"))
    }
    if(!use_outflows){
      qout_file <- format_flow_simstrat(flow_file = NULL, levels = NULL,
                                        surf_flow = NULL,
                                        in_or_out = "outflow", type = "discharge",
                                        sim_par = sim_par)
      writeLines(qout_file, file.path(folder, "Simstrat/Qout.dat"))
    }
  }

  ##---------------MyLake-------------

  if("MyLake" %in% model){
    # Load config file MyLake
    load(get_yaml_value(config_file, "config_files", "MyLake"))

    if(!use_inflows){
      mylake_config[["Inflw"]] <- matrix(rep(0, 8 * length(seq.POSIXt(from = floor_date(as.POSIXct(start_date),
                                                                                        unit = "day"),
                                                                      to = ceiling_date(as.POSIXct(stop_date),
                                                                                        unit = "day"),
                                                                      by = "day"))),
                                         ncol = 8)

      # save lake-specific config file for MyLake
      temp_fil <- gsub(".*/", "", get_yaml_value(config_file, "config_files", "MyLake"))
      save(mylake_config, file = file.path(folder, "MyLake", temp_fil))
    }
  }

  ##-------------If inflow == TRUE---------------

  if(use_inflows == TRUE){
    inflow_file <- get_yaml_value(file = config_file, label = "inflows", key = "file")
    # Check if file exists
    if(!file.exists(inflow_file)){
      stop(inflow_file, " does not exist. Check filepath in ", config_file)
    }
    ### Import data
    message("Loading inflow data...")
    inflow <- read.csv(file.path(folder, inflow_file), stringsAsFactors = FALSE)
    inflow[, 1] <- as.POSIXct(inflow[, 1])
    start_date <- get_yaml_value(config_file, "time", "start")
    # Stop date
    stop_date <- get_yaml_value(config_file, "time", "stop")
    inflow_start <- which(inflow$datetime == as.POSIXct(start_date))
    inflow_stop <- which(inflow$datetime == as.POSIXct(stop_date))
    
    # In case the start and stop dates do not occur in the flow files
    if(length(inflow_start) == 0L){
      inflow_start <- max(1L, min(which(inflow$datetime >= as.POSIXct(start_date))) - 1L)
    }
    if(length(inflow_stop) == 0L){
      inflow_stop <- min(nrow(inflow), max(which(inflow$datetime <= as.POSIXct(stop_date))) + 1L)
    }
    
    inflow <- inflow[inflow_start:inflow_stop, ]

    ### Naming conventions standard input
    chk_names_flow(inflow, num_inflows, inflow_file)

    ##### FLake
    if("FLake" %in% model){
      # check if model specific scaling factors are present
      if(!is.null(yml_fl$scaling_factors$FLake$inflow)) {
        scale_param_tmp <- yml_fl$scaling_factors$FLake$inflow
      } else {
        scale_param_tmp <- scale_param_inf
      }
      inflow_tmp <-  scale_flow(inflow, num_inflows, scale_param_tmp)

      # if multiple inflows are present put them in a list
      if(num_inflows > 1) {
        inflow_ls <- list()
        for (i in 1:num_inflows) {
          inflow_ls[[paste0("inflow_", i)]] <-
            data.frame(datetime = inflow_tmp$datetime,
                       Flow_metersCubedPerSecond = inflow_tmp[[paste0("Flow_metersCubedPerSecond_", i)]],
                       Water_Temperature_celsius = inflow_tmp[[paste0("Water_Temperature_celsius_", i)]],
                       Salinity_practicalSalinityUnits =
                         inflow_tmp[[paste0("Salinity_practicalSalinityUnits_", i)]]
            )
        }
        inflow_tmp <- inflow_ls
        rm(inflow_ls)
      } else {
        inflow_tmp <- list(inflow_1 = inflow_tmp)
      }
      
      flake_inflow <- format_inflow(inflow = inflow_tmp, model = "FLake",
                                    config_file = config_file)

      flake_outfile <- "Tinflow"

      flake_outfpath <- file.path(folder, "FLake", flake_outfile)

      write.table(flake_inflow, flake_outfpath, quote = FALSE, row.names = FALSE,
                  sep = "\t",
                  col.names = FALSE)


      temp_fil <- get_yaml_value(config_file, "config_files", "FLake")
      input_nml(temp_fil, label = "inflow", key = "time_step_number", nrow(flake_inflow))

      message("FLake: Created file ", file.path(folder, "FLake", flake_outfile))


    }

    ###### GLM
    if("GLM" %in% model){
      if(!is.null(yml_fl$scaling_factors$GLM$inflow)) {
        scale_param_tmp <- yml_fl$scaling_factors$GLM$inflow
      } else {
        scale_param_tmp <- scale_param_inf
      }
      inflow_tmp <-  scale_flow(inflow, num_inflows, scale_param_tmp)
      
      # if multiple inflows are present put them in a list
      if(num_inflows > 1) {
        inflow_ls <- list()
        for (i in 1:num_inflows) {
          inflow_ls[[paste0("inflow_", i)]] <-
            data.frame(datetime = inflow_tmp$datetime,
                       Flow_metersCubedPerSecond = inflow_tmp[[paste0("Flow_metersCubedPerSecond_", i)]],
                       Water_Temperature_celsius = inflow_tmp[[paste0("Water_Temperature_celsius_", i)]],
                       Salinity_practicalSalinityUnits =
                         inflow_tmp[[paste0("Salinity_practicalSalinityUnits_", i)]]
            )
        }
        inflow_tmp <- inflow_ls
        rm(inflow_ls)
      } else {
        inflow_tmp <- list(inflow_1 = inflow_tmp)
      }
      
      for (i in 1:num_inflows) {
        glm_inflow <- format_inflow(inflow = inflow_tmp[[i]], model = "GLM",
                                    config_file = config_file)

        inflow_outfile <- file.path("GLM", paste0("inflow_", i, ".csv"))
        write.csv(glm_inflow, inflow_outfile, row.names = FALSE, quote = FALSE)
        message("GLM: Created file ", file.path(folder, "GLM",
                                                paste0("inflow_", i, ".csv")))
      }
    }

    ##### GOTM
    if("GOTM" %in% model){

      if(!is.null(yml_fl$scaling_factors$GOTM$inflow)) {
        scale_param_tmp <- yml_fl$scaling_factors$GOTM$inflow
      } else {
        scale_param_tmp <- scale_param_inf
      }
      inflow_tmp <-  scale_flow(inflow, num_inflows, scale_param_tmp)
      
      # if multiple inflows are present put them in a list
      if(num_inflows > 1) {
        inflow_ls <- list()
        for (i in 1:num_inflows) {
          inflow_ls[[paste0("inflow_", i)]] <-
            data.frame(datetime = inflow_tmp$datetime,
                       Flow_metersCubedPerSecond = inflow_tmp[[paste0("Flow_metersCubedPerSecond_", i)]],
                       Water_Temperature_celsius = inflow_tmp[[paste0("Water_Temperature_celsius_", i)]],
                       Salinity_practicalSalinityUnits =
                         inflow_tmp[[paste0("Salinity_practicalSalinityUnits_", i)]]
            )
        }
        inflow_tmp <- inflow_ls
        rm(inflow_ls)
      } else {
        inflow_tmp <- list(inflow_1 = inflow_tmp)
      }
      
      for (i in 1:num_inflows) {
        gotm_outfile <- paste0("inflow_file_", i, ".dat")
        gotm_outfpath <- file.path(folder, "GOTM", gotm_outfile)
        gotm_inflow <- format_inflow(inflow_tmp[[i]], model = "GOTM",
                                     config_file = config_file)

        write.table(gotm_inflow, gotm_outfpath, quote = FALSE, row.names = FALSE,
                    sep = "\t",
                    col.names = TRUE)

        message("GOTM: Created file ", file.path(folder, "GOTM", gotm_outfile))
      }

    }

    ## Simstrat
    if("Simstrat" %in% model){
      # Set inflow mode to manually-chosen depths
      input_json(sim_par, label = "ModelConfig", key = "InflowMode", value = 1)
      
      if(!is.null(yml_fl$scaling_factors$Simstrat$inflow)) {
        scale_param_tmp <- yml_fl$scaling_factors$Simstrat$inflow
      } else {
        scale_param_tmp <- scale_param_inf
      }
      inflow_tmp <-  scale_flow(inflow, num_inflows, scale_param_tmp)
      
      # if multiple inflows are present put them in a list
      if(num_inflows > 1) {
        inflow_ls <- list()
        for (i in 1:num_inflows) {
          inflow_ls[[paste0("inflow_", i)]] <-
            data.frame(datetime = inflow_tmp$datetime,
                       Flow_metersCubedPerSecond = inflow_tmp[[paste0("Flow_metersCubedPerSecond_", i)]],
                       Water_Temperature_celsius = inflow_tmp[[paste0("Water_Temperature_celsius_", i)]],
                       Salinity_practicalSalinityUnits =
                         inflow_tmp[[paste0("Salinity_practicalSalinityUnits_", i)]]
            )
        }
        inflow_tmp <- inflow_ls
        rm(inflow_ls)
      } else {
        inflow_tmp <- list(inflow_1 = inflow_tmp)
      }
      
      # output file paths
      inflow_outfpath <- file.path(folder, "Simstrat", "Qin.dat")
      temp_outfpath <- file.path(folder, "Simstrat", "Tin.dat")
      salt_outfpath <- file.path(folder, "Simstrat", "Sin.dat")

      sim_inflow <- format_inflow(inflow = inflow_tmp, model = "Simstrat",
                                  config_file = config_file)

      ## inflow files
      qin_file <- format_flow_simstrat(flow_file = sim_inflow,
                                       levels = lvl_inflows,
                                       surf_flow = (lvl_inflows == -1),
                                       in_or_out = "inflow",
                                       type = "discharge",
                                       sim_par = sim_par)
      writeLines(qin_file, inflow_outfpath)
      tin_file <-  format_flow_simstrat(flow_file = sim_inflow,
                                        levels = lvl_inflows,
                                        surf_flow = (lvl_inflows == -1),
                                        in_or_out = "inflow",
                                        type = "temp",
                                        sim_par = sim_par)
      writeLines(tin_file, temp_outfpath)
      sin_file <-  format_flow_simstrat(flow_file = sim_inflow,
                                        levels = lvl_inflows,
                                        surf_flow = (lvl_inflows == -1),
                                        in_or_out = "inflow",
                                        type = "salt",
                                        sim_par = sim_par)
      writeLines(sin_file, salt_outfpath)
      
      message("Simstrat: Created file ", inflow_outfpath)
    }

    ## MyLake
    if("MyLake" %in% model){

      if(!is.null(yml_fl$scaling_factors$MyLake$inflow)) {
        scale_param_tmp <- yml_fl$scaling_factors$MyLake$inflow
      } else {
        scale_param_tmp <- scale_param_inf
      }
      inflow_tmp <-  scale_flow(inflow, num_inflows, scale_param_tmp)
      
      # if multiple inflows are present put them in a list
      if(num_inflows > 1) {
        inflow_ls <- list()
        for (i in 1:num_inflows) {
          inflow_ls[[paste0("inflow_", i)]] <-
            data.frame(datetime = inflow_tmp$datetime,
                       Flow_metersCubedPerSecond = inflow_tmp[[paste0("Flow_metersCubedPerSecond_", i)]],
                       Water_Temperature_celsius = inflow_tmp[[paste0("Water_Temperature_celsius_", i)]],
                       Salinity_practicalSalinityUnits =
                         inflow_tmp[[paste0("Salinity_practicalSalinityUnits_", i)]]
            )
        }
        inflow_tmp <- inflow_ls
        rm(inflow_ls)
      } else {
        inflow_tmp <- list(inflow_1 = inflow_tmp)
      }
      
      temp_fil <- get_yaml_value(config_file, "config_files", "MyLake")
      load(temp_fil)

      mylake_inflow <- format_inflow(inflow = inflow_tmp, model = "MyLake",
                                     config_file = config_file)

      # discharge [m3/d], temperature [deg C], conc of passive tracer [-], conc of passive
      # sediment tracer [-], TP [mg/m3], DOP [mg/m3], Chla [mg/m3], DOC [mg/m3]
      dummy_inflow <- matrix(rep(1e-10, 8 *
                                   length(seq.POSIXt(from = floor_date(as.POSIXct(start_date), unit = "day"),
                                                     to = ceiling_date(as.POSIXct(stop_date), unit = "day"),
                                                     by = "day"))),
                             ncol = 8)
      dummy_inflow[, 1] <- mylake_inflow$Flow_metersCubedPerDay
      dummy_inflow[, 2] <- mylake_inflow$Water_Temperature_celsius
      dummy_inflow[, 5] <- dummy_inflow[, 5] * 1e7
      dummy_inflow[, 6] <- dummy_inflow[, 6] * 1e1


      mylake_config[["Inflw"]] <- dummy_inflow

      temp_fil <- gsub(".*/", "", temp_fil)
      # save lake-specific config file for MyLake
      save(mylake_config, file = file.path(folder, "MyLake", temp_fil))

      message("MyLake: Created file ", file.path(folder, "MyLake", temp_fil))

    }
  }

  ##-------------If outflow == TRUE---------------

  if(use_outflows == TRUE) {

    outflow_file <- get_yaml_value(file = config_file, label = "outflows", key = "file")
    # Check if file exists
    if(!file.exists(outflow_file)){
      stop(outflow_file, " does not exist. Check filepath in ", config_file)
    }

    ### Import data
    message("Loading outflow data...")
    outflow <- read.csv(file.path(folder, outflow_file), stringsAsFactors = FALSE)
    outflow[, 1] <- as.POSIXct(outflow[, 1])

    start_date <- get_yaml_value(config_file, "time", "start")
    # Stop date
    stop_date <- get_yaml_value(config_file, "time", "stop")

    outflow_start <- which(outflow$datetime == as.POSIXct(start_date))
    outflow_stop <- which(outflow$datetime == as.POSIXct(stop_date))
    
    # In case the start and stop dates do not occur in the flow files
    if(length(outflow_start) == 0L){
      outflow_start <- max(1L, min(which(outflow$datetime >= as.POSIXct(start_date))) - 1L)
    }
    if(length(outflow_stop) == 0L){
      outflow_stop <- min(nrow(outflow), max(which(outflow$datetime <= as.POSIXct(stop_date))) + 1L)
    }

    outflow <- outflow[outflow_start:outflow_stop, ]

    ### Naming conventions standard input
    chk_names_flow(outflow, num_outflows, outflow_file)
    
    # FLake
    ###
    if("FLake" %in% model) {
        message("FLake does not need outflows, as mass fluxes are not considered.")
    }

    # GLM
    #####
    if("GLM" %in% model) {
      if(!is.null(yml_fl$scaling_factors$GLM$outflow)) {
        scale_param_tmp <- yml_fl$scaling_factors$GLM$outflow
      } else {
        scale_param_tmp <- scale_param_out
      }
      outflow_tmp <-  scale_flow(outflow, num_outflows, scale_param_tmp)
      
      # if multiple outflows are present put them in a list
      if(num_outflows > 1) {
        outflow_ls <- list()
        for (i in 1:num_outflows) {
          outflow_ls[[paste0("outflow_", i)]] <-
            data.frame(datetime = outflow_tmp$datetime,
                       Flow_metersCubedPerSecond = outflow_tmp[[paste0("Flow_metersCubedPerSecond_", i)]])
        }
        outflow_tmp <- outflow_ls
        rm(outflow_ls)
      } else {
        outflow_tmp <- list(outflow_1 = outflow_tmp)
      }
      
      
      for (i in 1:num_outflows) {
        glm_outflow <- format_outflow(outflow = outflow_tmp[[i]], model = "GLM",
                                    config_file = config_file)

        outflow_outfile <- file.path("GLM", paste0("outflow_", i, ".csv"))
        write.csv(glm_outflow, outflow_outfile, row.names = FALSE, quote = FALSE)
        message("GLM: Created file ", file.path(folder, "GLM",
                                                paste0("outflow_", i, ".csv")))
      }
    }

    ## GOTM
    if("GOTM" %in% model) {

      if(!is.null(yml_fl$scaling_factors$GOTM$outflow)) {
        scale_param_tmp <- yml_fl$scaling_factors$GOTM$outflow
      } else {
        scale_param_tmp <- scale_param_out
      }
      outflow_tmp <-  scale_flow(outflow, num_outflows, scale_param_tmp)
      
      # if multiple outflows are present put them in a list
      if(num_outflows > 1) {
        outflow_ls <- list()
        for (i in 1:num_outflows) {
          outflow_ls[[paste0("outflow_", i)]] <-
            data.frame(datetime = outflow_tmp$datetime,
                       Flow_metersCubedPerSecond = outflow_tmp[[paste0("Flow_metersCubedPerSecond_", i)]])
        }
        outflow_tmp <- outflow_ls
        rm(outflow_ls)
      } else {
        outflow_tmp <- list(outflow_1 = outflow_tmp)
      }
      
      for (i in 1:num_outflows) {

        gotm_outfile <- paste0("outflow_file_", i, ".dat")

        gotm_outfpath <- file.path(folder, "GOTM", gotm_outfile)

        gotm_outflow <- format_outflow(outflow_tmp[[i]], model = "GOTM",
                                       config_file = config_file)
        write.table(gotm_outflow, gotm_outfpath, quote = FALSE, row.names = FALSE,
                    sep = "\t",
                    col.names = TRUE)

        message("GOTM: Created file ", file.path(folder, "GOTM", gotm_outfile))
      }
    }

    ## Simstrat
    if("Simstrat" %in% model) {
      # Set inflow mode to manually-chosen depths
      input_json(sim_par, label = "ModelConfig", key = "InflowMode", value = 1)
      
      
      if(!is.null(yml_fl$scaling_factors$Simstrat$outflow)) {
        scale_param_tmp <- yml_fl$scaling_factors$Simstrat$outflow
      } else {
        scale_param_tmp <- scale_param_out
      }
      outflow_tmp <-  scale_flow(outflow, num_outflows, scale_param_tmp)
      
      if(num_outflows > 1) {
        outflow_ls <- list()
        for (i in 1:num_outflows) {
          outflow_ls[[paste0("outflow_", i)]] <-
            data.frame(datetime = outflow_tmp$datetime,
                       Flow_metersCubedPerSecond = outflow_tmp[[paste0("Flow_metersCubedPerSecond_", i)]])
        }
        outflow_tmp <- outflow_ls
        rm(outflow_ls)
      } else {
        outflow_tmp <- list(outflow_1 = outflow_tmp)
      }
      
      outf_surf <- rep(FALSE, num_outflows)
      outf_surf[lvl_outflows == -1] <- TRUE
      #!! outflow elevations need to be relative to initial water level!!
      lvl_outflows_simstrat <- round(lvl_outflows - init_lvl)
      lvl_outflows_simstrat[outf_surf] <- -1

      # if there are several outflows with the same depth throw an error
      if(any(duplicated(lvl_outflows_simstrat))){
        stop(paste0("Duplicated outflow levels in simstrat detected, please merge outflows",
                    " at 1m intervals."))
      }
      
      outflow_outfpath <- file.path(folder, "Simstrat", "Qout.dat")
      sim_outflow <- format_outflow(outflow_tmp, "Simstrat", config_file, folder)
      
      qout_file <- format_flow_simstrat(flow_file = sim_outflow,
                                        levels = lvl_outflows_simstrat,
                                        surf_flow = outf_surf,
                                        in_or_out = "outflow",
                                        type = "discharge",
                                        sim_par = sim_par)
      writeLines(qout_file, outflow_outfpath)
      
      message("Simstrat: Created outflow file ", outflow_outfpath)
    }

    ## MyLake
    if("MyLake" %in% model) {
      message("MyLake does not need specific outflows, as it employs automatic overflow.")
    }
  }
  
  message("export_flow complete!")
}