Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
30 changes: 15 additions & 15 deletions models/dalec/R/model2netcdf.DALEC.R
Original file line number Diff line number Diff line change
Expand Up @@ -83,22 +83,22 @@ model2netcdf.DALEC <- function(outdir, sitelat, sitelon, start_date, end_date) {

## Setup outputs for netCDF file in appropriate units
output <- list()
## Fluxes
output[[1]] <- (sub.DALEC.output[, 1] * 0.001)/timestep.s # Autotrophic Respiration in kgC/m2/s
output[[2]] <- (sub.DALEC.output[, 21] + sub.DALEC.output[, 23]) * 0.001 / timestep.s # Heterotrophic Resp kgC/m2/s
output[[3]] <- (sub.DALEC.output[, 31] * 0.001)/timestep.s # GPP in kgC/m2/s
output[[4]] <- (sub.DALEC.output[, 33] * 0.001)/timestep.s # NEE in kgC/m2/s
output[[5]] <- (sub.DALEC.output[, 3] + sub.DALEC.output[, 5] + sub.DALEC.output[, 7]) * 0.001/timestep.s # NPP kgC/m2/s
output[[6]] <- (sub.DALEC.output[, 9] * 0.001) / timestep.s # Leaf Litter Flux, kgC/m2/s
output[[7]] <- (sub.DALEC.output[, 11] * 0.001) / timestep.s # Woody Litter Flux, kgC/m2/s
output[[8]] <- (sub.DALEC.output[, 13] * 0.001) / timestep.s # Root Litter Flux, kgC/m2/s
## Fluxes (convert g/m2/day to kg/m2/s)
output[[1]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 1], "g/m2/d", "kg/m2/s") # Autotrophic Respiration
output[[2]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 21] + sub.DALEC.output[, 23], "g/m2/d", "kg/m2/s") # Heterotrophic Resp
output[[3]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 31], "g/m2/d", "kg/m2/s") # GPP
output[[4]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 33], "g/m2/d", "kg/m2/s") # NEE
output[[5]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 3] + sub.DALEC.output[, 5] + sub.DALEC.output[, 7], "g/m2/d", "kg/m2/s") # NPP
output[[6]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 9], "g/m2/d", "kg/m2/s") # Leaf Litter Flux
output[[7]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 11], "g/m2/d", "kg/m2/s") # Woody Litter Flux
output[[8]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 13], "g/m2/d", "kg/m2/s") # Root Litter Flux

## Pools
output[[9]] <- (sub.DALEC.output[, 15] * 0.001) # Leaf Carbon, kgC/m2
output[[10]] <- (sub.DALEC.output[, 17] * 0.001) # Wood Carbon, kgC/m2
output[[11]] <- (sub.DALEC.output[, 19] * 0.001) # Root Carbon, kgC/m2
output[[12]] <- (sub.DALEC.output[, 27] * 0.001) # Litter Carbon, kgC/m2
output[[13]] <- (sub.DALEC.output[, 29] * 0.001) # Soil Carbon, kgC/m2
## Pools (convert g/m2 to kg/m2)
output[[9]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 15], "g/m2", "kg/m2") # Leaf Carbon
output[[10]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 17], "g/m2", "kg/m2") # Wood Carbon
output[[11]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 19], "g/m2", "kg/m2") # Root Carbon
output[[12]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 27], "g/m2", "kg/m2") # Litter Carbon
output[[13]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 29], "g/m2", "kg/m2") # Soil Carbon

## standard composites
output[[14]] <- output[[1]] + output[[2]] # Total Respiration
Expand Down
8 changes: 4 additions & 4 deletions models/fates/R/write.configs.FATES.R
Original file line number Diff line number Diff line change
Expand Up @@ -307,25 +307,25 @@ write.config.FATES <- function(defaults, trait.values, settings, run.id){
# Ha activation energy for vcmax - FATES units: J/mol
if(var == "Ha_Modified_Arrhenius_Vcmax"){
ncdf4::ncvar_put(nc=fates.param.nc, varid='fates_vcmaxha', start = ipft, count = 1,
vals=pft[v]*1000) ## convert from kj/mol to J/mol (FATES units)
vals=PEcAn.utils::ud_convert(pft[v], "kJ/mol", "J/mol")) ## convert from kJ/mol to J/mol (FATES units)
}

# Hd deactivation energy for vcmax - FATES units: J/mol
if(var == "Hd_Modified_Arrhenius_Vcmax"){
ncdf4::ncvar_put(nc=fates.param.nc, varid='fates_vcmaxhd', start = ipft, count = 1,
vals=pft[v]*1000) ## convert from kj/mol to J/mol (FATES units)
vals=PEcAn.utils::ud_convert(pft[v], "kJ/mol", "J/mol")) ## convert from kJ/mol to J/mol (FATES units)
}

# Ha activation energy for Jmax - FATES units: J/mol
if(var == "Ha_Modified_Arrhenius_Jmax"){
ncdf4::ncvar_put(nc=fates.param.nc, varid='fates_jmaxha', start = ipft, count = 1,
vals=pft[v]*1000) ## convert from kj/mol to J/mol (FATES units)
vals=PEcAn.utils::ud_convert(pft[v], "kJ/mol", "J/mol")) ## convert from kJ/mol to J/mol (FATES units)
}

# Hd deactivation energy for Jmax - FATES units: J/mol
if(var == "Hd_Modified_Arrhenius_Jmax"){
ncdf4::ncvar_put(nc=fates.param.nc, varid='fates_jmaxhd', start = ipft, count = 1,
vals=pft[v]*1000) ## convert from kj/mol to J/mol (FATES units)
vals=PEcAn.utils::ud_convert(pft[v], "kJ/mol", "J/mol")) ## convert from kJ/mol to J/mol (FATES units)
}

# deltaS Vcmax - BETY units:J/mol/K; FATES units: J/mol/K
Expand Down
23 changes: 10 additions & 13 deletions models/gday/R/model2netcdf.GDAY.R
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,6 @@
model2netcdf.GDAY <- function(outdir, sitelat, sitelon, start_date, end_date) {


G_2_KG <- 0.001
TONNES_PER_HA_TO_G_M2 <- 100
THA_2_KG_M2 <- TONNES_PER_HA_TO_G_M2 * 0.001

### Read in model output in GDAY format
GDAY.output <- utils::read.csv(file.path(outdir, "gday_out.csv"), header = TRUE, sep = ",", skip = 1)
Expand All @@ -44,17 +41,17 @@ model2netcdf.GDAY <- function(outdir, sitelat, sitelon, start_date, end_date) {
output <- list()

## standard variables: C-Fluxes
output[[1]] <- (sub.GDAY.output[, "auto_resp"] * THA_2_KG_M2) / timestep.s
output[[2]] <- (sub.GDAY.output[, "hetero_resp"] * THA_2_KG_M2) / timestep.s
output[[3]] <- (sub.GDAY.output[, "auto_resp"] + sub.GDAY.output[, "hetero_resp"] *
THA_2_KG_M2) / timestep.s
output[[4]] <- (sub.GDAY.output[, "gpp"] * THA_2_KG_M2) / timestep.s
output[[5]] <- (sub.GDAY.output[, "nep"] * -1 * THA_2_KG_M2) / timestep.s
output[[6]] <- (sub.GDAY.output[, "npp"] * THA_2_KG_M2) / timestep.s
# GDAY outputs in Mg/ha/year, convert directly to kgC/m2/s (ud_convert returns per-second)
output[[1]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "auto_resp"], "Mg/ha/yr", "kg/m2/s")
output[[2]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "hetero_resp"], "Mg/ha/yr", "kg/m2/s")
output[[3]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "auto_resp"] + sub.GDAY.output[, "hetero_resp"], "Mg/ha/yr", "kg/m2/s")
output[[4]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "gpp"], "Mg/ha/yr", "kg/m2/s")
output[[5]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "nep"] * -1, "Mg/ha/yr", "kg/m2/s")
output[[6]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "npp"], "Mg/ha/yr", "kg/m2/s")

## standard variables: C-State
output[[7]] <- (sub.GDAY.output[, "stem"] + sub.GDAY.output[, "branch"] * THA_2_KG_M2) / timestep.s
output[[8]] <- (sub.GDAY.output[, "soilc"] * THA_2_KG_M2) / timestep.s
## standard variables: C-State (pools) - divide by timestep as in original code
output[[7]] <- (PEcAn.utils::ud_convert(sub.GDAY.output[, "stem"], "Mg/ha", "kg/m2") + sub.GDAY.output[, "branch"] * THA_2_KG_M2) / timestep.s
output[[8]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "soilc"], "Mg/ha", "kg/m2") / timestep.s
output[[9]] <- (sub.GDAY.output[, "lai"])

## standard variables: water fluxes
Expand Down
41 changes: 20 additions & 21 deletions models/sipnet/R/model2netcdf.SIPNET.R
Original file line number Diff line number Diff line change
Expand Up @@ -134,21 +134,20 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
bounds <- round(bounds,4)

## Setup outputs for netCDF file in appropriate units
output <- list(
"GPP" = (sub.sipnet.output$gpp * 0.001) / timestep.s, # GPP in kgC/m2/s
"NPP" = (sub.sipnet.output$gpp * 0.001) / timestep.s - ((sub.sipnet.output$rAboveground *
0.001) / timestep.s + (sub.sipnet.output$rRoot * 0.001) / timestep.s), # NPP in kgC/m2/s. Post SIPNET calculation
"TotalResp" = (sub.sipnet.output$rtot * 0.001) / timestep.s, # Total Respiration in kgC/m2/s
"AutoResp" = (sub.sipnet.output$rAboveground * 0.001) / timestep.s + (sub.sipnet.output$rRoot *
0.001) / timestep.s, # Autotrophic Respiration in kgC/m2/s
"HeteroResp" = ((sub.sipnet.output$rSoil - sub.sipnet.output$rRoot) * 0.001) / timestep.s, # Heterotrophic Respiration in kgC/m2/s
"SoilResp" = (sub.sipnet.output$rSoil * 0.001) / timestep.s, # Soil Respiration in kgC/m2/s
"NEE" = (sub.sipnet.output$nee * 0.001) / timestep.s, # NEE in kgC/m2/s
"AbvGrndWood" = (sub.sipnet.output$plantWoodC * 0.001), # Above ground wood kgC/m2
"leaf_carbon_content" = (sub.sipnet.output$plantLeafC * 0.001), # Leaf C kgC/m2
"TotLivBiom" = (sub.sipnet.output$plantWoodC * 0.001) + (sub.sipnet.output$plantLeafC * 0.001) +
(sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC) * 0.001, # Total living C kgC/m2
"TotSoilCarb" = (sub.sipnet.output$soil * 0.001) + (sub.sipnet.output$litter * 0.001) # Total soil C kgC/m2
output <- list(
"GPP" = PEcAn.utils::ud_convert(sub.sipnet.output$gpp, "g/m2", "kg/m2") / timestep.s,
"NPP" = (PEcAn.utils::ud_convert(sub.sipnet.output$gpp, "g/m2", "kg/m2") -
PEcAn.utils::ud_convert(sub.sipnet.output$rAboveground + sub.sipnet.output$rRoot, "g/m2", "kg/m2")) / timestep.s,
"TotalResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rtot, "g/m2", "kg/m2") / timestep.s,
"AutoResp" = (PEcAn.utils::ud_convert(sub.sipnet.output$rAboveground + sub.sipnet.output$rRoot, "g/m2", "kg/m2")) / timestep.s,
"HeteroResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rSoil - sub.sipnet.output$rRoot, "g/m2", "kg/m2") / timestep.s,
"SoilResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rSoil, "g/m2", "kg/m2") / timestep.s,
"NEE" = PEcAn.utils::ud_convert(sub.sipnet.output$nee, "g/m2", "kg/m2") / timestep.s,
"AbvGrndWood" = PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC, "g/m2", "kg/m2"),
"leaf_carbon_content" = PEcAn.utils::ud_convert(sub.sipnet.output$plantLeafC, "g/m2", "kg/m2"),
"TotLivBiom" = (PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC +
sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC, "g/m2", "kg/m2")),
"TotSoilCarb" = PEcAn.utils::ud_convert(sub.sipnet.output$soil + sub.sipnet.output$litter, "g/m2", "kg/m2")
)
if (revision == "unk") {
## *** NOTE : npp in the sipnet output file is actually evapotranspiration, this is due to a bug in sipnet.c : ***
Expand All @@ -164,8 +163,8 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
output[["SoilMoist"]] <- (sub.sipnet.output$soilWater * 10) # Soil moisture kgW/m2
output[["SoilMoistFrac"]] <- (sub.sipnet.output$soilWetnessFrac) # Fractional soil wetness
output[["SWE"]] <- (sub.sipnet.output$snow * 10) # SWE
output[["litter_carbon_content"]] <- sub.sipnet.output$litter * 0.001 ## litter kgC/m2
output[["litter_mass_content_of_water"]] <- (sub.sipnet.output$litterWater * 10) # Litter water kgW/m2
output[["litter_carbon_content"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$litter, "g/m2", "kg/m2")
output[["litter_mass_content_of_water"]] <- (sub.sipnet.output$litterWater * 10) # * 10 converts cm to mm
#calculate LAI for standard output
param <- utils::read.table(file.path(gsub(pattern = "/out/",
replacement = "/run/", x = outdir),
Expand All @@ -174,10 +173,10 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
leafC <- 0.48
SLA <- 1000 / param[id, 2] #SLA, m2/kgC
output[["LAI"]] <- output[["leaf_carbon_content"]] * SLA # LAI
output[["fine_root_carbon_content"]] <- sub.sipnet.output$fineRootC * 0.001 ## fine_root_carbon_content kgC/m2
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
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")
output[["GWBI"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$woodCreation, "g/m2/d", "kg/m2/s")
output[["AGB"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC, "g/m2", "kg/m2")
output[["time_bounds"]] <- c(rbind(bounds[,1], bounds[,2]))

# ******************** Declare netCDF variables ********************#
Expand Down
3 changes: 1 addition & 2 deletions models/sipnet/R/write_restart.SIPNET.R
Original file line number Diff line number Diff line change
Expand Up @@ -47,12 +47,11 @@ write_restart.SIPNET <- function(outdir, runid, start.time, stop.time, settings,

## Converting to sipnet units
prior.sla <- new.params[[which(!names(new.params) %in% c("soil", "soil_SDA", "restart"))[1]]]$SLA
unit.conv <- 2 * (10000 / 1) * (1 / 1000) * (3.154 * 10^7) # kgC/m2/s -> Mg/ha/yr

analysis.save <- list()

if ("NPP" %in% variables) {
analysis.save[[length(analysis.save) + 1]] <- PEcAn.utils::ud_convert(new.state$NPP, "kg/m^2/s", "Mg/ha/yr") #*unit.conv -> Mg/ha/yr
analysis.save[[length(analysis.save) + 1]] <- PEcAn.utils::ud_convert(new.state$NPP, "kg/m^2/s", "Mg/ha/yr")
names(analysis.save[[length(analysis.save)]]) <- c("NPP")
}

Expand Down
2 changes: 1 addition & 1 deletion models/stics/R/met2model.STICS.R
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ met2model.STICS <- function(in.path, in.prefix, outfolder, start_date, end_date,
# column 10: rainfall (mm.j-1)
Rain <- ncdf4::ncvar_get(nc, "precipitation_flux") # kg m-2 s-1
Rain <- Rain[ydays %in% simdays]
raini <- tapply(Rain * 86400, ind, mean, na.rm = TRUE)
raini <- tapply(PEcAn.utils::ud_convert(Rain, "kg/m2/s", "kg/m2/d"), ind, mean, na.rm = TRUE) # 1 kg/m2 = 1 mm water
weather_df[ ,10] <- round(raini, digits = 2) # precipitation (mm d-1)

# column 11: wind (m.s-1)
Expand Down
49 changes: 49 additions & 0 deletions tests/test_unit_conversions.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,49 @@
# Test script to validate unit conversion changes
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for including tests! Unit tests like this should be placed in the tests/testthat folder of the appropriate package and not to the top-level tests folder, which is used for whole-system test material (and is due for a reorg/cleanout, but that's a separate issue). For tests of ud_convert itself, the appropriate package is PEcAn.utils so they should be added to the existing base/utils/tests/testthat/test-ud_convert.R.

I also request using standard testthat expectations rather than rolling your own comparison machinery -- this entire script is basically equivalent to

expect_equal(ud_convert(100, "g/m2", "kg/m2"), 0.1) # eg DALEC/SIPNET C pools
expect_equal(ud_convert(86400, "g/m2/d", "kg/m2/s"), 0.001) # eg DALEC/SIPNET C fluxes
expect_equal(ud_convert(10, "Mg/ha", "kg/m2"), 1) # eg GDAY pools
expect_equal(ud_convert(1000, "J/mol", "kJ/mol"), 1) # eg photosynthesis params

# This verifies that ud_convert works with the new unit strings

library(PEcAn.utils)

cat("Testing unit string validity with UDUNITS2...\n\n")

# Test cases from the modified files
test_conversions <- list(
list(value = 100, u1 = "g/m2", u2 = "kg/m2", expected = 0.1, desc = "DALEC/SIPNET pool: g/m2 to kg/m2"),
list(value = 86400, u1 = "g/m2/d", u2 = "kg/m2/s", expected = 0.001, desc = "DALEC/SIPNET flux: g/m2/d to kg/m2/s"),
list(value = 10, u1 = "Mg/ha", u2 = "kg/m2", expected = 1, desc = "GDAY: Mg/ha to kg/m2"),
list(value = 1000, u1 = "J/mol", u2 = "kJ/mol", expected = 1, desc = "FATES: J/mol to kJ/mol")
)

results <- list()
for (i in seq_along(test_conversions)) {
test <- test_conversions[[i]]
tryCatch({
result <- ud_convert(test$value, test$u1, test$u2)
passed <- abs(result - test$expected) < 1e-10
results[[i]] <- list(
desc = test$desc,
passed = passed,
result = result,
expected = test$expected
)
cat(sprintf("[%s] %s\n", if(passed) "PASS" else "FAIL", test$desc))
if (!passed) {
cat(sprintf(" Expected: %f, Got: %f\n", test$expected, result))
}
}, error = function(e) {
results[[i]] <<- list(desc = test$desc, passed = FALSE, error = e$message)
cat(sprintf("[FAIL] %s\n", test$desc))
cat(sprintf(" Error: %s\n", e$message))
})
}

cat("\n")
passed_count <- sum(sapply(results, function(x) x$passed))
total_count <- length(results)
cat(sprintf("Summary: %d/%d tests passed\n", passed_count, total_count))

if (passed_count < total_count) {
cat("\nFAILED TESTS DETECTED - These must be fixed before proceeding\n")
quit(status = 1)
} else {
cat("\nAll unit strings are valid!\n")
}
4 changes: 4 additions & 0 deletions tests/testthat.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
# Regression tests for unit conversion refactoring
# Run with: testthat::test_dir("tests/testthat", filter = "test_model_conversions")

source("test_model_conversions.R")
Comment on lines +1 to +4
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

delete this file -- test_model_conversions.R belongs in base/utils/tests/testthat, which already has an existing testthat.R

66 changes: 66 additions & 0 deletions tests/testthat/test_model_conversions.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@

# These tests validate the refactoring to use PEcAn.utils::ud_convert()
# for consistent unit conversions across SIPNET, DALEC, GDAY, and FATES models.

test_that("UDUNITS2 recognizes new unit strings", {

# Pool conversions (mass only)
expect_equal(ud_convert(100, "g/m2", "kg/m2"), 0.1)
expect_equal(ud_convert(1000, "g/m2", "kg/m2"), 1.0)
Comment on lines +7 to +9
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do not think we need to check the same unit conversion at multiple values -- Most of the actual conversion machinery is provided by the external units package and we trust them to test their own implementation, so our goal here should be to verify support for units that are important to us rather than to exhaustively test specific input values.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or in less words: "Let's only check each unit pair once"


# Flux conversions (mass and time)
expect_equal(ud_convert(86400, "g/m2/d", "kg/m2/s"), 0.001, tolerance = 1e-10)

# Large-scale conversions (GDAY)
expect_equal(ud_convert(10, "Mg/ha", "kg/m2"), 1.0)

# Energy conversions (FATES)
expect_equal(ud_convert(1000, "J/mol", "kJ/mol"), 1.0)
})

test_that("SIPNET flux conversions produce realistic values", {

# GPP: ~5 g C/m2/d --> kg/m2/s
expect_equal(ud_convert(5, "g/m2/d", "kg/m2/s"), 5.79e-8)

# Pool: 1000 g C/m2 → 1 kg/m2
wood_result <- ud_convert(1000, "g/m2", "kg/m2")
expect_equal(wood_result, 1.0)
})

test_that("DALEC flux conversions are mathematically sound", {
library(PEcAn.utils)

# Autotrophic respiration: 5 g C/m2/d → kg/m2/s
ar_result <- ud_convert(5, "g/m2/d", "kg/m2/s")
expect_true(ar_result > 0)
expect_true(ar_result < 5)

# Leaf carbon: 200 g C/m2 → 0.2 kg/m2
leaf_result <- ud_convert(200, "g/m2", "kg/m2")
expect_equal(leaf_result, 0.2)
})

test_that("GDAY conversions from Mg/ha to kg/m2 work correctly", {
library(PEcAn.utils)

# 5 Mg/ha = 0.5 kg/m2
stem_result <- ud_convert(5, "Mg/ha", "kg/m2")
expect_equal(stem_result, 0.5)

# 50 Mg/ha = 5 kg/m2
soil_result <- ud_convert(50, "Mg/ha", "kg/m2")
expect_equal(soil_result, 5.0)
})

test_that("No 'C' prefix in unit strings (UDUNITS2 requirement)", {
library(PEcAn.utils)

# These should NOT work (they include 'C' for carbon)
expect_error(ud_convert(100, "gC/m2", "kgC/m2"))
expect_error(ud_convert(86400, "gC/m2/d", "kgC/m2/s"))
Comment on lines +60 to +61
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like that you're taking a problem seen during development and trying to turn it into a test! But I don't think these expectations are useful to us -- they will only fail if udunits does learn how to parse C, and alas having the test case here won't keep users from making the same typo elsewhere.

It's possible that C is a common enough special case that it'd be nice for PEcAn's ud_convert to handle it explicitly (either converting it correctly or throwing a more informative error), but I suspect that would be a lot of tricky work for little practical gain.


# But these (without 'C') should work
expect_equal(ud_convert(100, "g/m2", "kg/m2"), 0.1)
expect_equal(ud_convert(86400, "g/m2/d", "kg/m2/s"), 0.001, tolerance = 1e-10)
Comment on lines +63 to +65
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

already checked above

})
Loading