My personal function works independently but not when I call it in my overall code

[matthewlombardo02-dr]

Hello everyone,

I am using R to create my calibration code for my research project using the model called SWAT+. My code takes samples 27 different parameters and then replaces the old values for those parameters with the new ones and then runs SWAT+, and does this 1,000 times. SWAT+ holds it's data in different text files, and for it to run correctly the code must replace the original values with new values without changing the formatting of the text files. If the formatting is off by even a little, SWAT won't run.
My function of interest, "SWAT_SOILSK_Mod", successfully replaced the soil_k value in the text file called soils.sol when I test it independently. Soils.sol contains all the information for my different soil types. I've attached an example of it being successful. However, when I call on the function in my overall code, called "SWAT_Calibration" it doesn't work and messes up the file. I've also attached an example of this. In the overall code there is another function that edits the same text file called "SWAT_SOILSAWC_Mod" and when I ran the functions myself sequentially there was no issue. There are no other functions that edit this text file. My SOILSAWC_MOD replaces the the soils.sol with a template that is the exact same as soils.sol, but with different awc values, so the formatting is exactly the same. Does anyone know why my function does not work when called on by the overall code? There are no mixups in variable arguments. Since I'm a new user and I can only add one piece of media, I've pasted my functions below and attached a picture of the correct and broken soils.sol files. Since my overall code is so large, I only included the section that calls on the function of interest and a little extra. I've tried to use ChatGPT to help solve my problem, but I've been at it for a week and made no progress.

Thank you for you attention and help.
Please let me know if more information is needed.

SWAT_SOILSK_MOD:

SWAT_SOILSK_Mod <- function(path,template,k_mod_value) {

  #path = paste0("D:/SWAT+/Calibration/SWAT_Watershed_Tool_NoLake/", watershed, "/TxtInOut")
  soils_file <- file.path(path, "soils.sol")
 
  
  # 1. Read soils.sol
  soils <- readLines(soils_file, warn = FALSE)
  
  # Layer lines = indented numeric rows that start with dp (e.g. "    25.40000")
  # Soil header lines start with "33", "15", etc. (no decimal) and are skipped.
  layer_idx <- grep("^\\s*[0-9]+\\.[0-9]{5}", soils)
  
  # 2. Read original soil_k table and keep only layer rows
  #    (rows with numeric dp; header rows have NA in dp, isolates where the soil_k will be inserted)
  soil_tbl  = read.csv(paste(template,"SWAT_SOILSK_Table.csv",sep=""))
  
  # Trim column names so they become "dp", "bd", "awc", "soil_k", etc.
  names(soil_tbl) <- trimws(names(soil_tbl))
  
  # Keep only layer rows (dp not NA)
  soil_layers <- subset(soil_tbl, !is.na(dp))
  
  # Sanity check: number of CSV layer rows vs number of layer lines in soils.sol
  if (nrow(soil_layers) != length(layer_idx)) {
    stop(sprintf(
      "Mismatch: %d layer lines in soils.sol but %d rows in SWAT_SoilSK_Table.csv",
      length(layer_idx), nrow(soil_layers)
    ))
  }
  
  # 3. Build new soil_k values
  k_new <- soil_layers$soil_k * k_mod_value
  
  # Enforce 12–1607 max and min bounds
  k_new <- pmax(12, pmin(1607, k_new))
  
  # Format exactly 10 characters wide with 5 decimals
  k_new_fmt <- sprintf("%10.5f", k_new)
  

  # 4. Replace ONLY the soil_k field in-place on each layer line
  #   soil_k = cols 169–178 (10 characters wide)
  soil_start <- 169L            # first character of soil_k
  soil_width <- 10L
  soil_end   <- soil_start + soil_width - 1L
  
  for (i in seq_along(layer_idx)) {
    line <- soils[layer_idx[i]]
    
    # Replace only the soil_k field; everything else stays identical
    substr(line, soil_start, soil_end) <- k_new_fmt[i]
    
    soils[layer_idx[i]] <- line
  }
  
  # 5. Backup and write modified soils.sol-
  backup_file <- file.path(path, "soils_before_SOILK_MOD.sol")
  file.copy(soils_file, backup_file, overwrite = TRUE)
  
  writeLines(soils, soils_file, useBytes = TRUE)
  
  invisible(NULL)
}

SWAT_SOILSAWC_MOD:

SWAT_SOILSAWC_Mod <- function(path,template,AWC_Value) {
  
  awc <- sprintf("%07.5f", AWC_Value)
  
  #Review Note: put paths in dynamically. Will avoid simple mistakes
  awctable <- readLines(paste(template,"soils-awc_template.txt",sep="")) |>
    stringr::str_replace_all(c("changme" = awc))
  writeLines(awctable, con = paste(path,"/soils.sol",sep=""))
}

Part of SWAT_Calibration:

#Run for loop to run calibration
for(j in restart:iterations){
  
  #1a. Resample AWC
  Output$awc[j] <- runif(1, min = 0, max = 0.7)
  SWAT_SOILSAWC_Mod(txtinoutpath,template,Output$awc[j])

  #1b. Resample Soil K
  Output$k_mult[j] <- runif(1, min = 0.5, max = 1.5)
  SWAT_SOILSK_Mod(txtinoutpath, template,Output$soilk_mult[j])
  
  
  #2. Resample bf_max, flo_min, NO3_N, and HL_NO3N
  Output$no3_n[j] <- runif(1, min = 0, max = 1000)
  Output$bf_max[j] <- runif(1, min = 0.5, max = 2)
  Output$hl_no3n[j] <- runif(1, min = 0, max = 200)
  Output$flo_min[j]  <- runif(1, min = 0.1, max = 10)
  SWAT_AQUIFER_Mod(txtinoutpath,template,Output$no3_n[j],Output$bf_max[j],Output$hl_no3n[j],Output$flo_min[j])
  
  #3a. Resample bed_K
  Output$bed_K[j] <- runif(1, min = 0, max = 500)
  SWAT_CHANHYD_Mod(txtinoutpath,template,Output$bed_K[j])
  
  #3b. Resample channel nutrient parameters
  Output$plt_n[j] <- runif(1, min = 0, max = 100)
  Output$ben_nh3n[j] <- runif(1, min = 0, max = 1)
  Output$ptln_stl[j] <- runif(1, min = 0.001, max = 0.1)
  Output$nh3n_no2n[j]  <- runif(1, min = 0.1, max = 1)
  Output$no2n_no3n[j] <- runif(1, min = 0.2, max = 2)
  Output$ptln_nh3n[j] <- runif(1, min = 0.2, max = 0.4)
  SWAT_CHANNUTR_Mod(txtinoutpath,template,Output$plt_n[j],Output$ben_nh3n[j],Output$ptln_stl[j],Output$nh3n_no2n[j],Output$no2n_no3n[j],Output$ptln_nh3n[j])
  
  #4. Resample SWAT Curve Number
  Output$cn2_mult[j] <- runif(1, min = 0.5, max = 1.5)
  SWAT_CN_Mod(txtinoutpath,template,Output$cn2_mult[j])
  
  #5. Resample SWAT Hydrology
  Output$lat_time[j] <- runif(1,1,180)
  Output$esco[j] <- runif(1, min = 0.01, max = 1)
  Output$epco[j] <- runif(1, min = 0.01, max = 1)
  Output$orgn_enrich[j] <- runif(1, min = 0, max = 1)
  SWAT_HYDROLOGY_Mod(txtinoutpath,template,Output$lat_time[j],Output$esco[j],Output$epco[j],Output$orgn_enrich[j])
  
  #6. Resample snowmelt parameters
  Output$fall_tmp[j] <- runif(1, min = -5, max = 5)
  Output$melt_tmp[j] <- runif(1, min = -5, max = 5)
  Output$melt_max[j] <- runif(1, min = 0, max = 10)
  Output$melt_min[j] <- runif(1, min = 0, max = Output$melt_max[j])
  Output$melt_lag[j] <- runif(1, min = 0.01, max = 1)
  SWAT_SNOW_Mod(txtinoutpath,template,Output$fall_tmp[j],Output$melt_tmp[j],Output$melt_max[j],Output$melt_min[j],Output$melt_lag[j])
  
  #7. Resample surface lag and nutrient parameters
  Output$surlag[j] <- runif(1, min = 0.05, max = 24)
  Output$orgn_min[j] <- runif(1, min = 0.00001, max = 0.06)
  Output$n_uptake[j] <- runif(1, min = 0, max = 100)
  Output$n_perc[j] <- runif(1, min = 0.01, max = 1)
  Output$rsd_decomp[j] <- runif(1, min = 0.02, max = 0.1)
  Output$denit_exp[j] <- runif(1, min = 0, max = 3)
  Output$denit_frac[j] <- runif(1, min = 0, max = 1)
  SWAT_BASIN_Mod(txtinoutpath,template,Output$surlag[j],Output$orgn_min[j],Output$n_uptake[j],Output$n_perc[j],
                 Output$rsd_decomp[j],Output$denit_exp[j],Output$denit_frac[j])

  #7. Run SWAT Model
  system(paste(txtinoutpath,"/swatplus-61.0.2-ifx-win_amd64-Rel.exe",sep=""), intern = TRUE)
  Sys.sleep(0.5)