fix(SIPNET): apply plantStorageNInit default conditionally (#4114)

models/sipnet/NEWS.md

@@ -20,15 +20,8 @@
     trait mappings alongside the existing `leafOnReallocFrac`.
     `plantStorageNInit` is deliberately not trait mapped; it is an initial
     condition, set from the IC layer like the other Init parameters.
-* For SIPNET v2, `write.config.SIPNET` now sets `plantStorageNInit` to the N
-    needed for one carbon-limited leaf flush, computed from the initialized
-    wood pool, instead of leaving the template value. With the nitrogen cycle
-    on, a zero initial storage N suppressed the first leaf-on event entirely
-    for vegetation initialized with standing wood (e.g. established orchards).
-    Sipnet takes the smaller of its carbon and nitrogen limiters, so sizing N
-    to the carbon-capped transfer removes the suppression without adding
-    surplus N. A `plantStorageNInit` supplied through the `IC` argument still
-    takes precedence, and zero-wood (annual) initializations still get 0.
+* For SIPNET v2, `write.config.SIPNET` now sets `plantStorageNInit` from initial conditions
+    (`IC` argument or poolinitcond file) if provided (#4114, @coder-Yash886).
 * More flexible handling of Sipnet "restart" files (checkpoints containing full
   model state at the end of the run):
     - New setting `settings$model$copy.restart` is a logical runtime flag

models/sipnet/R/write.configs.SIPNET.R

@@ -62,10 +62,14 @@
 #' defines an option specified in settings$model$options, its value will be
 #' updated in place; options not already in the file will be added to the bottom.
 #'
-#' For SIPNET v2, plantStorageNInit is set to the additional N needed for one
-#' carbon-limited leaf flush, computed from the initialized wood pool, so the
-#' first leaf-on event is not nitrogen suppressed. A plantStorageNInit supplied
-#' through the IC argument takes precedence over the computed value.
+#' # Parameter precedence
+#'
+#' Parameters are initialized from template parameter files (e.g. `inst/template.param_v2`
+#' for SIPNET v2). Values in
+#' the template or `defaults$constants` are used unless overridden. PFT trait values
+#' in `trait.values` override constants and template defaults. Initial condition
+#' values passed via `IC` or `poolinitcond` take highest precedence, including
+#' explicit zero values.
 #'
 #' @param defaults nested list of named constant parameter values. The
 #' structure is `list(list(constants = list(trait1 = <value>, trait2 = <value>, ...)))`.
@@ -118,6 +122,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   # Version-specific capabilities (kept separate from input availability)
   caps <- list(
     has_runtime_flags = rev_str == "v2",
+    has_n_cycle = rev_str == "v2",
     has_param_spatial = rev_str == "v1",
     has_m_ballBerry = rev_str == "v1",
     has_cold_soil_resp = rev_str == "v1",
@@ -894,6 +899,10 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     if (caps$has_microbeInit && "microbe" %in% ic.names && "microbeInit" %in% param[, 1]) {
       param[which(param[, 1] == "microbeInit"), 2] <- IC$microbe
     }
+    ## plantStorageNInit gN/m2 (v2 only)
+    if (caps$has_n_cycle && "plantStorageNInit" %in% ic.names) {
+      param[which(param[, 1] == "plantStorageNInit"), 2] <- IC$plantStorageNInit
+    }
 
   } else if (length(settings$run$inputs$poolinitcond$path) > 0) {
     IC.path <- settings$run$inputs$poolinitcond$path
@@ -920,7 +929,8 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
         "SWE",
         "date_of_budburst",
         "date_of_senescence",
-        "Microbial Biomass C"
+        "Microbial Biomass C",
+        "plantStorageNInit"
       )
       ic_has_ncvars <- ic_ncvars_to_try %in% names(IC.nc$var)
       names(ic_has_ncvars) <- ic_ncvars_to_try
@@ -1020,6 +1030,12 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
           param[param[, 1] == "microbeInit", 2] <- PEcAn.utils::ud_convert(microbe, "mg kg-1", "mg g-1") #BETY: mg microbial C kg-1 soil
         }
       }
+      if (caps$has_n_cycle && ic_has_ncvars[["plantStorageNInit"]]) {
+        storageN <- ncdf4::ncvar_get(IC.nc, "plantStorageNInit")
+        if (!is.na(storageN) && is.numeric(storageN)) {
+          param[param[, 1] == "plantStorageNInit", 2] <- storageN
+        }
+      }
 
       ncdf4::nc_close(IC.nc)
     } else {
@@ -1028,32 +1044,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   } else {
     #some stuff about IC file that we can give in lieu of actual ICs
   }
-
-  if ("plantStorageNInit" %in% param[, 1]) {
-    # If not supplied, set initial storage N to support leaf flush (g N/m2).
-    # Sipnet takes the smaller of its C and N limiters, so N sized to the
-    # carbon-limited flush removes N suppression without adding surplus N.
-    # check carbon available for initial leaf on
-    available_carbon <- param[param[, 1] == "plantWoodInit", 2] *
-      (1 - param[param[, 1] == "fineRootFrac", 2]) *
-      param[param[, 1] == "leafOnReallocFrac", 2]
 
-    # actual amount of carbon used for leaf on
-    leaf_carbon <- min(
-      param[param[, 1] == "leafGrowth", 2],
-      available_carbon
-    )
-
-    # N required to support leaf growth
-    storage_n <- leaf_carbon * max(
-      0,
-      1 / param[param[, 1] == "leafCN", 2] -
-        1 / param[param[, 1] == "woodCN", 2]
-    )
-
-    param[param[, 1] == "plantStorageNInit", 2] <-
-      IC[["plantStorageNInit"]] %||% storage_n
-  }
 
   if (!is.null(settings$run$inputs$soilmoisture)) {
     #read soil moisture netcdf file, grab closet date to start_date, set equal to soilWFrac

models/sipnet/inst/template.param_v2

@@ -4,7 +4,7 @@ litterInit		280
 soilInit		10000
 soilWFracInit	0.5
 snowInit		1.0
-plantStorageNInit 0.0
+plantStorageNInit 5.0
 fineRootFrac	0.200000
 coarseRootFrac	0.200000
 aMax			112.0

models/sipnet/tests/testthat/test-write.config.SIPNET.R

@@ -68,6 +68,66 @@ test_that("write.config.SIPNET", {
 
 })
 
+test_that("write.config.SIPNET plantStorageNInit precedence in v2", {
+  pth <- withr::local_tempdir()
+  dir.create(file.path(pth, "run", "run_v2_default"), recursive = TRUE)
+
+  s_v2 <- PEcAn.settings::as.Settings(
+    list(
+      outdir = file.path(pth, "out"),
+      rundir = file.path(pth, "run"),
+      pfts = list(pft1 = list()),
+      model = list(binary = "", revision = "2.0"),
+      run = list(
+        site = list(name = "site1", lat = 40, lon = -88),
+        inputs = list(met = list(path = "")),
+        start.date = "2025-01-01",
+        end.date = "2025-01-02"
+      ),
+      host = list(
+        name = "",
+        outdir = file.path(pth, "out"),
+        rundir = file.path(pth, "run")
+      )
+    )
+  )
+
+  # 1. Default template value (5.0) is preserved without events file or special computation
+  write.config.SIPNET(
+    defaults = list(pft1 = list(constants = list())),
+    trait.values = list(pft1 = list()),
+    settings = s_v2,
+    run.id = "run_v2_default"
+  )
+  param_default <- readLines(file.path(pth, "run", "run_v2_default", "sipnet.param"))
+  expect_match(param_default, "plantStorageNInit 5", fixed = TRUE, all = FALSE)
+
+  # 2. Overriding via IC (including an intentional zero) is respected
+  dir.create(file.path(pth, "run", "run_v2_ic_zero"), recursive = TRUE)
+  write.config.SIPNET(
+    defaults = list(pft1 = list(constants = list())),
+    trait.values = list(pft1 = list()),
+    settings = s_v2,
+    run.id = "run_v2_ic_zero",
+    IC = list(plantStorageNInit = 0)
+  )
+  param_ic_zero <- readLines(file.path(pth, "run", "run_v2_ic_zero", "sipnet.param"))
+  expect_match(param_ic_zero, "plantStorageNInit 0", fixed = TRUE, all = FALSE)
+
+  # 3. Explicit IC value non-zero
+  dir.create(file.path(pth, "run", "run_v2_ic_val"), recursive = TRUE)
+  write.config.SIPNET(
+    defaults = list(pft1 = list(constants = list())),
+    trait.values = list(pft1 = list()),
+    settings = s_v2,
+    run.id = "run_v2_ic_val",
+    IC = list(plantStorageNInit = 12.5)
+  )
+  param_ic_val <- readLines(file.path(pth, "run", "run_v2_ic_val", "sipnet.param"))
+  expect_match(param_ic_val, "plantStorageNInit 12.5", fixed = TRUE, all = FALSE)
+})
+
+
 
 test_that("update_flag_lines", {
   txt <- c("!comment", "NITROG = 0", "GDD = 1")
@@ -160,84 +220,3 @@ test_that("stem and root respiration are written as per year rates", {
   expect_equal(param_value("baseVegResp"), 1 * g_per_year * 2.0^(-2.5))
   expect_equal(param_value("baseFineRootResp"), 10 * g_per_year * 2.6^(-2.5))
 })
-
-test_that("plantStorageNInit is sized to the carbon-limited leaf flush", {
-  # scaffold for a v2 run; the irrig-only event file keeps internal phenology
-  make_settings <- function(pth) {
-    event_src_path <- file.path(pth, "events-a.in")
-    dir.create(file.path(pth, "run", "run1"), recursive = TRUE)
-    writeLines("2025 1 irrig 0 1", con = event_src_path)
-    PEcAn.settings::as.Settings(
-      list(
-        outdir = file.path(pth, "out"),
-        rundir = file.path(pth, "run"),
-        pfts = list(pft1 = list()),
-        model = list(binary = "", revision = "v2.2.0"),
-        run = list(
-          site = list(name = "site1", lat = 40, lon = -88),
-          inputs = list(
-            met = list(path = ""),
-            events = list(path = event_src_path)
-          ),
-          start.date = "2025-01-01",
-          end.date = "2025-01-02"
-        ),
-        host = list(
-          name = "",
-          outdir = file.path(pth, "out"),
-          rundir = file.path(pth, "run")
-        )
-      )
-    )
-  }
-  read_param <- function(pth) {
-    utils::read.table(file.path(pth, "run", "run1", "sipnet.param"),
-                      stringsAsFactors = FALSE)
-  }
-  pval <- function(p, name) p[p[[1]] == name, 2]
-
-  # expected value from the same template the config writer reads
-  tmpl <- utils::read.table(
-    system.file("template.param_v2", package = "PEcAn.SIPNET"),
-    stringsAsFactors = FALSE
-  )
-  tval <- function(name) tmpl[tmpl[[1]] == name, 2]
-  expected_n <- min(tval("leafGrowth"),
-                    tval("plantWoodInit") * (1 - tval("fineRootFrac")) *
-                      tval("leafOnReallocFrac")) *
-    max(0, 1 / tval("leafCN") - 1 / tval("woodCN"))
-
-  # computed default from the template wood pool
-  pth <- withr::local_tempdir()
-  write.config.SIPNET(
-    defaults = list(pft1 = list(constants = list(SLA = 2.0))),
-    trait.values = list(pft1 = list()),
-    settings = make_settings(pth), run.id = "run1"
-  )
-  expect_equal(pval(read_param(pth), "plantStorageNInit"), expected_n)
-  expect_gt(expected_n, 0)
-
-  # an IC supplied value takes precedence over the computed default
-  pth <- withr::local_tempdir()
-  write.config.SIPNET(
-    defaults = list(pft1 = list(constants = list(SLA = 2.0))),
-    trait.values = list(pft1 = list()),
-    settings = make_settings(pth), run.id = "run1",
-    IC = data.frame(plantStorageNInit = 3.3)
-  )
-  expect_equal(pval(read_param(pth), "plantStorageNInit"), 3.3)
-
-  # zero initial wood means no flush to support: storage N stays 0.
-  # the computation must use the IC adjusted wood pool, not the template's.
-  pth <- withr::local_tempdir()
-  write.config.SIPNET(
-    defaults = list(pft1 = list(constants = list(SLA = 2.0))),
-    trait.values = list(pft1 = list()),
-    settings = make_settings(pth), run.id = "run1",
-    IC = data.frame(AbvGrndWood = 0, abvGrndWoodFrac = 0.6,
-                    coarseRootFrac = 0.2, fineRootFrac = 0.2)
-  )
-  p <- read_param(pth)
-  expect_equal(pval(p, "plantWoodInit"), 0)
-  expect_equal(pval(p, "plantStorageNInit"), 0)
-})