R/.scenarios_parallel_all-substances.R

In KWB-R/flextreat.hydrus1d: R Package for Soil Water Balance and Solute Transport Modelling Scenarios for Project Flextreat

# Preface ----------------------------------------------------------------------
if (FALSE) {
  remotes::install_github("kwb-r/kwb.hydrus1d@dev")
  remotes::install_github("kwb-r/flextreat.hydrus1d@dev")
  remotes::install_github("kwb-r/kwb.db")
}

library(magrittr)
library(flextreat.hydrus1d)

# get_atm ----------------------------------------------------------------------
get_atm <- function(atm, extreme_rain = NULL)
{
  `%>%` <- magrittr::`%>%`

  if (is.null(extreme_rain)) {
    return(atm)
  }

  if (!extreme_rain %in% c("dry", "wet")) {
    stop("extreme_rain has to be either 'NULL', 'dry' or 'wet")
  }

  atm <- atm %>%
    dplyr::mutate(
      hydrologic_year = flextreat.hydrus1d::get_hydrologic_years(date),
      year = as.integer(format(date, format = "%Y")),
      day_of_year = lubridate::yday(date)
    )

  atm_stats <- atm %>%
    dplyr::group_by(year) %>%
    dplyr::summarise(
      rain_mm = sum(rain_mm, na.rm = TRUE),
      evapo_p_mean_mm = sum(evapo_p_mean_mm, na.rm = TRUE)
    )

  if (extreme_rain == "dry") {

    atm_dry <- atm_stats %>%
      dplyr::filter(rain_mm == min(rain_mm))

    atm_dry_sel <- atm %>%
      dplyr::filter(year == atm_dry$year) %>%
      dplyr::select(day_of_year, rain_mm)

    atm_dry_input <- atm %>%
      dplyr::select(- rain_mm, - hydrologic_year, - year) %>%
      dplyr::left_join(atm_dry_sel) %>%
      dplyr::mutate(rain_mm = dplyr::if_else(is.na(rain_mm), 0, rain_mm)) %>%
      dplyr::select(- day_of_year) %>%
      dplyr::relocate(rain_mm, .after = clearwater.mmPerDay)

    return(atm_dry_input)
  }

  if (extreme_rain == "wet") {
    atm_wet <- atm_stats %>%
      dplyr::filter(rain_mm == max(rain_mm))

    atm_wet_sel <- atm %>%
      dplyr::filter(year == atm_wet$year) %>%
      dplyr::select(day_of_year, rain_mm)

    atm_wet_input <- atm %>%
      dplyr::select(- rain_mm, - hydrologic_year, - year) %>%
      dplyr::left_join(atm_wet_sel) %>%
      dplyr::mutate(rain_mm = dplyr::if_else(is.na(rain_mm), 0, rain_mm)) %>%
      dplyr::select(- day_of_year) %>%
      dplyr::relocate(rain_mm, .after = clearwater.mmPerDay)

    return(atm_wet_input)
  }
}

# prepare_solute_input ---------------------------------------------------------
prepare_solute_input <- function(
    dat,
    selector,
    Ks = NULL,
    SnkL1 = NULL,
    diff_w = 0, diff_g = 0,
    kd = NULL,
    halftime_to_firstorderrate = NULL
)
{
  `%>%` <- magrittr::`%>%`

  stopifnot(nrow(dat) <= 10)

  if (selector$solute$No.Solutes != nrow(dat)) {
    selector$solute$No.Solutes <- nrow(dat)
  }

  solute_names <- sprintf("solute_%d", seq_len(nrow(dat)))

  solutes_new <- setNames(nm = solute_names, lapply(seq_len(nrow(dat)), function(i) {

    dat_sel <- dat[i,]

    ks <- kd(
      porosity = selector$waterflow$soil$ths - selector$waterflow$soil$thr,
      retardation = dat_sel$retard,
      bulk_density = selector$solute$transport$Bulk.d.
    )

    cols <- c(
      "Ks", "Nu", "Beta", "Henry", "SnkL1", "SnkS1", "SnkG1", "SnkL1'",
      "SnkS1'", "SnkG1'", "SnkL0",  "SnkS0", "SnkG0",  "Alfa"
    )

    reaction <- matrix(data = 0, ncol = length(cols), nrow = length(ks)) %>%
      as.data.frame() %>%
      tibble::as_tibble()

    names(reaction) <- cols
    reaction$Beta <- 1

    reaction$Ks <- if (is.null(Ks)) {
      round(ks, 2)
    } else {
      Ks
    }

    reaction$SnkL1 <- if (is.null(SnkL1)) {
      round(halftime_to_firstorderrate(dat_sel$half_life_days), 5)
    } else {
      SnkL1
    }

    list(
      diffusion = tibble::tibble(DifW = diff_w, DifG = diff_g),
      reaction = reaction
    )

  }))

  sel_tmp <- selector$solute[!names(selector$solute) %in% solute_names]

  solutes_new_list <- c(
    sel_tmp[1:which(names(sel_tmp) == "transport")],
    solutes_new,
    sel_tmp[which(names(sel_tmp) == "kTopSolute"):length(sel_tmp)]
  )

  c(
    selector[names(selector) != "solute"],
    list(solute = solutes_new_list)
  )
}

# get_mean ---------------------------------------------------------------------
get_mean <- function(col)
{
  x <- stringr::str_split_fixed(col, "-", n = 2)
  mode(x) <- "numeric"
  round(rowMeans(x), digits = 2)
}

# kd ---------------------------------------------------------------------------
kd <- function(porosity, retardation, bulk_density)
{
  #https://www3.epa.gov/ceampubl/learn2model/part-two/onsite/retard.html
  (retardation - 1) * porosity / bulk_density
}

# halftime_to_firstorderrate ---------------------------------------------------
halftime_to_firstorderrate <- function(half_time)
{
  if(half_time != 0) {
    #https://chem.libretexts.org/Courses/Bellarmine_University/BU%3A_Chem_104_(Christianson)/Phase_2%3A_Understanding_Chemical_Reactions/4%3A_Kinetics%3A_How_Fast_Reactions_Go/4.5%3A_First_Order_Reaction_Half-Life#mjx-eqn-21.4.2
    0.693 / half_time
  } else {
    0
  }
}

# provide_soil_columns ---------------------------------------------------------
provide_soil_columns <- function(path)
{
  `%>%` <- magrittr::`%>%`

  x <- kwb.db::hsGetTable(path, "my_results2", stringsAsFactors = FALSE) %>%
    janitor::clean_names() %>%
    dplyr::mutate(substanz_name = stringr::str_trim(substanz_name),
                  half_life_days = dplyr::case_when(
      grepl(">", hwz_tage) ~ hwz_tage %>%
        stringr::str_remove(">") %>%
        as.numeric() %>%
        round(digits = 2),
      grepl("<", hwz_tage) ~ hwz_tage %>%
        stringr::str_remove("<") %>%
        as.numeric() %>%
        round(digits = 2),
      grepl("-", hwz_tage) ~ get_mean(hwz_tage),
      .default = as.numeric(hwz_tage) %>%
        round(digits = 2)),
      class_id = dplyr::case_when(
        half_life_days < 20 ~ "I",
        half_life_days >= 20 & half_life_days < 200 ~ "II",
        half_life_days >= 200 ~ "III"),
      class_hlt = dplyr::case_when(
        half_life_days < 20 ~ "< 20 Tage",
        half_life_days >= 20 & half_life_days < 200 ~ "20 - 200 Tage",
        half_life_days >= 200 ~ ">= 200 Tage"),
      retard = dplyr::case_when(
        grepl("-", retardation) ~ get_mean(retardation),
        is.na(retardation) ~ 1L,
        .default = as.numeric(retardation) %>%
          round(digits = 2))) %>%
    dplyr::filter(!is.na(half_life_days)) %>%
    dplyr::mutate(id = 1:dplyr::n()) %>%
    dplyr::relocate(id) # %>% dplyr::mutate(retard = 1, half_life_days = 0)

  x %>%
    dplyr::left_join(x %>%
                       dplyr::count(class_id) %>%
    dplyr::rename(class_nsubstances = n)) %>%
    dplyr::mutate(class_label = sprintf("%s (%s, %2d Stoffe)",
                                        class_id,
                                        class_hlt,
                                        class_nsubstances))

}

# MAIN -------------------------------------------------------------------------

#path <- "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_2_Boden-Grundwasser/daten_karten/Sickerwasserprognose/column-studies/Stoffeigenschaften_Säulen.xlsx"
#path <- "Y:/WWT_Department/Projects/FlexTreat/Work-packages/AP3/3_1_4_Prognosemodell/StofflicheModellrandbedingungen.xlsx"
root_path <- "D:/hydrus1d/all-substances"
path <- file.path(root_path, "StofflicheModellrandbedingungen.xlsx")

soil_columns <- provide_soil_columns(path)

### Select 1 substance for 5 different half life classes defined in this table
selected_substances <- readr::read_csv("inst/extdata/input-data/substance_classes.csv")

soil_columns_selected <- soil_columns  %>%
  dplyr::filter(substanz_nr %in% selected_substances$substance_id) %>%
  dplyr::arrange(substanz_nr) %>%
  dplyr::mutate(id = 1:dplyr::n())

arg_combis <- kwb.utils::expandGrid(
  extreme_rain = c("", "wet", "dry"), # "" needs to be treated as NULL!
  treatment = c("tracer", "ka", "o3"), # "tracer" # c("ka", "o3")
  scenario = unlist(lapply(c(1,10), function(x) {
    paste0("soil-", 1:3, sprintf("m_irrig-%02ddays", x))
  })),
  irrig_only_growing_season = c(TRUE, FALSE),
  duration_string = "long", #c("short", "long"),
  retardation_scenario = c("retardation_yes", "retardation_no")
)


# generate_solute_ids ----------------------------------------------------------
generate_solute_ids <- function(n)
{
  seq_start <- seq(1, n, 10)

  seq_end <- if (n > 10) {
    seq(10, n, 10)
  } else {
    n
  }

  if (length(seq_end) < length(seq_start)) {
    seq_end <- c(seq_end, n)
  }

  tibble::tibble(
    start = as.integer(seq_start),
    end = as.integer(seq_end)
  )
}

# generate_periods -------------------------------------------------------------
generate_periods <- function(n)
{
  tibble::tibble(
    start = seq(1, n, 10),
    end = if (n %% 10 != 0) {
      c(seq(10, n, 10), n)
    } else {
      seq(10, n, 10)
    }
  )
}

# prepare_files_for_irrig_int --------------------------------------------------
prepare_files_for_irrig_int <- function(paths)
{
  `%>%` <- magrittr::`%>%`

  p <- kwb.utils::createAccessor(paths)

  copy <- function(fun, from, to) {
    kwb.utils::catAndRun(
      sprintf("Copying from\n  %s to\n  %s", from, to),
      fun(path = from, new_path = to)
    )
  }

  if (!fs::dir_exists(p("model_dir"))) {
    copy(fs::dir_copy, from = p("model_dir_org"), to = p("model_dir"))
  }

  if (!fs::file_exists(p("model_gui_path"))) {
    copy(fun = fs::file_copy, from = p("model_gui_path_org"), to = p("model_gui_path"))
  }

  model_out_files <- list.files(
    p("model_dir"),
    pattern = "\\.out$",
    full.names = TRUE
  )

  if (length(model_out_files) > 0L) {
    fs::file_delete(model_out_files)
  }

  soil_depth_cm <- 100 *
    stringr::str_extract(p("model_name"), "soil-[0-9]+?m") %>%
    stringr::str_extract("[0-9]") %>%
    as.numeric()

  if (soil_depth_cm != 200) {
    soil_profile <- kwb.hydrus1d::read_profile(paths$profile)
    profile_extended <- kwb.hydrus1d::extend_soil_profile(
      soil_profile$profile,
      x_end = -soil_depth_cm
    )
    soil_profile_extended <- soil_profile
    soil_profile_extended$profile <- profile_extended
    kwb.hydrus1d::write_profile(soil_profile_extended, path = p("profile"))
  }

  string_irrig_int <- stringr::str_extract(p("model_dir"), "[0-9][0-9]?days")

  # Return the string that is used as "irrig_interval"
  paste(
    as.integer(stringr::str_extract(string_irrig_int, "\\d+")),
    stringr::str_extract(string_irrig_int, "[a-z]+")
  )
}

# sum_per_interval -------------------------------------------------------------
sum_per_interval <- function(data, interval)
{
  `%>%` <- magrittr::`%>%`

  data_org <- data
  data <- dplyr::select(data, tidyselect::all_of(
    c("date", "groundwater.mmPerDay", "clearwater.mmPerDay")
  ))
  cols_sum <- setdiff(names(data), "date")
  data_summed <- data %>%
    dplyr::mutate(
      group = lubridate::floor_date(date, unit = interval)
    ) %>%  # Konvertiere date in datetime-Format
    dplyr::group_by(group) %>%  # Gruppiere nach Zeitintervallen
    dplyr::summarise_at(
      .vars = tidyselect::all_of(cols_sum),
      .funs = sum
    ) %>%   # Berechne die Summe für jedes Intervall
    dplyr::rename(date = group)
  data_org[, cols_sum] <- 0
  data_org[data_org$date %in% data_summed$date, cols_sum] <- data_summed[, cols_sum]
  data_org
}

# get_valid_exe_path -----------------------------------------------------------
get_valid_exe_path <- function(exe_dir)
{
  if (file.exists(file.path(exe_dir, "H1D_CALC.exe"))) {
    file.path(exe_dir, "H1D_CALC.exe")
  } else {
    fs::dir_create(exe_dir)
    kwb.hydrus1d::check_hydrus_exe(dir = exe_dir)
  }
}

# provide_paths ----------------------------------------------------------------
provide_paths <- function(config, start, end)
{
  #Y:\WWT_Department\Projects\FlexTreat\Work-packages\AP3\3_1_4_Prognosemodell\Hydrus1D\irrig_fixed\irrig-period_status-quo\long_dry\retardation_no
  tracer <- config$treatment == "tracer"
  # Define a path grammar
  PATH_GRAMMAR <- list(
    exe_dir = sprintf("D:/hydrus1d/all-substances/<irrig_dir_string>/<duration_string><extreme_rain_string>/%s", config$retardation_scenario),
    model_name_org = "model_to_copy",
    model_name = "<location>_<scenario>_soil-column_<solute_id_start><solute_id_end>",
    ###model_gui_path_org =  "<exe_dir>/<model_name_org>.h1d",
    model_gui_path_org =  "<model_dir_org>/<model_name_org>.h1d",

    model_gui_path = "<exe_dir>/<model_name>.h1d",
    modelvs_gui_path = "<exe_dir>/<model_name>_vs.h1d",

    model_dir_org = "D:/hydrus1d/<model_name_org>",

    model_dir = "<exe_dir>/<model_name>",
    model_dir_vs = "<exe_dir>/<model_name>_vs",
    atmosphere = "<model_dir>/ATMOSPH.IN",
    atmosphere_vs = "<model_dir_vs>/ATMOSPH.IN",
    a_level = "<model_dir>/A_LEVEL.out",
    hydrus1d = "<model_dir>/HYDRUS1D.DAT",
    selector = "<model_dir>/SELECTOR.IN",
    profile = "<model_dir>/PROFILE.DAT",
    obs_node = "<model_dir>/Obs_Node.out",
    balance = "<model_dir>/BALANCE.out",
    t_level = "<model_dir>/T_LEVEL.out",
    t_level_vs = "<model_dir_vs>/T_LEVEL.out",
    runinf = "<model_dir>/Run_Inf.out",
    solute_id = "1",
    solute = "<model_dir>/solute<solute_id>.out",
    solute_vs = "<model_dir_vs>/solute<solute_id>.out",
    soil_data = "<extdata>/input-data/soil/soil_geolog.csv",
    solute_id_start = sprintf("%02d", start),
    solute_id_end = sprintf("%02d", end),
    location = if (tracer) {
      "tracer"
    } else {
      sprintf("ablauf_%s_median", config$treatment)
    } #"ablauf_ka_median"
  )

  # Resolve the path grammar by replacing the placeholders recursively
  kwb.utils::resolve(
    PATH_GRAMMAR,
    irrig_dir_string = ifelse(
      config$irrig_only_growing_season,
      "irrig-period_growing-season",
      "irrig-period_status-quo"
    ),
    duration_string = config$duration_string,
    extreme_rain_string = if (config$extreme_rain != "") {
      paste0("_", config$extreme_rain)
    } else {
      ""
    },
    #final_subdir = ifelse(tracer, "tracer", config$retardation_scenario),
    scenario = config$scenario
  )
}

# inner_function ---------------------------------------------------------------
inner_function <- function(config, atm_data, soil_columns, helper)
{
  `%>%` <- magrittr::`%>%`

  {
    # Define constants
    IRRIGATION_COLUMNS <- c("groundwater.mmPerDay", "clearwater.mmPerDay")

    # Provide variables from config
    extreme_rain <- if (config$extreme_rain == "") {
      NULL
    } else {
      config$extreme_rain
    }
    tracer <- config$treatment == "tracer"
    no_retardation <- config$retardation == "retardation_no"

    if (no_retardation) {
      soil_columns$retard <- 1
    }

    atm <- helper("get_atm")(atm_data, extreme_rain)

    if (config$irrig_only_growing_season) {
      atm[which(!lubridate::month(atm$date) %in% 4:9), IRRIGATION_COLUMNS] <- 0
    }

    if (config$duration_string == "test") {
      atm <- dplyr::filter(atm, date >= "2017-05-01" & date <= "2018-04-30")
    } else if (config$duration_string == "short") {
      atm <- dplyr::filter(atm, date >= "2017-05-01" & date <= "2020-04-30")
    } else {
      atm <- dplyr::filter(atm, date >= "2017-05-01" & date <= "2023-12-31")
    }

    days_monthly <- lubridate::days_in_month(
      seq.Date(from = min(atm$date), to = max(atm$date), by = "month")
    )

    loop_df <- if (tracer) {
      helper("generate_periods")(n = length(days_monthly))
    } else {
      helper("generate_solute_ids")(n = nrow(soil_columns))
    }
  }

  kwb.utils::catAndRun(
    messageText = sprintf(
      "Running '%s' period for treatment '%s'",
      extreme_rain,
      config$treatment
    ),
    expr = {
      sapply(seq_len(nrow(loop_df)), function(i) {

        #i <- 1L
        {
          paths <- helper("provide_paths")(
            config,
            start = loop_df$start[i],
            end = loop_df$end[i]
          )

          solute_start_id <- as.numeric(paths$solute_id_start)
          solute_end_id <- as.numeric(paths$solute_id_end)
          n_solutes <- solute_end_id - (solute_start_id - 1)

          no_irrig <- stringr::str_detect(paths$model_dir, "no-irrig")
          irrig_int <- stringr::str_detect(paths$model_dir, "irrig-[0-9][0-9]?days")
        }

        if (irrig_int) {
          irrig_interval <- helper("prepare_files_for_irrig_int")(paths)
        }

        # no-irrigation
        if (no_irrig) {
          atm[, IRRIGATION_COLUMNS] <- 0
        }

        if (irrig_int) {
          atm <- helper("sum_per_interval")(data = atm, interval = irrig_interval)
        }

        days_total <- cumsum(days_monthly)

        indices <- seq_along(days_total)

        c_tops <- lapply(indices, function(i) {

          x <- rep(0, nrow(atm))
          x_min <- ifelse(i == 1, 1, days_total[i - 1] + 1)
          x[x_min:days_total[i]] <- rep(100, days_monthly[i])

          tib <- data.frame(x)

          colnames(tib) <- if (i == indices[1]) {
            "cTop"
          } else {
            sprintf("cTop%d", which(indices %in% i))
          }

          tib
        }) %>% dplyr::bind_cols()

        c_tops_sel <- c_tops[, paths$solute_id_start:paths$solute_id_end]

        atm_prep <- if (tracer) {
          flextreat.hydrus1d::prepare_atmosphere(
            atm = atm,
            conc_irrig_clearwater = c_tops_sel,
            conc_irrig_groundwater = c_tops_sel,
            conc_rain = c_tops_sel
          )
        } else {
          flextreat.hydrus1d::prepare_atmosphere(
            atm = atm,
            conc_irrig_clearwater = soil_columns[solute_start_id:solute_end_id, paths$location],
            conc_irrig_groundwater = 0,
            conc_rain = 0
          )
        }

        n_tsteps <- nrow(atm_prep)

        writeLines(
          kwb.hydrus1d::write_atmosphere(atm = atm_prep),
          paths$atmosphere
        )

        selector <- kwb.hydrus1d::read_selector(path = paths$selector)

        selector$time$tMax <- n_tsteps
        selector$time$MPL <- 250

        selector$time$TPrint <- seq(0, n_tsteps, n_tsteps/selector$time$MPL)

        if (selector$time$TPrint[1] == 0) {
          selector$time$TPrint[1] <- 1
        }

        solutes_new <- helper("prepare_solute_input")(
          dat = soil_columns[solute_start_id:solute_end_id,],
          Ks = if (tracer) 0, # else NULL
          SnkL1 = if (tracer) 0, # else NULL
          selector = selector,
          diff_w = 0,
          diff_g = 0,
          kd = helper("kd"),
          halftime_to_firstorderrate = helper("halftime_to_firstorderrate")
        )

        kwb.hydrus1d::write_selector(solutes_new, paths$selector)

        hydrus1d <- kwb.hydrus1d::read_hydrus1d(paths$hydrus1d)
        hydrus1d$Main$NumberOfSolutes <- n_solutes
        kwb.hydrus1d::write_hydrus1d(hydrus1d, paths$hydrus1d)

        exe_path <- helper("get_valid_exe_path")(paths$exe_dir)

        kwb.hydrus1d::run_model(
          exe_path = exe_path,
          model_path = paths$model_dir,
          print_output = FALSE
        )

        atmos <- kwb.hydrus1d::read_atmosph(paths$atmosphere)

        #atmos$data[names(c_tops)] <- c_tops

        atm_default <- atmos

        tlevel <- kwb.hydrus1d::read_tlevel(paths$t_level)

        vs_atm <- flextreat.hydrus1d::recalculate_ctop_with_virtualstorage(
          atm = atm_default$data,
          tlevel = tlevel,
          crit_v_top = - 0.05
        )

        atmos$data[names(vs_atm$df)] <- vs_atm$df

        fs::dir_copy(paths$model_dir, paths$model_dir_vs, overwrite = TRUE)
        fs::file_copy(paths$model_gui_path, paths$modelvs_gui_path, overwrite = TRUE)

        model_vs_out_files <- list.files(paths$model_dir_vs, pattern = "\\.out$", full.names = TRUE)

        if (length(model_vs_out_files) > 0) {
          fs::file_delete(model_vs_out_files)
        }

        writeLines(
          kwb.hydrus1d::write_atmosphere(atm = atmos$data),
          paths$atmosphere_vs
        )

        kwb.hydrus1d::run_model(
          exe_path = exe_path,
          model_path = paths$model_dir_vs,
          print_output = FALSE
        )

      })
    },
    dbg = TRUE
  )
}

# helper -----------------------------------------------------------------------
helper <- kwb.utils::createAccessor(list(
  get_atm = get_atm,
  generate_solute_ids = generate_solute_ids,
  provide_paths = provide_paths,
  prepare_files_for_irrig_int = prepare_files_for_irrig_int,
  sum_per_interval = sum_per_interval,
  prepare_solute_input = prepare_solute_input,
  halftime_to_firstorderrate = halftime_to_firstorderrate,
  get_valid_exe_path = get_valid_exe_path,
  generate_periods = generate_periods,
  kd = kd
))

# Main loop --------------------------------------------------------------------

if (FALSE) {

atm_data <- flextreat.hydrus1d::prepare_atmosphere_data()

arg_combis <- kwb.utils::expandGrid(
  extreme_rain = c("", "wet", "dry"), # "" needs to be treated as NULL!
  treatment = c("ka", "o3"), # "tracer" # c("ka", "o3")
  scenario = unlist(lapply(c(1,10), function(x) {
    paste0("soil-", 1:3, sprintf("m_irrig-%02ddays", x))
  })),
  irrig_only_growing_season = c(TRUE, FALSE),
  duration_string = "long", #c("short", "long"),
  retardation_scenario = c("retardation_yes", "retardation_no")
)

  configs <- lapply(seq_len(nrow(arg_combis)), function(i) {
    as.list(arg_combis[i, ])
  })

  system.time({
  # Sequential loop
  for (config in configs) {
    #config <- configs[[1L]]
    inner_function(
      config = config,
      atm_data = atm_data,
      soil_columns = soil_columns,
      helper = helper
    )
  }
  })


  system.time(expr = {
  # Parallel loop
  ncores <- parallel::detectCores() - 1
  #ncores <- 8
  cl <- parallel::makeCluster(ncores)


  parallel::parLapply(
    cl = cl,
    X = configs,
    fun = inner_function,
    atm_data = atm_data,
    soil_columns = soil_columns,
    helper = helper
  )

  parallel::stopCluster(cl)
  }
  )
}

# After main loop --------------------------------------------------------------
if (FALSE)
{
  solute <- kwb.hydrus1d::read_solute(paths$solute_vs) %>%
    dplyr::mutate(difftime = c(0,diff(time)))

  max(solute$sum_cv_top) + min(solute$sum_cv_bot) + min(solute$cv_ch1)

  plot(solute$time, solute$c_top)
  points(solute$c_bot, col = "red")
  (1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop)) * 100

  paths$solute_vs2 <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/H1D/1a2a_tracer_vs/solute2.out"

  solute <- kwb.hydrus1d::read_solute(paths$solute) %>%
    dplyr::mutate(difftime = c(0,diff(time)))

  (1 - max(solute$sum_cv_top)/sum(atmos$data$Prec*atmos$data$cTop2)) * 100

  sum(atmos$data$Prec)
  max(tlevel$sum_infil)
  max(tlevel$sum_evap)

  balance <- kwb.hydrus1d::read_balance(paths$balance)

  balance %>%
    dplyr::filter(
      subdomain_id == 0,
      time < 400,
      parameter == "CncBalT"
    ) %>%
    ggplot2::ggplot(
      mapping = ggplot2::aes(x = time, y = value, col = as.factor(solute_id))
    ) + ggplot2::geom_point()

  balance %>%
    dplyr::filter(
      subdomain_id == 0,
      time < 400,
      parameter == "CncBalR"
    ) %>%
    ggplot2::ggplot(mapping = ggplot2::aes(x = time, y = value, col = as.factor(solute_id))) +
    ggplot2::geom_point()

  sum(solute$difftime * as.numeric(solute$c_top) * as.numeric(tlevel$r_top))

  solute_day <- flextreat.hydrus1d::aggregate_solute(solute)

  plot(solute_day$date, solute_day$c_top)
  abline(h = 7)

  plot(atmos$data$tAtm, atmos$data$cTop)

  sum(solute_day$cv_top[solute_day$cv_top > 0])
  sum(atmos$data$Prec*atmos$data$cTop)
  sum(atmos$data$Prec*atmos$data$cTop)/sum(solute_day$cv_top[solute_day$cv_top > 0])
  sum(solute$cv_top)

  sum((as.numeric(solute$cv_top) * solute$difftime))
  sum((solute$cv_ch1 * c(0,diff(solute$time))))

  max(solute$sum_cv_top)

  condition <- solute$cv_top > 0
  sum(solute$cv_top[condition])
  condition <- solute$cv_top < 0
  sum(solute$cv_top[condition])

  solute_aggr_date <- flextreat.hydrus1d::aggregate_solute(solute)

  obsnode <- kwb.hydrus1d::read_obsnode(paths$obs_node)

  obsnode %>%
    dplyr::filter(variable == "flux") %>%
    dplyr::group_by(node_id) %>%
    dplyr::summarise(sum = sum(value))

  profile <- kwb.hydrus1d::read_profile(paths$profile)

  p <- obsnode %>%
    dplyr::left_join(profile[,c("node_id", "x")]) %>%
    dplyr::filter(variable == "conc1") %>%
    ggplot2::ggplot(mapping = ggplot2::aes(
      x = time,
      y = value,
      col = as.factor(x))
    ) +
    ggplot2::geom_point() +
    ggplot2::theme_bw()

  plotly::ggplotly(p)

  solute_aggr_date
  View(tlevel_aggr_date)
  View(solute_aggr_date)

  t_level <- kwb.hydrus1d::read_tlevel(paths$t_level)
  t_level

  ## t_level aggregate
  tlevel_aggr_date <- flextreat.hydrus1d::aggregate_tlevel(t_level)

  tlevel_aggr_yearmonth <- flextreat.hydrus1d::aggregate_tlevel(
    t_level,
    col_aggr = "yearmonth"
  )

  tlevel_aggr_year_hydrologic <- flextreat.hydrus1d::aggregate_tlevel(
    t_level,
    col_aggr = "year_hydrologic"
  ) %>%
    dplyr::filter(.data$diff_time >= 364) ### filter out as only may-october

  wb_date_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_date)
  wb_yearmonth_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_yearmonth)
  wb_yearhydrologic_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_year_hydrologic)

  wb_date_plot
  wb_yearmonth_plot
  wb_yearhydrologic_plot

  solute$time[min(which(solute$sum_cv_top == max(solute$sum_cv_top)))]

  paths$model_dir_vs

  #scenarios_solutes <- paste0(scenarios, "_soil-column")
  scenarios_solutes <- paste0("ablauf_", c("o3", "ka"), "_median")

  #root_path <- "D:/hydrus1d/irrig_fixed_01"
  #root_path <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed"
  root_path <- "D:/hydrus1d/all-substances"

  scenario_dirs <- fs::dir_ls(
    path =   root_path,
    recurse = TRUE,
    regexp = "retardation_no.*vs$",
    #regexp = "irrig-period_status-quo/long/retardation_no.*vs$",
    #regexp = "irrig-period_status-quo/long/retardation_no/ablauf_ka_median_soil-1m_irrig-01days_soil-column_0105.*vs$",
    type = "directory"
  )

  # Set up parallel plan
  system.time(expr = {
    # future::plan(future::multisession)

    future.apply::future_sapply(scenario_dirs, function(scenario_dir) {

      solutes_list <- setNames(lapply(scenarios_solutes, function(scenario) {
        solute_files <- fs::dir_ls(scenario_dir,
                                   regexp = "solute\\d\\d?.out",
                                   recurse = TRUE)

        stopifnot(length(solute_files) > 0)

        solutes <- setNames(lapply(solute_files, function(file) {
          kwb.hydrus1d::read_solute(file, dbg = TRUE)
        }), nm = solute_files) %>% dplyr::bind_rows(.id = "path")


        solute_files_df <- tibble::tibble(path = solute_files,
                                          model_solute_id = path  %>%  basename() %>% stringr::str_extract(pattern = "[0-9][0-9]?") %>% as.integer(),
                                          soilcolumn_id_start = path  %>%  dirname() %>% stringr::str_extract(pattern = "[0-9]{4}") %>% stringr::str_sub(1,2) %>% as.integer(),
                                          soilcolumn_id_end = path  %>%  dirname() %>% stringr::str_extract(pattern = "[0-9]{4}") %>% stringr::str_sub(3,4) %>% as.integer(),
                                          soil_column_id = soilcolumn_id_start + model_solute_id - 1) %>%
          dplyr::left_join(soil_columns, by = c(soil_column_id = "id"))


        dplyr::left_join(solutes, solute_files_df)
      }), nm = scenarios_solutes)

      solutes_df <- solutes_list %>%
        dplyr::bind_rows(.id = "scenario")

      solutes_df_stats <- solutes_df %>%
        dplyr::bind_rows(.id = "scenario") %>%
        dplyr::mutate(scen = stringr::str_remove(basename(dirname(path)), "_soil-column.*")) %>%
        dplyr::group_by(path, scen,substanz_nr, substanz_name) %>%
        dplyr::summarise(sum_cv_top = max(sum_cv_top),
                         sum_cv_bot = min(sum_cv_bot),
                         cv_ch1 = min(cv_ch1)) %>%
        dplyr::mutate(mass_balance_error_percent = 100*(sum_cv_top + cv_ch1 + sum_cv_bot)/sum_cv_top) %>%
        dplyr::arrange(mass_balance_error_percent)


      solutes_df_stats$soil <- solutes_df_stats$sum_cv_top + solutes_df_stats$sum_cv_bot + solutes_df_stats$cv_ch1

      saveRDS(solutes_df, file = file.path(scenario_dir, "solutes.rds"))

      openxlsx::write.xlsx(solutes_df_stats, file = file.path(scenario_dir, "hydrus_scenarios.xlsx"))
    })
    # Close the parallel plan
    future::plan(future::sequential)
  })

  retardation_short <- "retardation_yes"

  soil_columns_plot <- if(retardation_short == "retardation_yes") {
    soil_columns
  } else {
    soil_columns %>%
      dplyr::mutate(retard = 1)
  }

  scenario_dirs <- fs::dir_ls(
    path =   root_path,
    recurse = TRUE,
    #regexp = "irrig-period_growing-season/long/retardation_no/.*vs/hydrus_scenarios.xlsx$",
    regexp = sprintf("%s.*vs/hydrus_scenarios.xlsx$", retardation_short),
    type = "file"
  )

  #fs::file_info(scenario_dirs) %>% View()


  res_stats <- lapply(
    stats::setNames(nm = scenario_dirs),
    function(scenario_dir) {
      readxl::read_excel(scenario_dir)
    }
  )

  load_default_paths <- fs::dir_ls("D:/hydrus1d/all-substances/irrig-period_status-quo", recurse = TRUE, regexp = sprintf("long/%s/ablauf_ka_median_soil-2m_irrig-10days_.*vs/hydrus_scenarios.xlsx",  retardation_short))

 # load_default <- res_stats$`D:/hydrus1d/irrig_fixed_01/irrig-period_status-quo/long/retardation_no/ablauf_ka_median_soil-2m_irrig-10days_soil-column_0105_vs/hydrus_scenarios.xlsx`
  load_default <-  res_stats %>%
    dplyr::bind_rows(.id = "path") %>%
    dplyr::filter(path %in% load_default_paths) %>%
    # dplyr::mutate(retardation = basename(dirname(dirname(path))),
    #               duration = basename(dirname(dirname(dirname(path)))),
    #               irrigation_period = basename(dirname(dirname(dirname(dirname((path))))))) %>%
    dplyr::select(- path, - scen, - mass_balance_error_percent, - soil)

  names(load_default)[3:5] <- paste0("default_",   names(load_default)[3:5])

  load_default_sum <- as.data.frame(t(abs(colSums(load_default[3:5]))))


  res_stats_df <- dplyr::bind_rows(res_stats, .id = "path") %>%
    dplyr::mutate(retardation = basename(dirname(dirname(path))),
                  duration = basename(dirname(dirname(dirname(path)))),
                  irrigation_period = basename(dirname(dirname(dirname(dirname((path))))))) %>%
    dplyr::select(- path, - mass_balance_error_percent, - soil) %>%
    tidyr::separate(col = scen, into = c("treatment", "soil_depth_irrig"), sep = "n_s", remove = FALSE) %>%
    dplyr::mutate(treatment = sprintf("%sn", treatment),
                  soil_depth_irrig = sprintf("s%s", soil_depth_irrig)) %>%
    tidyr::separate(col = soil_depth_irrig, into = c("soil_depth", "irrigation_intervall"), sep = "_", remove = FALSE) %>%
    dplyr::mutate(duration_irrigperiod = paste0(duration, "_", irrigation_period))


  substances_per_class <- soil_columns %>%
    dplyr::count(class_id)

  res_stats_df_class <- res_stats_df %>%
    dplyr::left_join(soil_columns %>%  dplyr::select(substanz_name, class_id, class_label)) %>%
    dplyr::group_by(scen, treatment, soil_depth, irrigation_intervall, duration, irrigation_period, class_id, class_label) %>%
    dplyr::summarise(dplyr::across(c(sum_cv_top, sum_cv_bot), ~abs(sum(.x, na.rm = TRUE))),
                     number_of_substances = dplyr::n()) %>%
    dplyr::mutate(scenario_label = dplyr::case_when(
       scen == "ablauf_ka_median_soil-2m_irrig-10days" & duration == "long" & irrigation_period == "irrig-period_status-quo" ~ "Status Quo",
       scen == "ablauf_ka_median_soil-2m_irrig-10days" & duration == "long_wet" & irrigation_period == "irrig-period_status-quo" ~ "Klima, feucht",
       scen == "ablauf_ka_median_soil-2m_irrig-10days" & duration == "long_dry" & irrigation_period == "irrig-period_status-quo" ~ "Klima, trocken",
       scen == "ablauf_ka_median_soil-2m_irrig-10days" & duration == "long" & irrigation_period == "irrig-period_growing-season" ~ "Bewässerung (nur Mai-Sep)",
       scen == "ablauf_o3_median_soil-2m_irrig-10days" & duration == "long" & irrigation_period == "irrig-period_status-quo" ~ "Ozone-UV",
       scen == "ablauf_ka_median_soil-1m_irrig-10days" & duration == "long" & irrigation_period == "irrig-period_status-quo" ~ "Boden, 1 m",
       scen == "ablauf_ka_median_soil-3m_irrig-10days" & duration == "long" & irrigation_period == "irrig-period_status-quo" ~ "Boden, 3 m",
       scen == "ablauf_ka_median_soil-2m_irrig-01days" & duration == "long" & irrigation_period == "irrig-period_status-quo" ~ "Bewässerungsintervall, 1 Tag",
       .default = ""
    ))

#  View(res_stats_df_class)

  res_stats_df_class_selected <- res_stats_df_class %>%
    dplyr::filter(scenario_label != "") #%>%
    #dplyr::mutate(sum_cv_bot = 2.7e+07*10*sum_cv_bot/10e+9/10e+3,
    #              sum_cv_top = 2.7e+07*10*sum_cv_top/10e+9/10e+3)

  load_gw <- res_stats_df_class_selected %>%
    dplyr::group_by(scenario_label) %>%
    dplyr::summarise(total_sum_cv_bot = sum(sum_cv_bot))

  load_default_per_class <- res_stats_df_class_selected %>%
    dplyr::ungroup() %>%
    dplyr::filter(scenario_label == "Status Quo") %>%
    dplyr::group_by(class_label) %>%
    dplyr::rename(default_class_sum_cv_top = sum_cv_top,
                  default_class_sum_cv_bot = sum_cv_bot)  %>%
    dplyr::summarize(default_class_sum_cv_top = sum(default_class_sum_cv_top),
                     default_class_sum_cv_bot = sum(default_class_sum_cv_top))  %>%
    dplyr::select(class_label, default_class_sum_cv_top, default_class_sum_cv_bot)



  load_per_scenario_and_class <- res_stats_df_class %>%
    dplyr::group_by(scenario_label, class_label) %>%
    dplyr::rename(class_sum_cv_top = sum_cv_top,
                  class_sum_cv_bot = sum_cv_bot)  %>%
    dplyr::summarize(class_sum_cv_top = sum(class_sum_cv_top),
                     class_sum_cv_bot = sum(class_sum_cv_bot))  %>%
    dplyr::select(scenario_label, class_label, class_sum_cv_top, class_sum_cv_bot)



    load_per_scenario <- res_stats_df_class %>%
    dplyr::group_by(scenario_label) %>%
    dplyr::rename(scenario_sum_cv_top = sum_cv_top,
                  scenario_sum_cv_bot = sum_cv_bot)  %>%
    dplyr::summarize(scenario_sum_cv_top = sum(scenario_sum_cv_top),
                     scenario_sum_cv_bot = sum(scenario_sum_cv_bot))  %>%
    dplyr::select(scenario_label, scenario_sum_cv_top, scenario_sum_cv_bot)


    x_desired_order <- c("Status Quo",
                         "Boden, 1 m",
                         "Boden, 3 m",
                         "Klima, trocken",
                         "Klima, feucht",
                         "Bewässerungsintervall, 1 Tag",
                         "Bewässerung (nur Mai-Sep)",
                         "Ozone-UV"
                         )

    gg0 <- res_stats_df_class_selected %>%
      #dplyr::filter(scenario_label %in% c("Bewässerung (nur Mai-Sep)", "Status Quo", "Ozone-UV")) %>%
      # dplyr::group_by(scenario_label,  class_label) %>%
      # dplyr::summarise(total_sum_cv_bot = sum(sum_cv_bot),
      #                  total_sum_cv_top = sum(sum_cv_top)) %>%
      dplyr::bind_cols(load_default_sum) %>%
      dplyr::left_join( load_per_scenario_and_class) %>%
      dplyr::mutate(percent_sum_cv_bot =  class_sum_cv_bot / default_sum_cv_top,
                    percent_sum_cv_top =  class_sum_cv_top / default_sum_cv_top,
                    percent_gw = percent_sum_cv_bot/percent_sum_cv_top) %>%
    #dplyr::mutate(percent_sum_cv_top_remainder = percent_sum_cv_top - percent_sum_cv_bot) %>%
    tidyr::pivot_longer(cols = tidyselect::starts_with("percent_sum"),
                        names_to = "variable",
                        values_to = "value") %>%
    dplyr::mutate(variable = dplyr::case_when(
      variable == "percent_sum_cv_bot" ~ "Grundwasser",
      variable == "percent_sum_cv_top" ~ "Boden",
      .default = variable),
      perc_gw = dplyr::if_else(variable == "Grundwasser",
                               100*percent_gw,
                               100 - 100*percent_gw)#,
      # value = dplyr::if_else(variable == "Grundwasser",
      #                          - value,
      #                          value)
      ) %>%
    dplyr::left_join(load_gw) %>%
    dplyr::left_join(load_per_scenario) %>%
    ggplot2::ggplot(ggplot2::aes(fill = class_label,
                                 group = variable,
                                 y = value,
                                 x = forcats::fct_relevel(scenario_label, x_desired_order))) +
    ggplot2::labs(x = "Szenario",
                  y = "Gesamtfracht normiert auf den Status-Quo (%)",
                  fill = "Stoffeintrag",
                  title = ifelse(retardation_short == "retardation_no",
                                 "ohne Retardation",
                                 "mit Retardation")) +
    ggplot2::geom_bar(stat="identity") +
    ggplot2::facet_wrap(~ variable, scales = "fixed", ncol = 1) +
    #ggplot2::scale_fill_manual(values = reversed_palette) +
    ggplot2::scale_y_continuous(labels = scales::percent, breaks=seq(0, 1, 0.1), limits = c(0, 1.05)) +
    ggplot2::geom_text(ggplot2::aes(label = dplyr::if_else(round(value*100, 0) == 0,
                                                           "",
                                                           sprintf("%2.1f %%",  round(value*100, 1)))),
                       position = ggplot2::position_stack(vjust = 0.5),
                       size = 3, color = "black") +
    ggplot2::theme_bw() +
    ggplot2::theme(legend.position = "top")

    # Extrahiere die Standardfarben von ggplot2
    default_colors <- ggplot2::ggplot_build(gg0)$data[[1]]$fill

    # Kehre die Reihenfolge der Farben um
    reversed_colors <- rev(default_colors)

    gg0 <- gg0 +
    ggplot2::scale_fill_manual(values = reversed_colors)

  png(sprintf("stoffeintrag_%s.png", retardation_short), width = 1200, height = 800)
  gg0
  dev.off()


  res_stats_df_class_selected_agg <- res_stats_df_class_selected
    dplyr::group_by(scenario_label) %>%
    dplyr::summarise(total_sum_cv_bot = sum(sum_cv_bot)*10*irrigation$irrigation_area_sqm[1]/10e+9/10e3,
                     total_sum_cv_top = sum(sum_cv_top)*10*irrigation$irrigation_area_sqm[1]/10e+9/10e3,
                     tot = total_sum_cv_bot + total_sum_cv_top,
                     percent_gw = total_sum_cv_bot/total_sum_cv_top)



  gg1 <- res_stats_df_class_selected_agg %>%
    dplyr::mutate(total_sum_cv_top_remainder = total_sum_cv_top - total_sum_cv_bot) %>%
    dplyr::select(-total_sum_cv_top) %>%
  tidyr::pivot_longer(cols = tidyselect::starts_with("total"),
                      names_to = "variable",
                      values_to = "value") %>%
  dplyr::mutate(variable = dplyr::case_when(
    variable == "total_sum_cv_bot" ~ "Grundwasser",
    variable == "total_sum_cv_top_remainder" ~ "Boden",
    .default = variable),
    perc_gw = dplyr::if_else(variable == "Grundwasser",
                             100*percent_gw,
                             100 - 100*percent_gw)) %>%
    dplyr::left_join(load_gw) %>%
  ggplot2::ggplot(ggplot2::aes(fill = variable,
                  y = value,
                  x = forcats::fct_reorder(scenario_label, total_sum_cv_bot))) +
  ggplot2::labs(x = "Szenario",
                y = "Stoffeintrag (% im Vergleich zum Status QUo)",
                fill = "Stoffeintrag",
                title = ifelse(retardation_short == "retardation_no",
                               "ohne Retardation",
                               "mit Retardation")) +
  ggplot2::geom_bar(stat = "identity") +
  ggplot2::scale_y_continuous(labels = scales::percent) +
  ggplot2::geom_text(ggplot2::aes(label = sprintf("%d %%", round(perc_gw, 0))),
                     nudge_y = -0.05) +
  ggplot2::theme_bw() +
  ggplot2::theme(legend.position = "top")

  gg1

  load_status_quo <- res_stats_df_class_selected %>%
    dplyr::group_by(scenario_label) %>%
    dplyr::summarise(total_sum_cv_bot = sum(sum_cv_bot)) %>%
    dplyr::filter(scenario_label == "Status Quo") %>%
    dplyr::pull(total_sum_cv_bot) %>%  sum()

  share_gw_load <- res_stats_df_class_selected %>%
    dplyr::group_by(class_label) %>%
    dplyr::summarise(load_gw = sum(sum_cv_bot)) %>%
    dplyr::mutate(class_load_gw_percent = 100*load_gw/sum(load_gw))

  share_sw_load <- res_stats_df_class_selected %>%
    dplyr::group_by(class_label) %>%
    dplyr::summarise(load_sw = sum(sum_cv_top)) %>%
    dplyr::mutate(class_load_sw_percent = 100*load_sw/sum(load_sw))

  gg2 <- res_stats_df_class_selected %>%
    dplyr::left_join(share_gw_load) %>%
    dplyr::left_join(share_sw_load) %>%
    dplyr::mutate(class_label = stringr::str_remove(class_label, "\\)")) %>%
    dplyr::mutate(class_label = sprintf("%s, %3.1f %% Oberflächenstoffeintrag, %3.1f %% GW-Fracht)",
                                        class_label,
                                        class_load_sw_percent,
                                        class_load_gw_percent)) %>%
    dplyr::left_join(res_stats_df_class_selected_agg) %>%
    ggplot2::ggplot(ggplot2::aes(fill = class_label,
                                 y = sum_cv_bot/load_status_quo,
                                 x = forcats::fct_reorder(scenario_label, total_sum_cv_bot))) +
    ggplot2::geom_bar(position = "stack", stat = "identity") +
    # ggplot2::scale_y_log10(labels = scales::percent) +
    ggplot2::geom_text(ggplot2::aes(label = sprintf("%3.1f %%", round(100 * sum_cv_bot/load_status_quo, 1))),
                       position = ggplot2::position_stack(vjust = 0.5),
                       size = 3, color = "black") +
    ggplot2::labs(fill = "Halbwertszeitsklasse",
                  x = "Szenario",
                  y = "Stoffeintrag ins Grundwasser (% im Vergleich zum Status Quo)"#,
                  # title = ifelse(retardation_short == "retardation_no",
                  #                "ohne Retardation",
                  #                "mit Retardation")
                  ) +
    ggplot2::theme_bw() +
    ggplot2::theme(legend.position = "top")

  gg2

  gridExtra::grid.arrange(gg1, gg2, ncol = 1)

  plotly::ggplotly(gg)


  res_stats_df <- res_stats_df %>%
    dplyr::left_join(load_default, by = c("substanz_nr", "substanz_name")) %>%
    dplyr::mutate(percental_load_gw = dplyr::if_else(abs(default_sum_cv_bot) < 10000 | abs(sum_cv_bot) < 10000,
                                                     NA_real_,
                                                     100 + 100 * (abs(sum_cv_bot) - abs(default_sum_cv_bot)) /  abs(default_sum_cv_bot)))


  de <- TRUE


  if(de) {
  res_stats_df$duration_irrigperiod <- res_stats_df$duration_irrigperiod %>%
      stringr::str_replace("long_wet", "DWD, trocken (2013: 380 mm/a)") %>%
      stringr::str_replace("long_dry", "DWD, nass    (2023: 835 mm/a)") %>%
      stringr::str_replace("long", "DWD   (2017-2023: 611 mm/a)") %>%
      stringr::str_replace("_irrig-period_status-quo", ", Bewässerung: ganzjährig") %>%
      stringr::str_replace("_irrig-period_growing-season", ", Bewässerung: nur Vegetationsperiode (Apr-Sep)")

    res_stats_df$soil_depth_irrig <- res_stats_df$soil_depth_irrig %>%
      stringr::str_replace("soil-", "Boden ") %>%
      stringr::str_replace("_irrig-10days", ", Bewässerung: zehntägig") %>%
      stringr::str_replace("_irrig-01days", ", Bewässerung: täglich")

    ci_scenarios <- unique(res_stats_df$duration_irrigperiod)[order(unique(res_stats_df$duration_irrigperiod))]
    si_scenarios <- unique(res_stats_df$soil_depth_irrig)[order(unique(res_stats_df$soil_depth_irrig))]

    res_stats_df <- res_stats_df %>%
      dplyr::left_join(tibble::tibble(soil_depth_irrig = si_scenarios,
                   soil_depth_irrig_scenarios = sprintf("Bo%01dBi%02d",
                                                        stringr::str_extract(si_scenarios, "[0-9]") %>% as.integer(),
                                                        ifelse(stringr::str_detect(si_scenarios, "zehn"),
                                                               10, 1))
                   )) %>%
      dplyr::left_join(tibble::tibble(duration_irrigperiod = ci_scenarios,
                                      duration_irrigperiod_scenarios = sprintf("N%03dB%03d",
                                                                           stringr::str_extract(ci_scenarios, " ([0-9]+) ") %>% as.integer(),
                                                                           ifelse(stringr::str_detect(ci_scenarios, "ganz"),
                                                                                  405,
                                                                                  289)))
      )

  lang_id <- "de"
  lab_x <- "Bodenm\u00E4chtigkeit (m) und Bew\u00E4sserungsintervall (Tage)"
  lab_y <- "Prozentuale Fracht ins Grundwasser im Vergleich zu Status Quo (%)"
  title <- "Stoff: %s"
  subtitle <- "Halbwertszeit: %3.1f Tage, Retardation: %2.1f"
  legend_col <- "Niederschlag und Bew\u00E4sserung (mm/Jahr)"
  legend_shape <- "Technische Aufbereitung"
  treatment_wwtp <- "Konventionell (Median:"
  treatment_o3 <-   "UV-Ozon       (Median:"
  treatment_unit <- "ng/l"


  } else {
    lang_id <- "en"
    lab_x <- "Scenario"
    lab_y <- "Percental Load to Groundwater\ncompared to Status Quo (%)"
    title <- "Substance: %s"
    subtitle <- "half life time: %3.1f days, retardation: %2.1f"
    legend_col <- "Scenario"
    legend_shape <- "Treatment"
    treatment_wwtp <- "WWTP effluent (median: "
    treatment_o3 <- "Ozone effluent (median: "
    treatment_unit <- "ng/l"

    si_scenarios <- unique(res_stats_df$soil_depth_irrig)[order(unique(res_stats_df$soil_depth_irrig))]

    res_stats_df <- res_stats_df %>%
      dplyr::left_join(tibble::tibble(soil_depth_irrig = si_scenarios,
                                      soil_depth_irrig_scenarios = sprintf("S%01dI%01d",
                                                                           stringr::str_extract(si_scenarios, "[0-9]") %>% as.integer(),
                                                                           ifelse(stringr::str_detect(si_scenarios, "ten"),
                                                                                  0, 1)))
      )
  }

  #View(res_stats_df)


  substances <- unique(res_stats_df$substanz_name)[order(unique(res_stats_df$substanz_name))]

  substances2 <- soil_columns_plot %>%
    dplyr::select(substanz_name, half_life_days, retard) %>%
    dplyr::filter(substanz_name %in% substances) %>%
    dplyr::arrange(half_life_days, retard) %>%
    dplyr::pull(substanz_name)

  substance <- substances[1]

  jpeg(filename = "cardensartan.jpeg", height = 500, width = 1200, quality = 100)
  substance_meta  <- soil_columns_plot[soil_columns_plot$substanz_name == substance,]

  res_stats_df_sel <-  res_stats_df %>%
    dplyr::filter(substanz_name == substance)

  res_stats_df_sel$treatment_conc <- ""

  res_stats_df_sel$treatment_conc[which(res_stats_df_sel$treatment == "ablauf_ka_median")] <- sprintf("%s %5d %s)", treatment_wwtp, round(substance_meta$ablauf_ka_median, 0),   treatment_unit)
  res_stats_df_sel$treatment_conc[which(res_stats_df_sel$treatment == "ablauf_o3_median")] <- sprintf("%s %5d %s)", treatment_o3, round(substance_meta$ablauf_o3_median, 0),   treatment_unit)

  gg <- res_stats_df_sel %>%
    ggplot2::ggplot(mapping = ggplot2::aes(x = soil_depth_irrig_scenarios,
                                           y = percental_load_gw,
                                           col = duration_irrigperiod_scenarios,
                                           shape = treatment_conc)) +
    ggplot2::geom_jitter(width = 0.15, size = 3, alpha = 0.5) +
    ggplot2::scale_y_continuous(limits = c(0, 200), labels = scales::percent_format(scale = 1)) +
    ggplot2::labs(y = lab_y, x = lab_x,
                  title = sprintf(title, substance_meta$substanz_name),
                  subtitle = sprintf(subtitle,
                                     substance_meta$half_life_days,
                                     substance_meta$retard),
                  col = legend_col,
                  shape = legend_shape) +
    ggplot2::theme_bw() +
    ggplot2::theme(legend.position=c(0.85,0.85),
                   legend.box.just = "left",
                   legend.direction = "vertical",
                   legend.margin = ggplot2::margin(),
                   legend.text = ggplot2::element_text(size = 14),# face = "bold"),  # Größer und dick
                   legend.title = ggplot2::element_text(size = 15, face = "bold"),
                   axis.text.x = ggplot2::element_text(size = 14, face = "bold"),
                   axis.text.y = ggplot2::element_text(size = 14, face = "bold"),  # Größer und dick
                   axis.title.x = ggplot2::element_text(size = 15, face = "bold"),  # Größer und dick
                   axis.title.y = ggplot2::element_text(size = 15, face = "bold"),
                   plot.title = ggplot2::element_text(size = 17, face = "bold"),
                   plot.subtitle = ggplot2::element_text(size = 15) #,   # Größer und dick
                   #plot.subtitle = ggplot2::element_text(size = 14, face = "bold")
    ) +
    ggplot2::guides(color = ggplot2::guide_legend(nrow = 2, byrow = FALSE))

  print(gg)
  dev.off()


  pdff <- sprintf("%02d_percental-load-to-groundwater_per-substance_%s_%s.pdf",
                  length(substances),
                  stringr::str_replace(retardation_short, "_", "-"),
                  lang_id)
  #kwb.utils::preparePdf(pdff, borderHeight.cm = 0,  borderWidth.cm = 0, width.cm = 32.67, height.cm = 23.1)
  kwb.utils::preparePdf(pdff)


  sapply(substances2, function(substance) {

   substance_meta  <- soil_columns_plot[soil_columns_plot$substanz_name == substance,]

    res_stats_df_sel <-  res_stats_df %>%
      dplyr::filter(substanz_name == substance)

    res_stats_df_sel$treatment_conc <- ""

    res_stats_df_sel$treatment_conc[which(res_stats_df_sel$treatment == "ablauf_ka_median")] <- sprintf("%s %5d %s)", treatment_wwtp, round(substance_meta$ablauf_ka_median, 0),   treatment_unit)
    res_stats_df_sel$treatment_conc[which(res_stats_df_sel$treatment == "ablauf_o3_median")] <- sprintf("%s %5d %s)", treatment_o3, round(substance_meta$ablauf_o3_median, 0),   treatment_unit)

    gg <- res_stats_df_sel %>%
      ggplot2::ggplot(mapping = ggplot2::aes(x = soil_depth_irrig_scenarios,
                                             y = percental_load_gw,
                                             col = duration_irrigperiod_scenarios,
                                             shape = treatment_conc)) +
      ggplot2::geom_jitter(width = 0.15, size = 3, alpha = 0.5) +
      ggplot2::scale_y_continuous(limits = c(0, 200), labels = scales::percent_format(scale = 1)) +
      ggplot2::labs(y = lab_y, x = lab_x,
                    title = sprintf(title, substance_meta$substanz_name),
                    subtitle = sprintf(subtitle,
                                       substance_meta$half_life_days,
                                       substance_meta$retard),
                    col = legend_col,
                    shape = legend_shape) +
      ggplot2::theme_bw() +
      ggplot2::theme(legend.position=c(0.80,0.85),
                     legend.box.just = "left",
                     legend.direction = "vertical",
                     legend.margin = ggplot2::margin(),
                     legend.text = ggplot2::element_text(size = 12),# face = "bold"),  # Größer und dick
                     legend.title = ggplot2::element_text(size = 14, face = "bold"),
                     axis.text.x = ggplot2::element_text(size = 12, face = "bold"),
                     axis.text.y = ggplot2::element_text(size = 12, face = "bold"),  # Größer und dick
                     axis.title.x = ggplot2::element_text(size = 14, face = "bold"),  # Größer und dick
                     axis.title.y = ggplot2::element_text(size = 14, face = "bold"),
                     plot.title = ggplot2::element_text(size = 16, face = "bold")#,   # Größer und dick
                     #plot.subtitle = ggplot2::element_text(size = 14, face = "bold")
                     ) +
      ggplot2::guides(color = ggplot2::guide_legend(nrow = 2, byrow = FALSE))

  print(gg)
})

kwb.utils::finishAndShowPdf(pdff)


pdff <- sprintf("percental-load-to-groundwater_per-scenario_%s.pdf",
                sprintf(retardation_short, "_", "-"))

kwb.utils::preparePdf(pdff)

sapply(duration_irrigperiods, function(duration_irrigper) {

  res_stats_df_sel <-  res_stats_df %>%
    dplyr::filter(duration_irrigperiod == duration_irrigper)

  gg <- res_stats_df_sel %>%
    ggplot2::ggplot(mapping = ggplot2::aes(x = soil_depth_irrig,
                                           y = percental_load_gw,
                                           col = substanz_name,
                                           shape = treatment
                                           #col = treatment
    )) +
    ggplot2::geom_point(size = 3) +
    ggplot2::labs(x = "Scenario",
                  y = "Percental Load to Groundwater compared to Status Quo (%)",
                  title = sprintf("Scenario: %s %s (%s)",
                                  res_stats_df_sel$irrigation_period[1],
                                  res_stats_df_sel$duration[1],
                                  res_stats_df_sel$retardation[1]),
                  col = "Substance Name") +
    #ggplot2::scale_y_log10(limits = c(0,200)) +
    ggplot2::ylim(c(0,200)) +
    ggplot2::theme_bw() +
    ggplot2::theme(legend.position = "top",
                   axis.text.x =  ggplot2::element_text(angle = 90, vjust = 0.5, hjust=1))
  print(gg)
})

kwb.utils::finishAndShowPdf(pdff)

 #root_path <- "D:/hydrus1d/irrig_fixed_01"
 root_path <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed"

  model_paths <- fs::dir_ls(
    path =   root_path,
    recurse = TRUE,
    regexp = "tracer$",
    type = "directory"
  )

  #model_paths <- model_paths[stringr::str_detect(model_paths, "wet|dry", negate = TRUE)]

  scenarios <- sapply(c(1,10), function(x) paste0("soil-", 1:3, sprintf("m_irrig-%02ddays", x))) %>%
    as.vector()

  #unique(configs$scenario)

  #traveltimes_list$`D:/hydrus1d/irrig_fixed_01/irrig-period_status-quo/long/tracer`[!sapply(traveltimes_list$`D:/hydrus1d/irrig_fixed_01/irrig-period_status-quo/long/tracer`, function(x) any(class(x) == "try-error"))]

  ### read traveltimes sequential
  system.time(
  traveltimes_list <- lapply(model_paths, function(model_path) {
    setNames(nm = (scenarios), lapply(scenarios, function(scenario) {
      try({
        message(sprintf("Scenario: %s", scenario))
        solute_files <- fs::dir_ls(
          path = model_path,
          recurse = TRUE,
          regexp = sprintf("tracer_%s_.*vs/solute\\d\\d?.out", scenario)
        )
        flextreat.hydrus1d::get_traveltimes(solute_files, dbg = TRUE)
      })
    }))
  }))



  ### read traveltimes in parallel
  library(future.apply)

  # Set up parallel plan
  system.time(expr = {
  future::plan(future::multisession)

  traveltimes_list <- future.apply::future_lapply(model_paths, function(model_path) {
    setNames(nm = (scenarios), future.apply::future_lapply(scenarios, function(scenario) {
      try({
        message(sprintf("Scenario: %s", scenario))
        solute_files <- fs::dir_ls(
          path = model_path,
          recurse = TRUE,
          regexp = sprintf("tracer_%s_.*vs/solute\\d\\d?.out", scenario)
        )
        flextreat.hydrus1d::get_traveltimes(solute_files, dbg = TRUE)
      })
    }))
  })
  }
  )

  # Close the parallel plan
  future::plan(future::sequential)
  traveltimes_list_backup <-  traveltimes_list

  sapply(seq_along(traveltimes_list), function(i) {

    label <- sprintf("%s", extrahiere_letzte_drei_teile(names(traveltimes_list)[i]))

    htmlwidgets::saveWidget(flextreat.hydrus1d::plot_traveltimes(traveltimes_list[[i]] %>% dplyr::bind_rows(),
                                                                 title = label ,
                                                                 ylim = c(0,650)),
                            file = sprintf("traveltimes_%s.html", label))
  })


  # extrahiere_letzte_drei_teile -------------------------------------------------
  extrahiere_letzte_drei_teile <- function(pfad)
  {
    sapply(pfad, function(pf) {
      # Teile den Pfad anhand des Schrägstrichs auf
      teile <- unlist(strsplit(pf, "/"))

      # Waehle die letzten drei Teile aus
      letzte_drei_teile <- tail(teile, 3)

      paste0(letzte_drei_teile, collapse = "_")
    })
  }

  # Plotting ---------------------------------------------------------------------
  if (FALSE)
  {
    pdff <- "traveltimes_per-scenario.pdf"
    kwb.utils::preparePdf(pdff)
    # sapply(names(traveltimes_list), function(path) {

      # traveltimes_sel <- traveltimes_list[[path]] %>% dplyr::bind_rows()
      #
      # label <-  extrahiere_letzte_drei_teile(path)

      # traveltime_bp <- traveltimes_sel %>%
      #       dplyr::bind_rows() %>%
      #       dplyr::filter(percentiles == 0.5) %>%
      #   dplyr::bind_rows(.id = "scenario") %>%
      #   dplyr::filter(!stringr::str_detect(scenario, "1.5")) %>%
      #   dplyr::mutate(quarter = lubridate::quarter(date) %>% as.factor(),
      #                 soil_depth =  stringr::str_extract(scenario, "soil-.*m") %>%
      #                   stringr::str_remove_all("soil-|m") %>%  as.factor())

      traveltime_bp <- lapply(traveltimes_list, function(x) {
        x %>%
          dplyr::bind_rows() %>%
          dplyr::filter(percentiles == 0.5)
      }) %>% dplyr::bind_rows(.id = "scenario") %>%
        dplyr::filter(!stringr::str_detect(scenario, "1.5")) %>%
        dplyr::mutate(quarter = lubridate::quarter(date) %>% as.factor(),
                      soil_depth =  stringr::str_extract(scenario, "soil-.*m") %>%
                        stringr::str_remove_all("soil-|m") %>%  as.factor())


      scenario_by_median_traveltime <- traveltime_bp %>%
        dplyr::group_by(scenario) %>%
        dplyr::summarise(median = median(time_diff, na.rm = TRUE)) %>%
        dplyr::arrange(median)

      traveltime_bp <- traveltime_bp %>%
        dplyr::left_join(scenario_by_median_traveltime)

      y_lim <- c(0,350)


      tt_bp_total <- traveltime_bp %>%
        dplyr::mutate(scenario_short = extrahiere_letzte_drei_teile(scenario)) %>%
        ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_short, median), y = time_diff)) +
        ggplot2::geom_boxplot(outliers = FALSE) +
        ggplot2::geom_jitter(position = ggplot2::position_jitter(width = 0.1),
                             col = "darkgrey",
                             alpha = 0.6) +
        ggplot2::ylim(y_lim) +
        ggplot2::labs(y = "Median Traveltime (days)", x = "Scenario",
                      title = "Boxplot: median traveltime total") +
        ggplot2::theme_bw() +
        ggplot2::theme(axis.text.x =  ggplot2::element_text(angle = 90, vjust = 0.5, hjust=1))

      print(tt_bp_total)


    kwb.utils::finishAndShowPdf(pdff)

    pdff <- "traveltimes_all.pdf"
    kwb.utils::preparePdf(pdff)

    traveltimes_df <- lapply(traveltimes_list, function(sublist) sublist %>% dplyr::bind_rows()) %>%
      dplyr::bind_rows(.id = "scenario_main_raw") %>%
      dplyr::mutate(
        scenario_name = stringr::str_remove_all(model_name, "_soil-column_.*vs$") %>% stringr::str_remove_all("tracer_"),
        scenario_main =  scenario_main_raw %>% extrahiere_letzte_drei_teile(),
        quarter = lubridate::quarter(date) %>% as.factor(),
        soil_depth =  stringr::str_extract(scenario_name, "soil-.*m") %>%
          stringr::str_remove_all("soil-|m") %>%  as.factor())

    scenario_by_median_traveltime <- traveltimes_df %>%
      dplyr::group_by(scenario_main, scenario_name) %>%
      dplyr::summarise(median = median(time_diff, na.rm = TRUE)) %>%
      dplyr::arrange(median)

    traveltimes_bp <- traveltimes_df %>%
      dplyr::left_join(scenario_by_median_traveltime)

    y_lim <- c(0,350)

    tt_bp_total <- traveltimes_bp %>%
      ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_name, median),
                                   y = time_diff,
                                   col = scenario_main)) +
      ggplot2::geom_boxplot(outliers = FALSE) +
      # ggplot2::geom_jitter(position = ggplot2::position_jitterdodge(
      #   jitter.width = 0.1,
      #   dodge.width = 0.75),
      #   alpha = 0.6) +
      ggplot2::ylim(y_lim) +
      ggplot2::labs(y = "Median Traveltime (days)",
                    x = "Scenario",
                    title = "Boxplot: median traveltime total") +
      ggplot2::theme_bw() +
      ggplot2::theme(legend.position = "top")

    print(tt_bp_total)

    kwb.utils::finishAndShowPdf(pdff)

    htmlwidgets::saveWidget(
      widget = plotly::ggplotly(tt_bp_total),
      file = "traveltimes_all.html"
    )

    jpeg(filename = "verweilzeiten.jpeg", height = 500, width = 1200, quality = 100)
    #pdff <- "traveltimes_all_percent_de.pdf"
    #kwb.utils::preparePdf(pdff)

    traveltimes_df <- lapply(traveltimes_list, function(sublist) sublist %>% dplyr::bind_rows()) %>%
      dplyr::bind_rows(.id = "scenario_main_raw") %>%
      dplyr::mutate(
        scenario_name = stringr::str_remove_all(model_name, "_soil-column_.*vs$") %>% stringr::str_remove_all("tracer_"),
        scenario_main =  scenario_main_raw %>% extrahiere_letzte_drei_teile() %>% stringr::str_remove("_tracer") %>%
          stringr::str_remove("irrig-period_") %>% stringr::str_remove("_long"),
        quarter = lubridate::quarter(date) %>% as.factor(),
        soil_depth =  stringr::str_extract(scenario_name, "soil-.*m") %>%
          stringr::str_remove_all("soil-|m") %>%  as.factor(),
        irrigation_interval = stringr::str_extract(scenario_name, pattern = "irrig.*") %>% stringr::str_extract("[0-9]+"),
        label_x = sprintf("Bo%sBi%s", soil_depth, irrigation_interval)
        ) %>%
      tidyr::separate(col = scenario_main, into = c("irrigation_period", "climate"), sep = "_", remove = FALSE) %>%
      dplyr::mutate(climate = dplyr::case_when(climate == "wet" ~ "N835",
                                               climate == "dry" ~ "N380",
                                               .default = "N606"),
                    irrigation_period = dplyr::if_else(irrigation_period == "growing-season", "B289", "B405"),
                    label_legend = sprintf("%s%s", climate, irrigation_period))


    scenario_base_median <- traveltimes_df %>%
      dplyr::filter(
        scenario_name == "soil-2m_irrig-10days",
        scenario_main == "status-quo",
        percentiles == 0.5
      ) %>%
      dplyr::select(- time_top, - time_bot) %>%
      dplyr::rename(time_diff_base = time_diff)

    traveltimes_bp <- traveltimes_df %>%
      dplyr::filter(percentiles == 0.5) %>%
      dplyr::left_join(scenario_base_median[, c("month_id", "time_diff_base")] %>% dplyr::mutate(percentiles = 0.5)) %>%
      dplyr::mutate(time_diff_percent = 100 + 100 * (time_diff - time_diff_base) / time_diff_base) %>%
      dplyr::filter(!is.na(time_diff_percent))

  # Plotting ---------------------------------------------------------------------

  y_lim <- c(0,350)

  tt_bp_percent <- traveltimes_bp %>%
    ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(as.factor(label_x), time_diff_percent),
                                 y = time_diff_percent,
                                 col = label_legend)) +
    ggplot2::geom_boxplot(outliers = FALSE) +
    # ggplot2::geom_jitter(position = ggplot2::position_jitterdodge(
    #   jitter.width = 0.1,
    #   dodge.width = 0.75),
    #   alpha = 0.6) +
    ggplot2::ylim(y_lim) +
    ggplot2::labs(y = "Mittlere Verweilzeit (%) im Vergleich zu Status Quo",
                  x = "Bodenmächtigkeit (m) und Bewässerungsintervall (Tage)",
                  col = "Niederschlag und Bew\u00E4sserung (mm/Jahr)",
                  title = "") +
    ggplot2::theme_bw() +
    ggplot2::theme(#legend.position=c(0.85,0.85),
                   legend.box.just = "left",
                   legend.direction = "vertical",
                   legend.position = "top",
                   legend.margin = ggplot2::margin(),
                   legend.text = ggplot2::element_text(size = 14),# face = "bold"),  # Größer und dick
                   legend.title = ggplot2::element_text(size = 15, face = "bold"),
                   axis.text.x = ggplot2::element_text(size = 14, face = "bold"),
                   axis.text.y = ggplot2::element_text(size = 14, face = "bold"),  # Größer und dick
                   axis.title.x = ggplot2::element_text(size = 15, face = "bold"),  # Größer und dick
                   axis.title.y = ggplot2::element_text(size = 15, face = "bold"),
                   plot.title = ggplot2::element_text(size = 17, face = "bold"),
                   plot.subtitle = ggplot2::element_text(size = 15) #,   # Größer und dick
                   #plot.subtitle = ggplot2::element_text(size = 14, face = "bold")
    ) +
    ggplot2::guides(color = ggplot2::guide_legend(nrow = 1, byrow = FALSE))

  print(tt_bp_percent)


  #kwb.utils::finishAndShowPdf(pdff)
  dev.off()

  pdff <- "traveltimes_all_percent_en.pdf"
  kwb.utils::preparePdf(pdff)

  traveltimes_df <- lapply(traveltimes_list, function(sublist) sublist %>% dplyr::bind_rows()) %>%
    dplyr::bind_rows(.id = "scenario_main_raw") %>%
    dplyr::mutate(
      scenario_name = stringr::str_remove_all(model_name, "_soil-column_.*vs$") %>% stringr::str_remove_all("tracer_"),
      scenario_main =  scenario_main_raw %>% extrahiere_letzte_drei_teile(),
      quarter = lubridate::quarter(date) %>% as.factor(),
      soil_depth =  stringr::str_extract(scenario_name, "soil-.*m") %>%
        stringr::str_remove_all("soil-|m") %>%  as.factor())

  scenario_base_median <- traveltimes_df %>%
    dplyr::filter(
      scenario_name == "soil-2m_irrig-10days",
      scenario_main == "irrig-period_status-quo_long_tracer",
      percentiles == 0.5
    ) %>%
    dplyr::select(- time_top, - time_bot) %>%
    dplyr::rename(time_diff_base = time_diff)

  traveltimes_bp <- traveltimes_df %>%
    dplyr::filter(percentiles == 0.5) %>%
    dplyr::left_join(scenario_base_median[, c("month_id", "time_diff_base")] %>% dplyr::mutate(percentiles = 0.5)) %>%
    dplyr::mutate(time_diff_percent = 100 + 100 * (time_diff - time_diff_base) / time_diff_base)

  # Plotting ---------------------------------------------------------------------

  y_lim <- c(0,350)

  tt_bp_percent <- traveltimes_bp %>%
    ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario_name, time_diff_percent),
                                 y = time_diff,
                                 col = scenario_main)) +
    ggplot2::geom_boxplot(outliers = FALSE) +
    # ggplot2::geom_jitter(position = ggplot2::position_jitterdodge(
    #   jitter.width = 0.1,
    #   dodge.width = 0.75),
    #   alpha = 0.6) +
    ggplot2::ylim(y_lim) +
    ggplot2::labs(y = "Median Traveltime (%) compared to Status Quo",
                  x = "Scenario",
                  title = "Boxplot: median traveltime percent") +
    ggplot2::theme_bw() +
    ggplot2::theme(legend.position = "top")

  print(tt_bp_percent)

  kwb.utils::finishAndShowPdf(pdff)

  htmlwidgets::saveWidget(
    widget = plotly::ggplotly(tt_bp_total),
    file = "traveltimes_all.html"
  )

  tt_bp_total_soil <- traveltime_bp %>%
    ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario, median), y = time_diff, col = soil_depth)) +
    ggplot2::geom_boxplot(outliers = FALSE) +
    ggplot2::geom_jitter(position = ggplot2::position_jitter(width = 0.1),
                         alpha = 0.6) +
    ggplot2::ylim(y_lim) +
    ggplot2::labs(y = "Median Traveltime (days)", x = "Scenario",
                  col = "Soil Depth (m)",
                  title = "Boxplot: median traveltime total") +
    ggplot2::theme_bw() +
    ggplot2::theme(legend.position = "top")

  tt_bp_total_soil



  tt_bp_total_quartal <- traveltime_bp %>%
    ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario, median), y = time_diff)) +
    ggplot2::geom_boxplot(outliers = FALSE) +
    ggplot2::geom_jitter(position = ggplot2::position_jitter(width = 0.1),
                         mapping = ggplot2::aes(col = quarter),
                         alpha = 0.6) +
    ggplot2::ylim(y_lim) +
    ggplot2::labs(y = "Median Traveltime (days)", x = "Scenario",
                  col = "Quartal",
                  title = "Boxplot: median traveltime total") +
    ggplot2::theme_bw() +
    ggplot2::theme(legend.position = "top")

  tt_bp_total_quartal


  tt_bp_quarter <- traveltime_bp %>%
    ggplot2::ggplot(ggplot2::aes(x = forcats::fct_reorder(scenario, median), y = time_diff, col = quarter)) +
    ggplot2::geom_boxplot(outliers = FALSE) +
    ggplot2::geom_jitter(position = ggplot2::position_jitterdodge(
      jitter.width = 0.1,
      dodge.width = 0.75),
      alpha = 0.6) +
    ggplot2::labs(y = "Median Traveltime (days)",
                  x = "Scenario",
                  col = "Quartal",
                  title = "Boxplot: median traveltime by quarter") +
    ggplot2::ylim(y_lim) +
    ggplot2::theme_bw() +
    ggplot2::theme(legend.position = "top")

  tt_bp_quarter

  # Save HTML widgets ------------------------------------------------------------
  htmlwidgets::saveWidget(
    widget = plotly::ggplotly(tt_bp_total),
    title = "Boxplot: median traveltime total",
    file = "boxplot_traveltimes-median_total.html"
  )

  htmlwidgets::saveWidget(
    plotly::ggplotly(tt_bp_quarter),
    title = "Boxplot: median traveltime by quarter",
    file = "boxplot_traveltimes-median_quarter.html"
  )

  }

### check model

  model_dir <- "C:/kwb/projects/flextreat/3_1_4_Prognosemodell/Vivian/Rohdaten/irrig_fixed"

  atm_files <- fs::dir_ls(model_dir, recurse = TRUE, type = "file", regexp = "tracer.*ATMOSPH.IN")

  atm_df <- lapply(atm_files, function(file) {
    atm <- kwb.hydrus1d::read_atmosph(file)
    tibble::tibble(path = file,
                   Prec_sum = sum(atm$data$Prec, na.rm = TRUE))
  }) %>% dplyr::bind_rows()

  path_split <- stringr::str_split_fixed(atm_df$path,"/", 13)

  atm_df %>%
    dplyr::mutate(irrig_period = path_split[,9],
                  duration_extreme = path_split[,10],
                  scenario= path_split[,12] %>% stringr::str_remove("_soil-column_.*")) %>%
    dplyr::group_by(irrig_period, duration_extreme, scenario) %>%
    dplyr::summarise(Prec_mean = sum(Prec_sum)/dplyr::n()) %>%
    View()


}