write n2o and ch4 flux variables to SIPNET netCDF

models/sipnet/R/model2netcdf.SIPNET.R

@@ -178,6 +178,15 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
     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
+    # columns only present in sipnet >= v2 with N and methane turned on
+    if ("n2o" %in% names(sub.sipnet.output)) {
+      output[["N2O_flux"]] <- (sub.sipnet.output$n2o * 0.001) / timestep.s
+      # convert g N m-2 per timestep -> kg N m-2 s-1
+    }
+    if ("ch4" %in% names(sub.sipnet.output)) {
+      output[["CH4_flux"]] <- (sub.sipnet.output$ch4 * 0.001) / timestep.s
+      # convert g C m-2 per timestep -> kg C m-2 s-1
+    }
     output[["time_bounds"]] <- c(rbind(bounds[,1], bounds[,2]))
 
     # ******************** Declare netCDF variables ********************#
@@ -236,6 +245,13 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
                                        longname = "history time interval endpoints", dim=list(time_interval,time = t), 
                                        prec = "double")              
     )
+
+    if ("N2O_flux" %in% names(output)) {
+      nc_var[["N2O_flux"]] <- PEcAn.utils::to_ncvar("N2O_flux", dims)
+    }
+    if ("CH4_flux" %in% names(output)) {
+      nc_var[["CH4_flux"]] <- PEcAn.utils::to_ncvar("CH4_flux", dims)
+    }
 
     # ******************** Create netCDF and output variables ********************#
     ### Output netCDF data
@@ -264,8 +280,8 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
     }else{
       nc      <- ncdf4::nc_create(file.path(outdir, paste(y, "nc", sep = ".")), nc_var)
       ncdf4::ncatt_put(nc, "time", "bounds", "time_bounds", prec=NA)
-      for (i in seq_along(nc_var)) {
-        ncdf4::ncvar_put(nc, nc_var[[i]], output[[i]])
+      for (key in names(nc_var)) {
+        ncdf4::ncvar_put(nc, nc_var[[key]], output[[key]])
       }
       ncdf4::nc_close(nc)
     }

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

@@ -0,0 +1,98 @@
+setup_sipnet_test <- function(sipnet_dat, delete.raw = FALSE) {
+  base <- withr::local_tempdir(pattern = "sipnet_test_", .local_envir = parent.frame())
+  outdir <- file.path(base, "out", "run1")
+  rundir <- file.path(base, "run", "run1")
+  dir.create(outdir, recursive = TRUE)
+  dir.create(rundir, recursive = TRUE)
+
+  writeLines("leafCSpWt\t32", file.path(rundir, "sipnet.param"))
+
+  out_path <- file.path(outdir, "sipnet.out")
+  writeLines("Notes: units in g/m2 per timestep; water in cm", out_path)
+  write.table(sipnet_dat, file = out_path, append = TRUE,
+              row.names = FALSE, quote = FALSE, sep = "\t")
+
+
+  model2netcdf.SIPNET(
+    outdir     = outdir,
+    sitelat    = 38.0,
+    sitelon    = -121.0,
+    start_date = "2002-01-01",
+    end_date   = "2002-12-31",
+    delete.raw = delete.raw,
+    revision   = "r136"
+  )
+
+  list(outdir = outdir, rundir = rundir, out_path = out_path)
+}
+
+make_base_sipnet <- function(n = 4L) {
+  data.frame(
+    year = 2002,
+    day = rep(c(1, 2), each = n / 2, length.out = n),
+    time = rep(c(6, 18), length.out = n),
+    plantWoodC = 5000, plantLeafC = 200, woodCreation = 0.5,
+    soil = 10000, microbeC = 8, coarseRootC = 1200, fineRootC = 800,
+    litter = 400, soilWater = 14, soilWetnessFrac = 0.85, snow = 0,
+    npp = 0.05, nee = 0.10, cumNEE = cumsum(rep(0.1, n)),
+    gpp = 0.30, rAboveground = 0.04, rSoil = 0.09, rRoot = 0.01,
+    ra = 0.05, rh = 0.08, rtot = 0.13,
+    evapotranspiration = 0.005, fluxestranspiration = 0.003
+  )
+}
+
+
+test_that("model2netcdf.SIPNET converts v2 output including N2O and CH4 fluxes", {
+  n <- 4L
+  ts_s <- 43200
+  sipnet_dat <- make_base_sipnet(n)
+  sipnet_dat$n2o <- 0.002
+  sipnet_dat$ch4 <- 0.001
+  paths <- setup_sipnet_test(sipnet_dat)
+  nc_file <- file.path(paths$outdir, "2002.nc")
+  expect_true(file.exists(nc_file))
+
+  nc <- ncdf4::nc_open(nc_file)
+  on.exit(ncdf4::nc_close(nc), add = TRUE)
+  vars <- names(nc$var)
+
+  expect_true("N2O_flux" %in% vars)
+  expect_true("CH4_flux" %in% vars)
+  expect_true(all(c("GPP", "NEE", "TotalResp", "TotSoilCarb") %in% vars))
+
+  n2o <- as.vector(ncdf4::ncvar_get(nc, "N2O_flux"))
+  ch4 <- as.vector(ncdf4::ncvar_get(nc, "CH4_flux"))
+  gpp <- as.vector(ncdf4::ncvar_get(nc, "GPP"))
+
+  expect_equal(n2o, rep(0.002 * 1e-3 / ts_s, n), tolerance = 1e-12)
+  expect_equal(ch4, rep(0.001 * 1e-3 / ts_s, n), tolerance = 1e-12)
+  expect_equal(gpp, rep(0.30  * 1e-3 / ts_s, n), tolerance = 1e-12)
+
+  expect_equal(nc$var$N2O_flux$units, "kg N m-2 s-1")
+  expect_equal(nc$var$CH4_flux$units, "kg C m-2 s-1")
+  expect_equal(nc$var$GPP$units,      "kg C m-2 s-1")
+
+  expect_equal(nc$dim$time$len, n)
+  expect_match(nc$dim$time$units, "days since 2002")
+})
+
+
+test_that("model2netcdf.SIPNET omits N2O/CH4 when columns absent", {
+  paths <- setup_sipnet_test(make_base_sipnet(n = 2L))
+
+  nc <- ncdf4::nc_open(file.path(paths$outdir, "2002.nc"))
+  on.exit(ncdf4::nc_close(nc), add = TRUE)
+  vars <- names(nc$var)
+
+  expect_false("N2O_flux" %in% vars)
+  expect_false("CH4_flux" %in% vars)
+  expect_true("GPP" %in% vars)
+})
+
+
+test_that("delete.raw removes sipnet.out after conversion", {
+  paths <- setup_sipnet_test(make_base_sipnet(n = 2L), delete.raw = TRUE)
+
+  expect_false(file.exists(paths$out_path))
+  expect_true(file.exists(file.path(paths$outdir, "2002.nc")))
+})