read LAI from sipnet.out if present

models/sipnet/NEWS.md

@@ -1,5 +1,8 @@
 # PEcAn.SIPNET 1.10.0.9000
 
+* `model2netcdf.SIPNET` takes `LAI` from sipnet.out if present (which requires
+    Sipnet > v2.2). If it is not present, LAI is calculated as
+    `plantLeafC / leafCSpWt` as previously.
 * `split_inputs.SIPNET` now avoids internal time format conversions, giving a
   substantial speedup and reduced memory use when processing multi-year files.
 * `model2netcdf.SIPNET` now detects the number of timesteps per day by taking the maximum count across all days in the first simulation year, rather than reading only from day 1. This prevents a factor-of-N error in flux unit conversions when the first day of output is partial (fewer timesteps than a complete day) (#3624, #3989).

models/sipnet/R/model2netcdf.SIPNET.R

@@ -150,16 +150,21 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
     )
 
 
-  # calculate LAI for standard output
-  # LAI = plantLeafC / leafCSpWt
-  # both operands are in carbon units (gC/m2 and gC/m2_leaf),
-  # so no carbon fraction conversion (e.g. cFracLeaf) is needed.
-  param <- utils::read.table(file.path(gsub(pattern = "/out/",
-                                            replacement = "/run/", x = outdir),
-                                       "sipnet.param"), stringsAsFactors = FALSE)
-  leafCSpWt <- param[param[, 1] == "leafCSpWt", 2]
-  SLA <- 1000 / leafCSpWt  # m2 leaf / kg C
-
+  if (!("LAI" %in% colnames(sipnet_output))) {
+    # When LAI not reported, calculate it from leaf C and leaf specific weight,
+    # with the latter read from the parameter file.
+    # Note the hardcoded + undocumented assumption that ../run/sipnet.param
+    # exists and contains the values that were used to generate this output.
+    # LAI = plantLeafC [gC/m2] / leafCSpWt [gC/m2 leaf] = m2 leaf / m2
+    # both operands are in carbon units (gC/m2 and gC/m2_leaf),
+    # so no carbon fraction conversion (e.g. cFracLeaf) is needed.
+    param <- utils::read.table(file.path(gsub(pattern = "/out/",
+                                              replacement = "/run/", x = outdir),
+                                         "sipnet.param"), stringsAsFactors = FALSE)
+    leafCSpWt <- param[param[, 1] == "leafCSpWt", 2]
+    SLA <- 1000 / leafCSpWt  # m2 leaf / kg C
+    sipnet_output$LAI = sipnet_output$plantLeafC * SLA
+  }
 
   ### Loop over years in SIPNET output to create separate netCDF outputs
   for (y in year_seq) {
@@ -208,7 +213,7 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
       "litter_carbon_content" = sub.sipnet.output$litter,
       "fine_root_carbon_content" = sub.sipnet.output$fineRootC,
       "coarse_root_carbon_content" = sub.sipnet.output$coarseRootC,
-      "LAI" = sub.sipnet.output$plantLeafC * SLA,
+      "LAI" = sub.sipnet.output$LAI,
       "TotLivBiom" = sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC +
                        sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC,
       "TotSoilCarb" = sub.sipnet.output$soil + sub.sipnet.output$litter,

models/sipnet/tests/testthat/test-model2netcdf.SIPNET.R

@@ -291,3 +291,22 @@ test_that("out_day is correct when the first day of output is partial", {
   ts_correct <- 86400 / 2
   expect_equal(gpp, rep(base_row$gpp * 1e-3 / ts_correct, 5), tolerance = 1e-10)
 })
+
+test_that("LAI taken from leafCSpWt if not present in output", {
+  dat <- make_v2_sipnet(2)
+  # setup_sipnet_test sets leafCSpWt = 32 => expect LAI of 1 and 10
+  dat$plantLeafC = c(32, 320)
+
+  # First without LAI in dat => looks for sipnet.param in run dir
+  paths <- setup_sipnet_test(dat, notes_line = NULL)
+  nc1 <- ncdf4::nc_open(file.path(paths$outdir, "2002.nc"))
+  on.exit(ncdf4::nc_close(nc1), add = TRUE)
+  expect_equal(as.vector(ncdf4::ncvar_get(nc1, "LAI")), c(1, 10))
+
+  # Now with LAI specified => uses it as-is
+  dat$LAI = c(2, 3)
+  paths <- setup_sipnet_test(dat, notes_line = NULL)
+  nc2 <- ncdf4::nc_open(file.path(paths$outdir, "2002.nc"))
+  on.exit(ncdf4::nc_close(nc2), add = TRUE)
+  expect_equal(as.vector(ncdf4::ncvar_get(nc2, "LAI")), c(2, 3))
+})