Skip to content

Refactor SIPNET to separate version logic from availability checks #3862

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Thomas-Sedhom wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Refactor SIPNET to separate version logic from availability checks #3862

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 46 additions & 11 deletions models/sipnet/R/model2netcdf.SIPNET.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,29 @@ mergeNC <- function(
##' @author Shawn Serbin, Michael Dietze
model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, delete.raw = FALSE, revision, prefix = "sipnet.out",
overwrite = FALSE, conflict = FALSE) {
# Version-specific capabilities (separate from output availability)
rev_raw <- revision
legacy_v1 <- c("102319", "136", "r136", "ssr", "git")
rev_known <- FALSE
sipnet_version <- NA
if (!is.null(rev_raw) && nzchar(rev_raw)) {
if (rev_raw %in% legacy_v1) {
sipnet_version <- numeric_version("1.0")
rev_known <- TRUE
} else {
rev_clean <- sub("^v", "", rev_raw, ignore.case = TRUE)
sipnet_version <- numeric_version(rev_clean, strict = FALSE)
if (!is.na(sipnet_version)) {
rev_known <- TRUE
}
}
}
rev_str <- if (rev_known && sipnet_version >= "2.0") "v2" else if (rev_known) "v1" else "unknown"
caps <- list(
has_litter_water = if (rev_known) rev_str == "v1" else TRUE,
has_n_cycle = if (rev_known) rev_str == "v2" else TRUE
)

### Read in model output in SIPNET format
sipnet_out_file <- file.path(outdir, prefix)
# SIPNET v1 had a "Notes" comment line before the header; v2 removed it.
Expand Down Expand Up @@ -106,9 +129,21 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
file.rename(file.path(outdir, paste(y, "nc", sep = ".")), file.path(outdir, "previous.nc"))
}
print(paste("---- Processing year: ", y)) # turn on for debugging

## Subset data for processing
sub.sipnet.output <- subset(sipnet_output, sipnet_output$year == y)
avail <- list(
litterWater = "litterWater" %in% names(sub.sipnet.output),
woodCreation = "woodCreation" %in% names(sub.sipnet.output),
minN = "minN" %in% names(sub.sipnet.output),
soilOrgN = "soilOrgN" %in% names(sub.sipnet.output),
litterN = "litterN" %in% names(sub.sipnet.output),
n2o = "n2o" %in% names(sub.sipnet.output),
nLeaching = "nLeaching" %in% names(sub.sipnet.output),
nFixation = "nFixation" %in% names(sub.sipnet.output),
nUptake = "nUptake" %in% names(sub.sipnet.output),
ch4 = "ch4" %in% names(sub.sipnet.output)
)

raw_time <- sub.sipnet.output[["time"]] # decimal hours (eg 13.75 = 1:45 PM)
doy <- sub.sipnet.output[["day"]] # day of year, not of month
Expand Down Expand Up @@ -167,7 +202,7 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
output[["SWE"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$snow, "cm", "mm") # Snow Water Equivalent
output[["litter_carbon_content"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$litter, "g/m2", "kg/m2")
# litterWater was removed in SIPNET v2; only extract if present
if ("litterWater" %in% names(sub.sipnet.output)) {
if (caps$has_litter_water && avail$litterWater) {
output[["litter_mass_content_of_water"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$litterWater, "cm", "mm")
}

Expand All @@ -183,34 +218,34 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
output[["LAI"]] <- output[["leaf_carbon_content"]] * SLA
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")
if ("woodCreation" %in% names(sub.sipnet.output)) {
if (avail$woodCreation) {
output[["GWBI"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$woodCreation, "g/m2/day", "kg/m2/s")
}
output[["AGB"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC, "g/m2", "kg/m2")
# columns only present in sipnet >= v2 with N and methane turned on
if ("minN" %in% names(sub.sipnet.output)) {
if (caps$has_n_cycle && avail$minN) {
output[["mineral_N"]] <- sub.sipnet.output$minN * 0.001 # gN/m2 -> kgN/m2
}
if ("soilOrgN" %in% names(sub.sipnet.output)) {
if (caps$has_n_cycle && avail$soilOrgN) {
output[["soil_organic_N"]] <- sub.sipnet.output$soilOrgN * 0.001 # kgN/m2
}
if ("litterN" %in% names(sub.sipnet.output)) {
if (caps$has_n_cycle && avail$litterN) {
output[["litter_N"]] <- sub.sipnet.output$litterN * 0.001 # kgN/m2
}
if ("n2o" %in% names(sub.sipnet.output)) {
if (caps$has_n_cycle && avail$n2o) {
output[["N2O_flux"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$n2o, "g/m2", "kg/m2") / timestep.s
# convert g N m-2 per timestep -> kg N m-2 s-1
}
if ("nLeaching" %in% names(sub.sipnet.output)) {
if (caps$has_n_cycle && avail$nLeaching) {
output[["N_leaching"]] <- (sub.sipnet.output$nLeaching * 0.001) / timestep.s # kgN/m2/s
}
if ("nFixation" %in% names(sub.sipnet.output)) {
if (caps$has_n_cycle && avail$nFixation) {
output[["N_fixation"]] <- (sub.sipnet.output$nFixation * 0.001) / timestep.s # kgN/m2/s
}
if ("nUptake" %in% names(sub.sipnet.output)) {
if (caps$has_n_cycle && avail$nUptake) {
output[["N_uptake"]] <- (sub.sipnet.output$nUptake * 0.001) / timestep.s # kgN/m2/s
}
if ("ch4" %in% names(sub.sipnet.output)) {
if (caps$has_n_cycle && avail$ch4) {
output[["CH4_flux"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$ch4, "g/m2", "kg/m2") / timestep.s
# convert g C m-2 per timestep -> kg C m-2 s-1
}
Expand Down
33 changes: 24 additions & 9 deletions models/sipnet/R/write.configs.SIPNET.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,16 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
}
}
rev_str <- if (sipnet_version >= "2.0") "v2" else "v1"

# Version-specific capabilities (kept separate from input availability)
caps <- list(
has_runtime_flags = rev_str == "v2",
has_param_spatial = rev_str == "v1",
has_m_ballBerry = rev_str == "v1",
has_cold_soil_resp = rev_str == "v1",
has_litWaterDrainRate = rev_str == "v1",
has_litterWFracInit = rev_str == "v1",
has_microbeInit = rev_str == "v1"
)

### WRITE sipnet.in
template.in <- system.file(
Expand All @@ -40,7 +49,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
writeLines(config.text, con = file.path(settings$rundir, run.id, "sipnet.in"))

# for v2, allow settings$model to override runtime flags in sipnet.in
if (rev_str == "v2") {
if (caps$has_runtime_flags) {
sipnet_in_path <- file.path(settings$rundir, run.id, "sipnet.in")
nc <- as.integer(settings$model$nitrogen_cycle %||% 1)
lp <- as.integer(settings$model$litter_pool %||% 1)
Expand Down Expand Up @@ -190,7 +199,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs


### WRITE *.param-spatial
if (rev_str == "v1") {
if (caps$has_param_spatial) {
template.paramSpatial <- system.file("template.param-spatial", package = "PEcAn.SIPNET")
file.copy(template.paramSpatial, file.path(settings$rundir, run.id, "sipnet.param-spatial"))
}
Expand Down Expand Up @@ -346,7 +355,9 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
}

# Ball-berry stomatal slope parameter m (v1 only; m_ballBerry removed in v2)
if ("stomatal_slope.BB" %in% pft.trait.names && "m_ballBerry" %in% param[, 1]) {
if (caps$has_m_ballBerry &&
"stomatal_slope.BB" %in% pft.trait.names &&
"m_ballBerry" %in% param[, 1]) {
id <- which(param[, 1] == "m_ballBerry")
param[id, 2] <- pft.traits[which(pft.trait.names == "stomatal_slope.BB")]
}
Expand Down Expand Up @@ -485,7 +496,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
# These results in improved winter soil respiration values
# they don't affect anything when the seasonal soil respiration functionality in SIPNET is turned-off
# 2025-07-22 CKB: soilRespQ10Cold and baseSoilRespCold were removed from Sipnet V2.0
if (rev_str == "v1") {
if (caps$has_cold_soil_resp) {
# assume soil resp Q10 cold == soil resp Q10
param[which(param[, 1] == "soilRespQ10Cold"), 2] <- param[which(param[, 1] == "soilRespQ10"), 2]
# default SIPNET prior of baseSoilRespCold was 1/4th of baseSoilResp
Expand Down Expand Up @@ -641,7 +652,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
}
if ("soil_hydraulic_conductivity_at_saturation" %in% names(soil_IC_list$vals)) {
#litwaterDrainrate in cm/day (v1 only; litWaterDrainRate removed in v2)
if ("litWaterDrainRate" %in% param[, 1]) {
if (caps$has_litWaterDrainRate && "litWaterDrainRate" %in% param[, 1]) {
param[which(param[, 1] == "litWaterDrainRate"), 2] <- PEcAn.utils::ud_convert(unlist(soil_IC_list$vals["soil_hydraulic_conductivity_at_saturation"])[1], "m s-1", "cm day-1")
}
}
Expand Down Expand Up @@ -692,7 +703,9 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
param[which(param[, 1] == "soilInit"), 2] <- IC$soil
}
## litterWFracInit fraction (v1 only; removed in v2)
if ("litter_mass_content_of_water" %in% ic.names && "litterWFracInit" %in% param[, 1]) {
if (caps$has_litterWFracInit &&
"litter_mass_content_of_water" %in% ic.names &&
"litterWFracInit" %in% param[, 1]) {
#here we use litterWaterContent/litterWHC to calculate the litterWFracInit
param[which(param[, 1] == "litterWFracInit"), 2] <- IC$litter_mass_content_of_water/(param[which(param[, 1] == "litterWHC"), 2]*10)
}
Expand All @@ -709,7 +722,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
param[which(param[, 1] == "snowInit"), 2] <- IC$SWE
}
## microbeInit mgC/g soil (v1 only; removed in v2)
if ("microbe" %in% ic.names && "microbeInit" %in% param[, 1]) {
if (caps$has_microbeInit && "microbe" %in% ic.names && "microbeInit" %in% param[, 1]) {
param[which(param[, 1] == "microbeInit"), 2] <- IC$microbe
}

Expand Down Expand Up @@ -826,7 +839,9 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
param[param[, 1] == "leafOffDay", 2] <- leafOffDay
}
}
if (ic_has_ncvars[["Microbial Biomass C"]] && "microbeInit" %in% param[, 1]) {
if (caps$has_microbeInit &&
ic_has_ncvars[["Microbial Biomass C"]] &&
"microbeInit" %in% param[, 1]) {
microbe <- ncdf4::ncvar_get(IC.nc, "Microbial Biomass C")
if (!is.na(microbe) && is.numeric(microbe)) {
param[param[, 1] == "microbeInit", 2] <- PEcAn.utils::ud_convert(microbe, "mg kg-1", "mg g-1") #BETY: mg microbial C kg-1 soil
Expand Down
Loading