diff --git a/base/utils/tests/testthat/test-ud_convert.R b/base/utils/tests/testthat/test-ud_convert.R
index bd994f1fb6a..799b79a6a73 100644
--- a/base/utils/tests/testthat/test-ud_convert.R
+++ b/base/utils/tests/testthat/test-ud_convert.R
@@ -36,4 +36,21 @@ test_that("ud_convert() warns with wrong input units for difftime", {
   expect_warning(ud_convert(as.difftime("12:00:00"), u1 = "years", u2 = "minutes"))
   #should still error if units are not convertible
   expect_error(ud_convert(as.difftime("12:00:00"), u1 = "kilograms", u2 = "minutes"))
+})
+
+test_that("model-specific pool conversions", {
+  # DALEC/SIPNET C pools
+  expect_equal(ud_convert(100, "g/m2", "kg/m2"), 0.1)
+  # GDAY pools
+  expect_equal(ud_convert(10, "Mg/ha", "kg/m2"), 1)
+})
+
+test_that("model-specific flux conversions", {
+  # DALEC/SIPNET C fluxes
+  expect_equal(ud_convert(86400, "g/m2/d", "kg/m2/s"), 0.001, tolerance = 1e-10)
+})
+
+test_that("photosynthesis parameters", {
+  # Photosynthesis energy parameters
+  expect_equal(ud_convert(1000, "J/mol", "kJ/mol"), 1)
 })
\ No newline at end of file
diff --git a/docker/depends/pecan_package_dependencies.csv b/docker/depends/pecan_package_dependencies.csv
index 73a474f78b4..a1b8db4ff8f 100644
--- a/docker/depends/pecan_package_dependencies.csv
+++ b/docker/depends/pecan_package_dependencies.csv
@@ -698,6 +698,7 @@
 "withr","*","base/workflow","Suggests",FALSE
 "withr","*","models/basgra","Suggests",FALSE
 "withr","*","models/ed","Suggests",FALSE
+"withr","*","models/gday","Suggests",FALSE
 "withr","*","models/rothc","Suggests",FALSE
 "withr","*","models/sibcasa","Suggests",FALSE
 "withr","*","models/sipnet","Suggests",FALSE
diff --git a/models/dalec/R/model2netcdf.DALEC.R b/models/dalec/R/model2netcdf.DALEC.R
index 4b40b49fc16..964d2243355 100644
--- a/models/dalec/R/model2netcdf.DALEC.R
+++ b/models/dalec/R/model2netcdf.DALEC.R
@@ -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
diff --git a/models/fates/R/write.configs.FATES.R b/models/fates/R/write.configs.FATES.R
index b58dc649861..871b324fa45 100644
--- a/models/fates/R/write.configs.FATES.R
+++ b/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"))
        }
        
        # 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"))
        }
        
        # 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"))
        }
        
        # 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"))
        }
        
        # deltaS Vcmax - BETY units:J/mol/K;  FATES units: J/mol/K
diff --git a/models/gday/DESCRIPTION b/models/gday/DESCRIPTION
index a7798ed610a..d4490d60efa 100644
--- a/models/gday/DESCRIPTION
+++ b/models/gday/DESCRIPTION
@@ -20,7 +20,8 @@ Imports:
     lubridate (>= 1.6.0),
     ncdf4 (>= 1.15)
 Suggests:
-    testthat (>= 1.0.2)
+    testthat (>= 1.0.2),
+    withr
 SystemRequirements: GDAY
 OS_type: unix
 License: BSD_3_clause + file LICENSE
diff --git a/models/gday/R/model2netcdf.GDAY.R b/models/gday/R/model2netcdf.GDAY.R
index 46ad0fb0670..9c3c6fab37d 100644
--- a/models/gday/R/model2netcdf.GDAY.R
+++ b/models/gday/R/model2netcdf.GDAY.R
@@ -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)
@@ -44,17 +41,20 @@ 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
+    ## NOTE: GDAY outputs daily accumulated values in Mg/ha/day
+    ## Transform all to kg/m2/sec for PEcAn output
+    output[[1]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "auto_resp"], "Mg/ha/day", "kg/m2/s")
+    output[[2]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "hetero_resp"], "Mg/ha/day", "kg/m2/s")
+    output[[3]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "auto_resp"] + sub.GDAY.output[, "hetero_resp"], "Mg/ha/day", "kg/m2/s")
+    output[[4]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "gpp"], "Mg/ha/day", "kg/m2/s")
+    output[[5]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "nep"] * -1, "Mg/ha/day", "kg/m2/s")
+    output[[6]] <- PEcAn.utils::ud_convert(sub.GDAY.output[, "npp"], "Mg/ha/day", "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)
+    # NOTE: GDAY outputs stocks in Mg/ha. 
+    # Transform to kg/m2 for PEcAn output
+    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
diff --git a/models/gday/tests/testthat/data/gday_out.csv b/models/gday/tests/testthat/data/gday_out.csv
new file mode 100644
index 00000000000..a76ea79c61a
--- /dev/null
+++ b/models/gday/tests/testthat/data/gday_out.csv
@@ -0,0 +1,7 @@
+# GDAY output
+year,day,auto_resp,hetero_resp,gpp,nep,npp,stem,branch,soilc,lai,et,transpiration
+2004,1,0.01,0.02,0.05,-0.03,0.03,100,50,200,2,0.5,0.3
+2004,2,0.01,0.02,0.05,-0.03,0.03,100,50,200,2,0.5,0.3
+2004,3,0.01,0.02,0.05,-0.03,0.03,100,50,200,2,0.5,0.3
+2004,4,0.01,0.02,0.05,-0.03,0.03,100,50,200,2,0.5,0.3
+2004,5,0.01,0.02,0.05,-0.03,0.03,100,50,200,2,0.5,0.3
diff --git a/models/gday/tests/testthat/test.model2netcdf.GDAY.R b/models/gday/tests/testthat/test.model2netcdf.GDAY.R
new file mode 100644
index 00000000000..68208b295d4
--- /dev/null
+++ b/models/gday/tests/testthat/test.model2netcdf.GDAY.R
@@ -0,0 +1,54 @@
+##' Test for GDAY model2netcdf unit conversions
+##' 
+##' This test verifies that the unit conversions in model2netcdf.GDAY are correct.
+##' 
+##' Reference: https://github.com/PecanProject/pecan/pull/3719
+##' GDAY outputs daily values in Mg/ha/day, which should be converted to kg/m2/s
+##'
+##' Conversion factors:
+##'   1 Mg/ha = 0.1 kg/m2 (area conversion)
+##'   1 day = 86400 seconds
+##'   Therefore: 1 Mg/ha/day = 0.1 / 86400 kg/m2/s = 1.157e-6 kg/m2/s
+
+context("GDAY model2netcdf unit conversions")
+
+test_that("model2netcdf.GDAY runs without error and produces netCDF", {
+  outdir <- withr::local_tempdir()
+  file.copy("data/gday_out.csv", outdir)
+  
+  # Run the function
+  expect_silent(
+    model2netcdf.GDAY(
+      outdir = outdir,
+      sitelat = 40,
+      sitelon = -88,
+      start_date = "2004-01-01",
+      end_date = "2004-12-31"
+    )
+  )
+  
+  # Check that netCDF file is created
+  nc_file <- file.path(outdir, "2004.nc")
+  expect_true(file.exists(nc_file))
+  
+  # Check that we can read the output
+  output <- PEcAn.utils::read.output(
+    ncfiles = nc_file,
+    variables = c("GPP", "AbvGrndWood"),
+    dataframe = TRUE,
+    verbose = FALSE,
+    print_summary = FALSE
+  )
+  # GPP should be in kg/m2/s (converted from Mg/ha/day)
+  secs_day <- 86400
+  kg_Mg <- 1000
+  m2_ha <- 10000
+  gday2pecan <- kg_Mg/m2_ha/secs_day
+  expect_equal(nrow(output), 5)
+  expect_equal(output$GPP, rep(0.5, 5) * gday2pecan, tolerance = 1e-6)
+  
+  # AbvGrndWood is a stock (Mg/ha), not a flux (Mg/ha/day).
+  # Conversion: 1 Mg/ha = 0.1 kg/m2
+  stock_conv <- 0.1
+  expect_equal(output$AbvGrndWood, rep(150, 5) * stock_conv, tolerance = 1e-6)
+})
diff --git a/models/sipnet/R/model2netcdf.SIPNET.R b/models/sipnet/R/model2netcdf.SIPNET.R
index 5e74de1d726..358adf3951c 100644
--- a/models/sipnet/R/model2netcdf.SIPNET.R
+++ b/models/sipnet/R/model2netcdf.SIPNET.R
@@ -134,21 +134,19 @@ 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 - (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 : ***
@@ -164,8 +162,12 @@ 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")
+    # litterWater was removed in SIPNET v2; only extract if present
+    if ("litterWater" %in% names(sub.sipnet.output)) {
+      # Units are labeled elsewhere as kg water m-2 (which is equivalent to mm, but ud_convert doesn't know that)
+      output[["litter_mass_content_of_water"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$litterWater, "cm", "mm")
+    }
     #calculate LAI for standard output
     # LAI = plantLeafC / leafCSpWt
     # both operands are in carbon units (gC/m2 and gC/m2_leaf),
@@ -176,17 +178,19 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
     leafCSpWt <- param[param[, 1] == "leafCSpWt", 2]
     SLA <- 1000 / leafCSpWt  # m2 leaf / kg C
     output[["LAI"]] <- output[["leaf_carbon_content"]] * SLA
-    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")
+    if ("woodCreation" %in% names(sub.sipnet.output)) {
+      output[["GWBI"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$woodCreation, "g/m2/day", "kg/m2/s")
+    }
+    output[["AGB"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC, "g/m2", "kg/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
+      output[["N2O_flux"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$n2o, "g/m2", "kg/m2") / 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
+      output[["CH4_flux"]] <- PEcAn.utils::ud_convert(sub.sipnet.output$ch4, "g/m2", "kg/m2") / timestep.s
       # convert g C m-2 per timestep -> kg C m-2 s-1
     }
     output[["time_bounds"]] <- c(rbind(bounds[,1], bounds[,2]))
@@ -235,19 +239,22 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
       "SoilMoistFrac" = PEcAn.utils::to_ncvar("SoilMoistFrac", dims),
       "SWE" = PEcAn.utils::to_ncvar("SWE", dims),
       "litter_carbon_content" = PEcAn.utils::to_ncvar("litter_carbon_content", dims),
-      "litter_mass_content_of_water" = PEcAn.utils::to_ncvar("litter_mass_content_of_water", dims),
       "LAI" = PEcAn.utils::to_ncvar("LAI", dims),
       "fine_root_carbon_content" = PEcAn.utils::to_ncvar("fine_root_carbon_content", dims),
       "coarse_root_carbon_content" = PEcAn.utils::to_ncvar("coarse_root_carbon_content", dims),
-      "GWBI" = ncdf4::ncvar_def("GWBI", units = "kg C m-2", dim = list(lon, lat, t), missval = -999,
-                                longname = "Gross Woody Biomass Increment"),
       "AGB" = ncdf4::ncvar_def("AGB", units = "kg C m-2", dim = list(lon, lat, t), missval = -999,
                                longname = "Total aboveground biomass"),
       "time_bounds" = ncdf4::ncvar_def(name="time_bounds", units='',
                                        longname = "history time interval endpoints", dim=list(time_interval,time = t), 
                                        prec = "double")              
     )
-
+    if ("litter_mass_content_of_water" %in% names(output)) {
+      nc_var[["litter_mass_content_of_water"]] <- PEcAn.utils::to_ncvar("litter_mass_content_of_water", dims)
+    }
+    if ("litter_mass_content_of_water" %in% names(output)) {
+      nc_var[["GWBI"]] <- ncdf4::ncvar_def("GWBI", units = "kg C m-2", dim = list(lon, lat, t), missval = -999,
+                                           longname = "Gross Woody Biomass Increment")
+    }
     if ("N2O_flux" %in% names(output)) {
       nc_var[["N2O_flux"]] <- PEcAn.utils::to_ncvar("N2O_flux", dims)
     }
diff --git a/models/sipnet/R/write_restart.SIPNET.R b/models/sipnet/R/write_restart.SIPNET.R
index a3ac00e8a3b..44a6ebf1ade 100755
--- a/models/sipnet/R/write_restart.SIPNET.R
+++ b/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")
     names(analysis.save[[length(analysis.save)]]) <- c("NPP")
   }
 
diff --git a/models/stics/R/met2model.STICS.R b/models/stics/R/met2model.STICS.R
index 1d53bdd13ea..70e2693e297 100644
--- a/models/stics/R/met2model.STICS.R
+++ b/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", "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)