Refactor BADM_IC_process() for Robust Single/Multi-site Handling

base/utils/tests/testthat/test-ud_convert.R

@@ -36,4 +36,21 @@ test_that("ud_convert() warns with wrong input units for difftime", {
   expect_warning(ud_convert(as.difftime("12:00:00"), u1 = "years", u2 = "minutes"))
   #should still error if units are not convertible
   expect_error(ud_convert(as.difftime("12:00:00"), u1 = "kilograms", u2 = "minutes"))
+})
+
+test_that("model-specific pool conversions", {
+  # DALEC/SIPNET C pools
+  expect_equal(ud_convert(100, "g/m2", "kg/m2"), 0.1)
+  # GDAY pools
+  expect_equal(ud_convert(10, "Mg/ha", "kg/m2"), 1)
+})
+
+test_that("model-specific flux conversions", {
+  # DALEC/SIPNET C fluxes
+  expect_equal(ud_convert(86400, "g/m2/d", "kg/m2/s"), 0.001, tolerance = 1e-10)
+})
+
+test_that("photosynthesis parameters", {
+  # Photosynthesis energy parameters
+  expect_equal(ud_convert(1000, "J/mol", "kJ/mol"), 1)
 })

models/dalec/R/model2netcdf.DALEC.R

@@ -83,22 +83,22 @@ model2netcdf.DALEC <- function(outdir, sitelat, sitelon, start_date, end_date) {
 
     ## Setup outputs for netCDF file in appropriate units
     output <- list()
-    ## Fluxes
-    output[[1]] <- (sub.DALEC.output[, 1] * 0.001)/timestep.s  # Autotrophic Respiration in kgC/m2/s
-    output[[2]] <- (sub.DALEC.output[, 21] + sub.DALEC.output[, 23]) * 0.001 / timestep.s  # Heterotrophic Resp kgC/m2/s
-    output[[3]] <- (sub.DALEC.output[, 31] * 0.001)/timestep.s  # GPP in kgC/m2/s    
-    output[[4]] <- (sub.DALEC.output[, 33] * 0.001)/timestep.s  # NEE in kgC/m2/s
-    output[[5]] <- (sub.DALEC.output[, 3] + sub.DALEC.output[, 5] + sub.DALEC.output[, 7]) * 0.001/timestep.s  # NPP kgC/m2/s
-    output[[6]] <- (sub.DALEC.output[, 9] * 0.001) / timestep.s  # Leaf Litter Flux, kgC/m2/s
-    output[[7]] <- (sub.DALEC.output[, 11] * 0.001) / timestep.s  # Woody Litter Flux, kgC/m2/s
-    output[[8]] <- (sub.DALEC.output[, 13] * 0.001) / timestep.s  # Root Litter Flux, kgC/m2/s
+    ## Fluxes (convert g/m2/day to kg/m2/s)
+    output[[1]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 1], "g/m2/d", "kg/m2/s") # Autotrophic Respiration
+    output[[2]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 21] + sub.DALEC.output[, 23], "g/m2/d", "kg/m2/s") # Heterotrophic Resp
+    output[[3]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 31], "g/m2/d", "kg/m2/s") # GPP
+    output[[4]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 33], "g/m2/d", "kg/m2/s") # NEE
+    output[[5]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 3] + sub.DALEC.output[, 5] + sub.DALEC.output[, 7], "g/m2/d", "kg/m2/s") # NPP
+    output[[6]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 9], "g/m2/d", "kg/m2/s") # Leaf Litter Flux
+    output[[7]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 11], "g/m2/d", "kg/m2/s") # Woody Litter Flux
+    output[[8]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 13], "g/m2/d", "kg/m2/s") # Root Litter Flux
 
-    ## Pools
-    output[[9]]  <- (sub.DALEC.output[, 15] * 0.001)  # Leaf Carbon, kgC/m2
-    output[[10]] <- (sub.DALEC.output[, 17] * 0.001)  # Wood Carbon, kgC/m2
-    output[[11]] <- (sub.DALEC.output[, 19] * 0.001)  # Root Carbon, kgC/m2
-    output[[12]] <- (sub.DALEC.output[, 27] * 0.001)  # Litter Carbon, kgC/m2
-    output[[13]] <- (sub.DALEC.output[, 29] * 0.001)  # Soil Carbon, kgC/m2
+    ## Pools (convert g/m2 to kg/m2)
+    output[[9]]  <- PEcAn.utils::ud_convert(sub.DALEC.output[, 15], "g/m2", "kg/m2") # Leaf Carbon
+    output[[10]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 17], "g/m2", "kg/m2") # Wood Carbon
+    output[[11]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 19], "g/m2", "kg/m2") # Root Carbon
+    output[[12]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 27], "g/m2", "kg/m2") # Litter Carbon
+    output[[13]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 29], "g/m2", "kg/m2") # Soil Carbon
 
     ## standard composites
     output[[14]] <- output[[1]] + output[[2]]  # Total Respiration

models/fates/R/write.configs.FATES.R

@@ -307,25 +307,25 @@ write.config.FATES <- function(defaults, trait.values, settings, run.id){
        # Ha activation energy for vcmax - FATES units: J/mol
        if(var == "Ha_Modified_Arrhenius_Vcmax"){
          ncdf4::ncvar_put(nc=fates.param.nc, varid='fates_vcmaxha', start = ipft, count = 1,
-                          vals=pft[v]*1000)  ## convert from kj/mol to J/mol (FATES units)
+                          vals=PEcAn.utils::ud_convert(pft[v], "kJ/mol", "J/mol"))  ## convert from kJ/mol to J/mol (FATES units)
        }
 
        # Hd deactivation energy for vcmax - FATES units: J/mol
        if(var == "Hd_Modified_Arrhenius_Vcmax"){
          ncdf4::ncvar_put(nc=fates.param.nc, varid='fates_vcmaxhd', start = ipft, count = 1,
-                          vals=pft[v]*1000)  ## convert from kj/mol to J/mol (FATES units)
+                          vals=PEcAn.utils::ud_convert(pft[v], "kJ/mol", "J/mol"))  ## convert from kJ/mol to J/mol (FATES units)
        }
 
        # Ha activation energy for Jmax - FATES units: J/mol
        if(var == "Ha_Modified_Arrhenius_Jmax"){
          ncdf4::ncvar_put(nc=fates.param.nc, varid='fates_jmaxha', start = ipft, count = 1,
-                          vals=pft[v]*1000)  ## convert from kj/mol to J/mol (FATES units)
+                          vals=PEcAn.utils::ud_convert(pft[v], "kJ/mol", "J/mol"))  ## convert from kJ/mol to J/mol (FATES units)
        }
 
        # Hd deactivation energy for Jmax - FATES units: J/mol
        if(var == "Hd_Modified_Arrhenius_Jmax"){
          ncdf4::ncvar_put(nc=fates.param.nc, varid='fates_jmaxhd', start = ipft, count = 1,
-                          vals=pft[v]*1000)  ## convert from kj/mol to J/mol (FATES units)
+                          vals=PEcAn.utils::ud_convert(pft[v], "kJ/mol", "J/mol"))  ## convert from kJ/mol to J/mol (FATES units)
        }
 
        # deltaS Vcmax - BETY units:J/mol/K;  FATES units: J/mol/K

models/gday/R/model2netcdf.GDAY.R

@@ -15,9 +15,6 @@
 model2netcdf.GDAY <- function(outdir, sitelat, sitelon, start_date, end_date) {
 
 
-  G_2_KG <- 0.001
-  TONNES_PER_HA_TO_G_M2 <- 100
-  THA_2_KG_M2 <- TONNES_PER_HA_TO_G_M2 * 0.001
 
   ### Read in model output in GDAY format
   GDAY.output <- utils::read.csv(file.path(outdir, "gday_out.csv"), header = TRUE, sep = ",", skip = 1)
@@ -44,17 +41,17 @@ model2netcdf.GDAY <- function(outdir, sitelat, sitelon, start_date, end_date) {
     output <- list()
 
     ## standard variables: C-Fluxes
-    output[[1]] <- (sub.GDAY.output[, "auto_resp"] * THA_2_KG_M2) / timestep.s
-    output[[2]] <- (sub.GDAY.output[, "hetero_resp"] * THA_2_KG_M2) / timestep.s
-    output[[3]] <- (sub.GDAY.output[, "auto_resp"] + sub.GDAY.output[, "hetero_resp"] *
-                      THA_2_KG_M2) / timestep.s
-    output[[4]] <- (sub.GDAY.output[, "gpp"] * THA_2_KG_M2) / timestep.s
-    output[[5]] <- (sub.GDAY.output[, "nep"] * -1 * THA_2_KG_M2) / timestep.s
-    output[[6]] <- (sub.GDAY.output[, "npp"] * THA_2_KG_M2) / timestep.s
+    # GDAY outputs in Mg/ha/year, convert directly to kgC/m2/s (ud_convert returns per-second)
+    output[[1]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "auto_resp"], "Mg/ha/yr", "kg/m2/s")
+    output[[2]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "hetero_resp"], "Mg/ha/yr", "kg/m2/s")
+    output[[3]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "auto_resp"] + sub.GDAY.output[, "hetero_resp"], "Mg/ha/yr", "kg/m2/s")
+    output[[4]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "gpp"], "Mg/ha/yr", "kg/m2/s")
+    output[[5]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "nep"] * -1, "Mg/ha/yr", "kg/m2/s")
+    output[[6]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "npp"], "Mg/ha/yr", "kg/m2/s")
 
-    ## standard variables: C-State
-    output[[7]] <- (sub.GDAY.output[, "stem"] + sub.GDAY.output[, "branch"] * THA_2_KG_M2) / timestep.s
-    output[[8]] <- (sub.GDAY.output[, "soilc"] * THA_2_KG_M2) / timestep.s
+    ## standard variables: C-State (pools) - divide by timestep as in original code
+    output[[7]] <- (PEcAn.utils::ud_convert(sub.GDAY.output[, "stem"], "Mg/ha", "kg/m2") + sub.GDAY.output[, "branch"] * THA_2_KG_M2) / timestep.s
+    output[[8]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "soilc"], "Mg/ha", "kg/m2") / timestep.s
     output[[9]] <- (sub.GDAY.output[, "lai"])
 
     ## standard variables: water fluxes

models/sipnet/R/model2netcdf.SIPNET.R

@@ -134,21 +134,20 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
     bounds <- round(bounds,4) 
 
     ## Setup outputs for netCDF file in appropriate units
-    output       <- list(
-      "GPP" = (sub.sipnet.output$gpp * 0.001) / timestep.s,  # GPP in kgC/m2/s
-      "NPP" = (sub.sipnet.output$gpp * 0.001) / timestep.s - ((sub.sipnet.output$rAboveground *
-                                                                 0.001) / timestep.s + (sub.sipnet.output$rRoot * 0.001) / timestep.s), # NPP in kgC/m2/s. Post SIPNET calculation
-      "TotalResp" = (sub.sipnet.output$rtot * 0.001) / timestep.s,  # Total Respiration in kgC/m2/s
-      "AutoResp" = (sub.sipnet.output$rAboveground * 0.001) / timestep.s + (sub.sipnet.output$rRoot *
-                                                                              0.001) / timestep.s,  # Autotrophic Respiration in kgC/m2/s
-      "HeteroResp" = ((sub.sipnet.output$rSoil - sub.sipnet.output$rRoot) * 0.001) / timestep.s,  # Heterotrophic Respiration in kgC/m2/s
-      "SoilResp" = (sub.sipnet.output$rSoil * 0.001) / timestep.s,  # Soil Respiration in kgC/m2/s
-      "NEE" = (sub.sipnet.output$nee * 0.001) / timestep.s,  # NEE in kgC/m2/s
-      "AbvGrndWood" = (sub.sipnet.output$plantWoodC * 0.001),  # Above ground wood kgC/m2
-      "leaf_carbon_content" = (sub.sipnet.output$plantLeafC * 0.001),  # Leaf C kgC/m2
-      "TotLivBiom" = (sub.sipnet.output$plantWoodC * 0.001) + (sub.sipnet.output$plantLeafC * 0.001) + 
-        (sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC) * 0.001, # Total living C kgC/m2
-      "TotSoilCarb" = (sub.sipnet.output$soil * 0.001) + (sub.sipnet.output$litter * 0.001)  # Total soil C kgC/m2
+    output <- list(
+      "GPP" = PEcAn.utils::ud_convert(sub.sipnet.output$gpp, "g/m2", "kg/m2") / timestep.s,
+      "NPP" = (PEcAn.utils::ud_convert(sub.sipnet.output$gpp, "g/m2", "kg/m2") -
+               PEcAn.utils::ud_convert(sub.sipnet.output$rAboveground + sub.sipnet.output$rRoot, "g/m2", "kg/m2")) / timestep.s,
+      "TotalResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rtot, "g/m2", "kg/m2") / timestep.s,
+      "AutoResp" = (PEcAn.utils::ud_convert(sub.sipnet.output$rAboveground + sub.sipnet.output$rRoot, "g/m2", "kg/m2")) / timestep.s,
+      "HeteroResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rSoil - sub.sipnet.output$rRoot, "g/m2", "kg/m2") / timestep.s,
+      "SoilResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rSoil, "g/m2", "kg/m2") / timestep.s,
+      "NEE" = PEcAn.utils::ud_convert(sub.sipnet.output$nee, "g/m2", "kg/m2") / timestep.s,
+      "AbvGrndWood" = PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC, "g/m2", "kg/m2"),
+      "leaf_carbon_content" = PEcAn.utils::ud_convert(sub.sipnet.output$plantLeafC, "g/m2", "kg/m2"),
+      "TotLivBiom" = (PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC +
+                                              sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC, "g/m2", "kg/m2")),
+      "TotSoilCarb" = PEcAn.utils::ud_convert(sub.sipnet.output$soil + sub.sipnet.output$litter, "g/m2", "kg/m2")
     )
     if (revision == "unk") {
       ## *** NOTE : npp in the sipnet output file is actually evapotranspiration, this is due to a bug in sipnet.c : ***
@@ -164,8 +163,8 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
     output[["SoilMoist"]] <- (sub.sipnet.output$soilWater * 10)  # Soil moisture kgW/m2
     output[["SoilMoistFrac"]] <- (sub.sipnet.output$soilWetnessFrac)  # Fractional soil wetness
     output[["SWE"]] <- (sub.sipnet.output$snow * 10)  # SWE
-    output[["litter_carbon_content"]] <- sub.sipnet.output$litter * 0.001  ## litter kgC/m2
-    output[["litter_mass_content_of_water"]] <- (sub.sipnet.output$litterWater * 10) # Litter water kgW/m2
+    output[["litter_carbon_content"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$litter, "g/m2", "kg/m2")
+    output[["litter_mass_content_of_water"]] <- (sub.sipnet.output$litterWater * 10)  # * 10 converts cm to mm
     #calculate LAI for standard output
     param <- utils::read.table(file.path(gsub(pattern = "/out/",
                                               replacement = "/run/", x = outdir),
@@ -174,10 +173,10 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
     leafC <- 0.48
     SLA <- 1000 / param[id, 2] #SLA, m2/kgC
     output[["LAI"]] <- output[["leaf_carbon_content"]] * SLA # LAI
-    output[["fine_root_carbon_content"]] <- sub.sipnet.output$fineRootC   * 0.001  ## fine_root_carbon_content kgC/m2
-    output[["coarse_root_carbon_content"]] <- sub.sipnet.output$coarseRootC * 0.001  ## coarse_root_carbon_content kgC/m2
-    output[["GWBI"]] <- (sub.sipnet.output$woodCreation * 0.001) / 86400 ## kgC/m2/s - this is daily in SIPNET
-    output[["AGB"]] <- (sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC) * 0.001 # Total aboveground biomass kgC/m2
+    output[["fine_root_carbon_content"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$fineRootC, "g/m2", "kg/m2")
+    output[["coarse_root_carbon_content"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$coarseRootC, "g/m2", "kg/m2")
+    output[["GWBI"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$woodCreation, "g/m2/d", "kg/m2/s")
+    output[["AGB"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC, "g/m2", "kg/m2")
     output[["time_bounds"]] <- c(rbind(bounds[,1], bounds[,2]))
 
     # ******************** Declare netCDF variables ********************#

models/sipnet/R/write_restart.SIPNET.R

@@ -47,12 +47,11 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,
 
   ## Converting to sipnet units
   prior.sla <- new.params[[which(!names(new.params) %in% c("soil", "soil_SDA", "restart"))[1]]]$SLA
-  unit.conv <- 2 * (10000 / 1) * (1 / 1000) * (3.154 * 10^7)  # kgC/m2/s -> Mg/ha/yr
 
   analysis.save <- list()
 
   if ("NPP" %in% variables) {
-    analysis.save[[length(analysis.save) + 1]] <- PEcAn.utils::ud_convert(new.state$NPP, "kg/m^2/s", "Mg/ha/yr")  #*unit.conv -> Mg/ha/yr
+    analysis.save[[length(analysis.save) + 1]] <- PEcAn.utils::ud_convert(new.state$NPP, "kg/m^2/s", "Mg/ha/yr")
     names(analysis.save[[length(analysis.save)]]) <- c("NPP")
   }
 

models/stics/R/met2model.STICS.R

@@ -140,7 +140,7 @@ met2model.STICS <- function(in.path, in.prefix, outfolder, start_date, end_date,
       # column 10: rainfall (mm.j-1)
       Rain  <- ncdf4::ncvar_get(nc, "precipitation_flux") # kg m-2 s-1
       Rain  <- Rain[ydays %in% simdays]
-      raini <- tapply(Rain * 86400, ind, mean, na.rm = TRUE) 
+      raini <- tapply(PEcAn.utils::ud_convert(Rain, "kg/m2/s", "kg/m2/d"), ind, mean, na.rm = TRUE) # 1 kg/m2 = 1 mm water
       weather_df[ ,10] <- round(raini, digits = 2) # precipitation (mm d-1)	
 
       # column 11: wind (m.s-1) 

modules/data.land/R/IC_BADM_Utilities.R

@@ -226,62 +226,45 @@ netcdf.writer.BADM <- function(lat, long, siteid, outdir, ens){
 
 #' BADM_IC_process
 #'
-#' @param settings pecan xml settings
-#' @param dir output dir which you want to store the IC netcdf file
+#' @param settings Single PEcAn Settings object (list with 'run' element)
+#' @param dir Output directory for storing IC netcdf file
 #' @param overwrite Flag for overwriting the IC file.
 #'
+#' @description Process initial conditions for a single site using BADM data.
+#' This function accepts only single Settings objects.
+#' For processing multiple sites, use a higher-level function.
+#'
 #' @return a list of paths to generated and stored IC files.
 #' @export
 #'
-BADM_IC_process <- function(settings, dir, overwrite=TRUE){
-
-  # check if this is a single-site or multi-site configuration
-  if ("run" %in% names(settings)) {
-    settings <- list(settings)
+BADM_IC_process <- function(settings, dir, overwrite = TRUE) {
+  # This function now only accepts single Settings objects
+  if (!("run" %in% names(settings))) {
+    stop("settings must be a single Settings object with 'run' element")
   }
 
-  # create site info.
-  new.site <- 
-    settings %>% 
-    purrr::map(~.x[['run']] ) %>% 
-    purrr::map('site')%>% 
-    purrr::map(function(site.list){
-      #conversion from string to number
-      site.list$lat <- as.numeric(site.list$lat)
-      site.list$lon <- as.numeric(site.list$lon)
-      list(id=site.list$id, lat=site.list$lat, lon=site.list$lon)
-    })%>% 
-    dplyr::bind_rows() %>% 
-    as.list()
-
-  # process each site configuration
-  out.ense <- list()
-
-  for (i in seq_along(settings)) {
-    site.settings <- settings[[i]]
-    ens.size <- ens.size <- max(1, site.settings$ensemble$size %||% 1)
-
-    # get site info for this specific site
-    site.info <- list(
-      id = new.site$id[i],
-      lat = new.site$lat[i], 
-      lon = new.site$lon[i]
-    )
-
-    site.outputs <- seq_len(ens.size) %>%
-      purrr::map(~ netcdf.writer.BADM(site.info$lat,
-                                      site.info$lon,
-                                      site.info$id,
-                                      outdir=dir,
-                                      ens=.x))
-
-    site.outputs <- site.outputs %>%
-      stats::setNames(rep("path", length(site.outputs)))
-
-    out.ense <- c(out.ense, site.outputs)
+  # extract and normalise site info
+  site.info <- settings[['run']][['site']]
+  site.id <- site.info$id
+  site.lat <- suppressWarnings(as.numeric(site.info$lat))[1]
+  site.lon <- suppressWarnings(as.numeric(site.info$lon))[1]
+
+  if (is.na(site.lat) || is.na(site.lon)) {
+    stop("Invalid site coordinates")
   }
-
-  return(out.ense)
+
+  ens.size <- max(1, (settings$ensemble$size %||% 1))
+
+  # generate ensemble IC files for this site
+  site.outputs <- seq_len(ens.size) %>%
+    purrr::map(~ netcdf.writer.BADM(site.lat,
+                                    site.lon,
+                                    site.id,
+                                    outdir = dir,
+                                    ens = .x))
+
+  # name each entry "path" to match previous behaviour
+  stats::setNames(site.outputs, rep("path", length(site.outputs)))
 }
 
 #' EPA_ecoregion_finder

modules/data.land/tests/testthat/test-IC_BADM_Utilities.R

@@ -43,51 +43,6 @@ test_that("netcdf.writer.BADM creates output directory if missing", {
   expect_true(dir.exists(tmp_outdir))
 })
 
-test_that("BADM_IC_process generates correct number of ensemble files for single-site", {
-  settings <- list(
-    run = list(site = list(id = "US-Ha1", lat = 42.5378, lon = -72.1715)),
-    ensemble = list(size = 3)
-  )
-  out_files <- BADM_IC_process(settings, dir = tempdir(), overwrite = TRUE)
-  expect_length(out_files, 3)
-  expect_true(all(file.exists(unlist(out_files))))
-})
-
-test_that("BADM_IC_process generates correct number of ensemble files for multi-site", {
-  settings <- list(
-    list(
-      run = list(site = list(id = "US-Ha1", lat = 42.5378, lon = -72.1715)),
-      ensemble = list(size = 2)
-    ),
-    list(
-      run = list(site = list(id = "US-WCr", lat = 45.805925, lon = -90.07961)),
-      ensemble = list(size = 3)
-    )
-  )
-  out_files <- BADM_IC_process(settings, dir = tempdir(), overwrite = TRUE)
-  expect_length(out_files, 5) 
-  expect_true(all(file.exists(unlist(out_files))))
-})
-
-test_that("BADM_IC_process handles missing or malformed settings gracefully", {
-  settings <- list(
-    list(
-      run = list(site = list(id = "US-Ha1", lat = NA, lon = -72.1715)),
-      ensemble = list(size = 1)
-    )
-  )
-  expect_error(BADM_IC_process(settings, dir = tempdir(), overwrite = TRUE))
-})
-
-test_that("BADM_IC_process handles missing ensemble size with fallback", {
-  settings <- list(
-    run = list(site = list(id = "US-Ha1", lat = 42.5378, lon = -72.1715)),
-    ensemble = list(size = 0)  
-  )
-  out_files <- BADM_IC_process(settings, dir = tempdir(), overwrite = TRUE)
-  expect_length(out_files, 1)  
-})
-
 test_that("Read.IC.info.BADM returns empty dataframe for invalid coordinates", {
   expect_error(Read.IC.info.BADM(999, 999))
 })