Replace manual unit conversions with ud_convert() for consistency across models

models/dalec/R/model2netcdf.DALEC.R

@@ -84,21 +84,23 @@ 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
+    # Inputs are in gC/m2/day; ud_convert(..., "gC/m2/d", "kgC/m2/s") returns kgC/m2/s
+    # Do NOT divide again by timestep.s (that would double-scale time)
+    output[[1]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 1], "gC/m2/d", "kgC/m2/s")  # Autotrophic Respiration in kgC/m2/s
+    output[[2]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 21] + sub.DALEC.output[, 23], "gC/m2/d", "kgC/m2/s")  # Heterotrophic Resp kgC/m2/s
+    output[[3]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 31], "gC/m2/d", "kgC/m2/s")  # GPP in kgC/m2/s    
+    output[[4]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 33], "gC/m2/d", "kgC/m2/s")  # NEE in kgC/m2/s
+    output[[5]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 3] + sub.DALEC.output[, 5] + sub.DALEC.output[, 7], "gC/m2/d", "kgC/m2/s")  # NPP kgC/m2/s
+    output[[6]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 9], "gC/m2/d", "kgC/m2/s")  # Leaf Litter Flux, kgC/m2/s
+    output[[7]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 11], "gC/m2/d", "kgC/m2/s")  # Woody Litter Flux, kgC/m2/s
+    output[[8]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 13], "gC/m2/d", "kgC/m2/s")  # Root Litter Flux, kgC/m2/s
 
     ## 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
+    output[[9]]  <- PEcAn.utils::ud_convert(sub.DALEC.output[, 15], "gC/m2", "kgC/m2")  # Leaf Carbon, kgC/m2
+    output[[10]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 17], "gC/m2", "kgC/m2")  # Wood Carbon, kgC/m2
+    output[[11]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 19], "gC/m2", "kgC/m2")  # Root Carbon, kgC/m2
+    output[[12]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 27], "gC/m2", "kgC/m2")  # Litter Carbon, kgC/m2
+    output[[13]] <- PEcAn.utils::ud_convert(sub.DALEC.output[, 29], "gC/m2", "kgC/m2")  # Soil Carbon, kgC/m2
 
     ## standard composites
     output[[14]] <- output[[1]] + output[[2]]  # Total Respiration

models/fates/R/write.configs.FATES.R

@@ -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

models/gday/R/model2netcdf.GDAY.R

@@ -15,9 +15,9 @@
 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
+# Unit conversions are handled via PEcAn.utils::ud_convert()
+# except for water fluxes, which use model-specific daily-to-second scaling
+
 
   ### Read in model output in GDAY format
   GDAY.output <- utils::read.csv(file.path(outdir, "gday_out.csv"), header = TRUE, sep = ",", skip = 1)
@@ -44,17 +44,18 @@ 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
+    ## C-State (pools) - convert Mg/ha to kg/m2 (no time division)
+    output[[7]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "stem"] + sub.GDAY.output[, "branch"], "Mg/ha", "kg/m2")
+    output[[8]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "soilc"], "Mg/ha", "kg/m2")
     output[[9]] <- (sub.GDAY.output[, "lai"])
 
     ## standard variables: water fluxes

models/sipnet/R/model2netcdf.SIPNET.R

@@ -135,20 +135,18 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
 
     ## 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
+      "GPP" = PEcAn.utils::ud_convert(sub.sipnet.output$gpp, "gC/m2", "kgC/m2") / timestep.s,  # GPP in kgC/m2/s
+      "NPP" = PEcAn.utils::ud_convert(sub.sipnet.output$gpp, "gC/m2", "kgC/m2") / timestep.s - (PEcAn.utils::ud_convert(sub.sipnet.output$rAboveground, "gC/m2", "kgC/m2") / timestep.s + PEcAn.utils::ud_convert(sub.sipnet.output$rRoot, "gC/m2", "kgC/m2") / timestep.s), # NPP in kgC/m2/s. Post SIPNET calculation
+      "TotalResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rtot, "gC/m2", "kgC/m2") / timestep.s,  # Total Respiration in kgC/m2/s
+      "AutoResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rAboveground, "gC/m2", "kgC/m2") / timestep.s + PEcAn.utils::ud_convert(sub.sipnet.output$rRoot, "gC/m2", "kgC/m2") / timestep.s,  # Autotrophic Respiration in kgC/m2/s
+      "HeteroResp" = PEcAn.utils::ud_convert((sub.sipnet.output$rSoil - sub.sipnet.output$rRoot), "gC/m2", "kgC/m2") / timestep.s,  # Heterotrophic Respiration in kgC/m2/s
+      "SoilResp" = PEcAn.utils::ud_convert(sub.sipnet.output$rSoil, "gC/m2", "kgC/m2") / timestep.s,  # Soil Respiration in kgC/m2/s
+      "NEE" = PEcAn.utils::ud_convert(sub.sipnet.output$nee, "gC/m2", "kgC/m2") / timestep.s,  # NEE in kgC/m2/s
+      "AbvGrndWood" = PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC, "gC/m2", "kgC/m2"),  # Above ground wood kgC/m2
+      "leaf_carbon_content" = PEcAn.utils::ud_convert(sub.sipnet.output$plantLeafC, "gC/m2", "kgC/m2"),  # Leaf C kgC/m2
+      "TotLivBiom" = PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC, "gC/m2", "kgC/m2") + PEcAn.utils::ud_convert(sub.sipnet.output$plantLeafC, "gC/m2", "kgC/m2") + 
+        PEcAn.utils::ud_convert(sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC, "gC/m2", "kgC/m2"), # Total living C kgC/m2
+      "TotSoilCarb" = PEcAn.utils::ud_convert(sub.sipnet.output$soil, "gC/m2", "kgC/m2") + PEcAn.utils::ud_convert(sub.sipnet.output$litter, "gC/m2", "kgC/m2")  # Total soil C kgC/m2
     )
     if (revision == "unk") {
       ## *** NOTE : npp in the sipnet output file is actually evapotranspiration, this is due to a bug in sipnet.c : ***
@@ -164,8 +162,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, "gC/m2", "kgC/m2")  ## litter kgC/m2
+    output[["litter_mass_content_of_water"]] <- (sub.sipnet.output$litterWater * 10) # Litter water kgW/m2 - Note: * 10 is cm->mm model-specific conversion, not unit conversion
     #calculate LAI for standard output
     param <- utils::read.table(file.path(gsub(pattern = "/out/",
                                               replacement = "/run/", x = outdir),
@@ -174,10 +172,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, "gC/m2", "kgC/m2")  ## fine_root_carbon_content kgC/m2
+    output[["coarse_root_carbon_content"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$coarseRootC, "gC/m2", "kgC/m2")  ## coarse_root_carbon_content kgC/m2
+    output[["GWBI"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$woodCreation, "gC/m2/d", "kgC/m2/s") ## kgC/m2/s - this is daily in SIPNET
+    output[["AGB"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC, "gC/m2", "kgC/m2") # Total aboveground biomass kgC/m2
     output[["time_bounds"]] <- c(rbind(bounds[,1], bounds[,2]))
 
     # ******************** Declare netCDF variables ********************#

models/sipnet/R/write_restart.SIPNET.R

@@ -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")  # Unit conversion: kgC/m2/s -> Mg/ha/yr
     names(analysis.save[[length(analysis.save)]]) <- c("NPP")
   }
 

models/stics/R/met2model.STICS.R

@@ -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", "mm/d"), ind, mean, na.rm = TRUE) 
       weather_df[ ,10] <- round(raini, digits = 2) # precipitation (mm d-1)	
 
       # column 11: wind (m.s-1)