diff --git a/models/ed/R/model2netcdf.ED2.R b/models/ed/R/model2netcdf.ED2.R
index 16a908f40c1..a7ad630dc7b 100644
--- a/models/ed/R/model2netcdf.ED2.R
+++ b/models/ed/R/model2netcdf.ED2.R
@@ -43,6 +43,7 @@ model2netcdf.ED2 <- function(outdir,
                              settings = NULL,
                              process_partial = FALSE) {
   if(!is.null(settings)) {
+    #extract arguments from `settings` if it is supplied
     if(!inherits(settings, "Settings")) {
       PEcAn.logger::logger.error("`settings` should be a PEcAn 'Settings' object")
     }
@@ -126,11 +127,22 @@ model2netcdf.ED2 <- function(outdir,
     # ----- read values from ED output files
     for (i in seq_along(out_list)) {
       rflag <- ed_res_flag[i]
-      # fcnx is either read_T_files() or read_E_files()
-      fcnx  <- paste0("read_", gsub("-", "", rflag), "_files")
-      out_list[[rflag]] <- do.call(fcnx, list(yr = y, ylist[[rflag]], flist[[rflag]],
-                                              outdir, start_date, end_date,
-                                              pfts, settings))
+      fcnx <- switch(
+        rflag,
+        "-T-" = read_T_files,
+        "-E-" = read_E_files
+      )
+      out_list[[rflag]] <-
+        fcnx(
+          yr = y,
+          yfiles = ylist[[rflag]],
+          h5_files = flist[[rflag]],
+          outdir = outdir,
+          start_date = start_date,
+          end_date = end_date,
+          pfts = pfts,
+          settings = settings
+        )
     }
 
     # generate start/end dates for processing
@@ -160,11 +172,21 @@ model2netcdf.ED2 <- function(outdir,
     nc_var <- list()
     for (i in seq_along(out_list)) {
       rflag   <- ed_res_flag[i]
-      #fcnx is either put_T_values() or put_E_values()
-      fcnx    <- paste0("put_", gsub("-", "", rflag), "_values")
-      put_out <- do.call(fcnx, list(yr = y, nc_var = nc_var, var_list = out_list[[rflag]],
-                                    lat = lat, lon = lon, start_date = start_date_real,
-                                    end_date = end_date_real))
+      fcnx <- switch(
+        rflag,
+        "-T-" = put_T_values,
+        "-E-" = put_E_values
+      )
+      put_out <-
+        fcnx(
+          yr = y,
+          nc_var = nc_var,
+          var_list = out_list[[rflag]],
+          lat = lat,
+          lon = lon,
+          start_date = start_date_real,
+          end_date = end_date_real
+        )
       nc_var            <- put_out$nc_var
       out_list[[rflag]] <- put_out$out
     }
@@ -187,7 +209,7 @@ model2netcdf.ED2 <- function(outdir,
     varfile <- file(file.path(outdir, paste(y, "nc", "var", sep = ".")), "w")
     # fill nc file with data
     for (i in seq_along(nc_var)) {
-      var_put(nc, varid = nc_var[[i]], vals = out[[i]])
+      varput(nc, varid = nc_var[[i]], vals = out[[i]]) #TODO: if this errors because of incorrect dimensions, make the whole function error?
       cat(paste(nc_var[[i]]$name, nc_var[[i]]$longname), file = varfile,
           sep = "\n")
     }
@@ -737,6 +759,7 @@ put_T_values <-
     warning("`out` is deprecated, using `var_list` instead")
     var_list <- out
   }
+  #TODO: use append() instead of s + index
   s <- length(nc_var)
   
   # create out list to be modified
@@ -810,118 +833,331 @@ put_T_values <-
 
   # ----- fill list
   
+  #TODO: why is mstmipvar only used for *non-standard* variables??
+  # out <- conversion(1, PEcAn.utils::ud_convert(1, "t ha-1", "kg m-2"))  ## tC/ha -> kg/m2
+  tvars <- list(
+    "AbvGrndWood" = ncdf4::ncvar_def(
+      "AbvGrndWood",
+      units = "kg C m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Above ground woody biomass"
+    ),
+    "AutoResp" = ncdf4::ncvar_def(
+      "AutoResp",
+      units = "kg C m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Autotrophic Respiration"
+    ),
+    "CarbPools" = ncdf4::ncvar_def(
+      "CarbPools",
+      units = "kg C m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Size of each carbon pool"
+    ),
+    "CO2CAS" = ncdf4::ncvar_def(
+      "CO2CAS",
+      units = "ppmv",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "CO2CAS"
+    ),
+    "CropYield" =  ncdf4::ncvar_def(
+      "CropYield",
+      units = "kg m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Crop Yield"
+    ),
+    "GPP" =  ncdf4::ncvar_def(
+      "GPP",
+      units = "kg C m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Gross Primary Productivity"
+    ),
+    "HeteroResp" = ncdf4::ncvar_def(
+      "HeteroResp",
+      units = "kg C m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Heterotrophic Respiration"
+    ),
+  # out <- conversion(8, yr2s)  ## kg C m-2 yr-1 -> kg C m-2 s-1
+    "NEE" =   ncdf4::ncvar_def(
+      "NEE",
+      units = "kg C m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Net Ecosystem Exchange"
+    ),
+  # out <- conversion(9, yr2s)  ## kg C m-2 yr-1 -> kg C m-2 s-1
+    "NPP" = ncdf4::ncvar_def(
+      "NPP",
+      units = "kg C m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Net Primary Productivity"
+    ),
+  # out <- conversion(10, yr2s)  ## kg C m-2 yr-1 -> kg C m-2 s-1
+    "TotalResp" = ncdf4::ncvar_def(
+      "TotalResp",
+      units = "kg C m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Total Respiration"
+    ),
+    "TotLivBiom" = ncdf4::ncvar_def(
+      "TotLivBiom",
+      units = "kg C m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Total living biomass"
+    ),
+    "TotSoilCarb" = ncdf4::ncvar_def(
+      "TotSoilCarb",
+      units = "kg C m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Total Soil Carbon"
+    ),
+    "Fdepth" = ncdf4::ncvar_def(
+      "Fdepth",
+      units = "m",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Frozen Thickness"
+    ),
+    "SnowDepth" = ncdf4::ncvar_def(
+      "SnowDepth",
+      units = "m",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Total snow depth"
+    ),
+    "SnowFrac" =  PEcAn.utils::mstmipvar("SnowFrac", lat, lon, t, zg), # not standard
+    "Tdepth" = ncdf4::ncvar_def(
+      "Tdepth",
+      units = "m",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Active Layer Thickness"
+    ),
+    "CO2air" =  ncdf4::ncvar_def(
+      "CO2air",
+      units = "umol mol-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Near surface CO2 concentration"
+    ),
+    "LWdown" =  ncdf4::ncvar_def(
+      "LWdown",
+      units = "W m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Surface incident longwave radiation"
+    ),
+    "Psurf" =  ncdf4::ncvar_def(
+      "Psurf",
+      units = "Pa",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Surface pressure"
+    ),
+    "Qair" =  ncdf4::ncvar_def(
+      "Qair",
+      units = "kg kg-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Near surface specific humidity"
+    ),
+    "Rainf" = ncdf4::ncvar_def(
+      "Rainf",
+      units = "kg m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Rainfall rate"
+    ),
+    "SWdown" = ncdf4::ncvar_def(
+      "SWdown",
+      units = "W m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Surface incident shortwave radiation"
+    ),
+  # out <- checkTemp(23)
+    "Tair" =  ncdf4::ncvar_def(
+      "Tair",
+      units = "K",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Near surface air temperature"
+    ),
+    "Wind" = ncdf4::ncvar_def(
+      "Wind",
+      units = "m s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Near surface module of the wind"
+    ),
+    "LWnet" =  ncdf4::ncvar_def(
+      "LWnet",
+      units = "W m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Net Longwave Radiation"
+    ),
+    "Qg" =  ncdf4::ncvar_def(
+      "Qg",
+      units = "W m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Ground heat"
+    ),
+    "Qh" =  ncdf4::ncvar_def(
+      "Qh",
+      units = "W m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Sensible heat"
+    ),
+  # out <- conversion(28, PEcAn.data.atmosphere::get.lv())  ## kg m-2 s-1 -> W m-2
+    "Qle" =  ncdf4::ncvar_def(
+      "Qle",
+      units = "W m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Latent heat"
+    ),
+    "SWnet" = ncdf4::ncvar_def(
+      "SWnet",
+      units = "W m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Net shortwave radiation"
+    ),
+    "RootMoist" = PEcAn.utils::mstmipvar("RootMoist", lat, lon, t, zg),   # not standard
+    "TVeg" =  ncdf4::ncvar_def(
+      "TVeg",
+      units = "kg m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Transpiration"
+    ),
+    "WaterTableD" = PEcAn.utils::mstmipvar("WaterTableD", lat, lon, t, zg), # not standard
+    "fPAR" = ncdf4::ncvar_def(
+      "fPAR",
+      units = "",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Absorbed fraction incoming PAR"
+    ),
+    "LAI" =  ncdf4::ncvar_def(
+      "LAI",
+      units = "m2 m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Leaf Area Index"
+    ),
+    "SMFrozFrac" = PEcAn.utils::mstmipvar("SMFrozFrac", lat, lon, t, zg),  # not standard
+    "SMLiqFrac" = PEcAn.utils::mstmipvar("SMLiqFrac", lat, lon, t, zg),   # not standard
+    "SoilMoist" =  ncdf4::ncvar_def(
+      "SoilMoist",
+      units = "kg m-2",
+      dim = list(lon, lat, zg, t),
+      missval = -999,
+      longname = "Average Layer Soil Moisture"
+    ),
+  # out <- checkTemp(38)
+    "SoilTemp" =  ncdf4::ncvar_def(
+      "SoilTemp",
+      units = "K",
+      dim = list(lon, lat, zg, t),
+      missval = -999,
+      longname = "Average Layer Soil Temperature"
+    ),
+    "SoilWet" = ncdf4::ncvar_def(
+      "SoilWet",
+      units = "",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Total Soil Wetness"
+    ),
+    "Albedo" =  PEcAn.utils::mstmipvar("Albedo", lat, lon, t, zg),   # not standard
+  # out <- checkTemp(41)
+    "SnowT" = PEcAn.utils::mstmipvar("SnowT", lat, lon, t, zg),       # not standard
+    "SWE" =  ncdf4::ncvar_def(
+      "SWE",
+      units = "kg m-2",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Snow Water Equivalent"
+    ),
+  # out <- checkTemp(43)
+    "VegT" =  PEcAn.utils::mstmipvar("VegT", lat, lon, t, zg),        # not standard
+    "Evap" = ncdf4::ncvar_def(
+      "Evap",
+      units = "kg m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Total Evaporation"
+    ),
+    "Qs" = ncdf4::ncvar_def(
+      "Qs",
+      units = "kg m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Surface runoff"
+    ),
+    "Qsb" = ncdf4::ncvar_def(
+      "Qsb",
+      units = "kg m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Subsurface runoff"
+    ),
+  # out <- conversion(47, yr2s)  ## kg C m-2 yr-1 -> kg C m-2 s-1
+    "SoilResp" =  ncdf4::ncvar_def(
+      "SoilResp",
+      units = "kg C m-2 s-1",
+      dim = list(lon, lat, t),
+      missval = -999,
+      longname = "Soil Respiration"
+    )
+  )
+  nc_var <- append(nc_var, tvars)
+  # Remove SLZ from output before finalizing list.  replace with time_bounds
+  if(!is.null(out[["SLZ"]])){
+    out[["SLZ"]] <- NULL
+  }
+  #TODO use named lists and append() here
+  out_length <- length(out)
+  out[[out_length + 1]] <- c(rbind(bounds[, 1], bounds[, 2]))
+  nc_var[[s + (out_length + 1)]] <- ncdf4::ncvar_def(name="time_bounds", units='', 
+                                    longname = "history time interval endpoints", 
+                                    dim=list(time_interval,time = t), 
+                                    prec = "double")
+  
+  
+  #TODO: change conversions to usse named vector instead of positions? Should happen in read_T_files, I think?
   out <- conversion(1, PEcAn.utils::ud_convert(1, "t ha-1", "kg m-2"))  ## tC/ha -> kg/m2
-  nc_var[[s + 1]] <- ncdf4::ncvar_def("AbvGrndWood", units = "kg C m-2", dim = list(lon, lat, t), missval = -999, 
-                                    longname = "Above ground woody biomass")
   out <- conversion(2, umol2kg_C)  ## umol/m2 s-1 -> kg/m2 s-1
-  nc_var[[s + 2]] <- ncdf4::ncvar_def("AutoResp", units = "kg C m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                    longname = "Autotrophic Respiration")
-  nc_var[[s + 3]] <- ncdf4::ncvar_def("CarbPools", units = "kg C m-2", dim = list(lon, lat, t), missval = -999, 
-                                    longname = "Size of each carbon pool")
-  nc_var[[s + 4]] <- ncdf4::ncvar_def("CO2CAS", units = "ppmv", dim = list(lon, lat, t), missval = -999, 
-                                    longname = "CO2CAS")
-  nc_var[[s + 5]] <- ncdf4::ncvar_def("CropYield", units = "kg m-2", dim = list(lon, lat, t), missval = -999, 
-                                    longname = "Crop Yield")
   out <- conversion(6, yr2s)  ## kg C m-2 yr-1 -> kg C m-2 s-1
-  nc_var[[s + 6]] <- ncdf4::ncvar_def("GPP", units = "kg C m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                    longname = "Gross Primary Productivity")
   out <- conversion(7, yr2s)  ## kg C m-2 yr-1 -> kg C m-2 s-1
-  nc_var[[s + 7]] <- ncdf4::ncvar_def("HeteroResp", units = "kg C m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                    longname = "Heterotrophic Respiration")
   out <- conversion(8, yr2s)  ## kg C m-2 yr-1 -> kg C m-2 s-1
-  nc_var[[s + 8]] <-  ncdf4::ncvar_def("NEE", units = "kg C m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Net Ecosystem Exchange")
   out <- conversion(9, yr2s)  ## kg C m-2 yr-1 -> kg C m-2 s-1
-  nc_var[[s + 9]] <- ncdf4::ncvar_def("NPP", units = "kg C m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                    longname = "Net Primary Productivity")
   out <- conversion(10, yr2s)  ## kg C m-2 yr-1 -> kg C m-2 s-1
-  nc_var[[s + 10]] <- ncdf4::ncvar_def("TotalResp", units = "kg C m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Total Respiration")
-  nc_var[[s + 11]] <- ncdf4::ncvar_def("TotLivBiom", units = "kg C m-2", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Total living biomass")
-  nc_var[[s + 12]] <- ncdf4::ncvar_def("TotSoilCarb", units = "kg C m-2", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Total Soil Carbon")
-  nc_var[[s + 13]] <- ncdf4::ncvar_def("Fdepth", units = "m", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Frozen Thickness")
-  nc_var[[s + 14]] <- ncdf4::ncvar_def("SnowDepth", units = "m", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Total snow depth")
-  nc_var[[s + 15]] <- PEcAn.utils::mstmipvar("SnowFrac", lat, lon, t, zg) # not standard
-  nc_var[[s + 16]] <- ncdf4::ncvar_def("Tdepth", units = "m", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Active Layer Thickness")
-  nc_var[[s + 17]] <- ncdf4::ncvar_def("CO2air", units = "umol mol-1", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Near surface CO2 concentration")
-  nc_var[[s + 18]] <- ncdf4::ncvar_def("LWdown", units = "W m-2", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Surface incident longwave radiation")
-  nc_var[[s + 19]] <- ncdf4::ncvar_def("Psurf", units = "Pa", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Surface pressure")
-  nc_var[[s + 20]] <- ncdf4::ncvar_def("Qair", units = "kg kg-1", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Near surface specific humidity")
-  nc_var[[s + 21]] <- ncdf4::ncvar_def("Rainf", units = "kg m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Rainfall rate")
-  nc_var[[s + 22]] <- ncdf4::ncvar_def("SWdown", units = "W m-2", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Surface incident shortwave radiation")
   out <- checkTemp(23)
-  nc_var[[s + 23]] <- ncdf4::ncvar_def("Tair", units = "K", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Near surface air temperature")
-  nc_var[[s + 24]] <- ncdf4::ncvar_def("Wind", units = "m s-1", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Near surface module of the wind")
-  nc_var[[s + 25]] <- ncdf4::ncvar_def("LWnet", units = "W m-2", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Net Longwave Radiation")
-  nc_var[[s + 26]] <- ncdf4::ncvar_def("Qg", units = "W m-2", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Ground heat")
-  nc_var[[s + 27]] <- ncdf4::ncvar_def("Qh", units = "W m-2", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Sensible heat")
   out <- conversion(28, PEcAn.data.atmosphere::get.lv())  ## kg m-2 s-1 -> W m-2
-  nc_var[[s + 28]] <- ncdf4::ncvar_def("Qle", units = "W m-2", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Latent heat")
-  nc_var[[s + 29]] <- ncdf4::ncvar_def("SWnet", units = "W m-2", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Net shortwave radiation")
-  nc_var[[s + 30]] <- PEcAn.utils::mstmipvar("RootMoist", lat, lon, t, zg)   # not standard
-  nc_var[[s + 31]] <- ncdf4::ncvar_def("TVeg", units = "kg m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Transpiration")
-  nc_var[[s + 32]] <- PEcAn.utils::mstmipvar("WaterTableD", lat, lon, t, zg) # not standard
-
-  nc_var[[s + 33]] <- ncdf4::ncvar_def("fPAR", units = "", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Absorbed fraction incoming PAR")
-  nc_var[[s + 34]] <- ncdf4::ncvar_def("LAI", units = "m2 m-2", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Leaf Area Index")
-  nc_var[[s + 35]] <- PEcAn.utils::mstmipvar("SMFrozFrac", lat, lon, t, zg)  # not standard
-  nc_var[[s + 36]] <- PEcAn.utils::mstmipvar("SMLiqFrac", lat, lon, t, zg)   # not standard
-  nc_var[[s + 37]] <- ncdf4::ncvar_def("SoilMoist", units = "kg m-2", dim = list(lon, lat, zg, t), missval = -999, 
-                                     longname = "Average Layer Soil Moisture")
   out <- checkTemp(38)
-  nc_var[[s + 38]] <- ncdf4::ncvar_def("SoilTemp", units = "K", dim = list(lon, lat, zg, t), missval = -999, 
-                                     longname = "Average Layer Soil Temperature")
-  nc_var[[s + 39]] <- ncdf4::ncvar_def("SoilWet", units = "", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Total Soil Wetness")
-  nc_var[[s + 40]] <- PEcAn.utils::mstmipvar("Albedo", lat, lon, t, zg)      # not standard
   out <- checkTemp(41)
-  nc_var[[s + 41]] <- PEcAn.utils::mstmipvar("SnowT", lat, lon, t, zg)       # not standard
-  nc_var[[s + 42]] <- ncdf4::ncvar_def("SWE", units = "kg m-2", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Snow Water Equivalent")
   out <- checkTemp(43)
-  nc_var[[s + 43]] <- PEcAn.utils::mstmipvar("VegT", lat, lon, t, zg)        # not standard
-  nc_var[[s + 44]] <- ncdf4::ncvar_def("Evap", units = "kg m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Total Evaporation")
-  nc_var[[s + 45]] <- ncdf4::ncvar_def("Qs", units = "kg m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Surface runoff")
-  nc_var[[s + 46]] <- ncdf4::ncvar_def("Qsb", units = "kg m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                     longname = "Subsurface runoff")
   out <- conversion(47, yr2s)  ## kg C m-2 yr-1 -> kg C m-2 s-1
-  nc_var[[s + 47]] <- ncdf4::ncvar_def("SoilResp", units = "kg C m-2 s-1", dim = list(lon, lat, t), missval = -999, 
-                                    longname = "Soil Respiration")
-  # Remove SLZ from output before finalizing list.  replace with time_bounds
-  if(!is.null(out[["SLZ"]])){
-    out[["SLZ"]] <- NULL
-  }
-  out_length <- length(out)
-  out[[out_length + 1]] <- c(rbind(bounds[, 1], bounds[, 2]))
-  nc_var[[s + (out_length + 1)]] <- ncdf4::ncvar_def(name="time_bounds", units='', 
-                                    longname = "history time interval endpoints", 
-                                    dim=list(time_interval,time = t), 
-                                    prec = "double")
   
   return(list(nc_var = nc_var, out = out))
+  #TODO: check that order of nc_var and out match
 
 } # put_T_values
 
@@ -1025,7 +1261,7 @@ read_E_files <- function(yr, yfiles, h5_files, outdir, start_date, end_date,
   
   # loop over the files for that year
   for(i in ysel) {
-    
+
     nc <- ncdf4::nc_open(file.path(outdir, h5_files[i]))
     on.exit(ncdf4::nc_close(nc), add = FALSE)
     allvars <- names(nc$var)
@@ -1065,6 +1301,7 @@ read_E_files <- function(yr, yfiles, h5_files, outdir, start_date, end_date,
     
   } # end ysel-loop
 
+  #TODO: warn if a pft in `pfts` isn't in the h5 file and remove the unused PFT from `pfts`
   
   # even if this is a SA run for soil, currently we are not reading any variable
   # that has a soil dimension. "soil" will be passed to read.output as pft.name
@@ -1164,7 +1401,7 @@ read_E_files <- function(yr, yfiles, h5_files, outdir, start_date, end_date,
 ##' 
 put_E_values <-
   function(yr,
-           nc_var,
+           nc_var, #TODO shouldn't append nc_var or var_list.  Appending should be done in main model2netcdf.ED2 function.
            var_list,
            lat,
            lon,
@@ -1266,50 +1503,51 @@ put_E_values <-
   # output by read_E_files.  Some day the whole model2netcdf.ED2 function should
   # probably be re-written to combine the read_*_files and put_*_values
   # functions to make this harder to accidentally screw up.
+  # TODO: find less fragile way to keep `var_list` and `evars` in the same order
  evars <- list(
-   ncdf4::ncvar_def(
+   "AGB_PFT"= ncdf4::ncvar_def(
      "AGB_PFT", #original ED2 name: AGB_CO
      units = "kgC m-2",
      dim = list(lon, lat, t, p),
      missval = -999,
      longname = "Above ground biomass by PFT"
    ),
-   ncdf4::ncvar_def(
+   "BSEEDS"= ncdf4::ncvar_def(
      "BSEEDS", #original ED2 name: BSEEDS_CO
      units = "kgC m-2",
      dim = list(lon, lat, t, p),
      missval = -999,
      longname = "Seed biomass by PFT"
    ),
-    ncdf4::ncvar_def(
+    "DBH" = ncdf4::ncvar_def(
       "DBH", #original ED2 name: DBH
       units = "cm",
       dim = list(lon, lat, t, p),
       missval = -999,
       longname = "Diameter at breast height by PFT"
     ),
-    ncdf4::ncvar_def(
+    "DDBH" = ncdf4::ncvar_def(
       "DDBH", #original ED2 name: DDBH_DT
       units = "cm yr-1",
       dim = list(lon, lat, t, p),
       missval = -999,
       longname = "Rate of change in dbh by PFT"
     ),
-   ncdf4::ncvar_def(
+   "NPP_PFT" = ncdf4::ncvar_def(
      "NPP_PFT", #original ED2 name: MMEAN_NPPDAILY_CO
      units = "KgC m-2 s-1",
      dim = list(lon, lat, t, p),
      missval = -999,
      longname = "Net primary productivity by PFT"
    ),
-    ncdf4::ncvar_def(
+   "TRANSP_PFT" = ncdf4::ncvar_def(
       "TRANSP_PFT", #original ED2 name: MMEAN_TRANSP_CO
       units = "kg m-2 s-1",
       dim = list(lon, lat, t, p),
       missval = -999,
       longname = "Leaf transpiration by PFT"
     ),
-   ncdf4::ncvar_def(
+   "DENS" = ncdf4::ncvar_def(
      "DENS", #original ED2 name: NPLANT
      units = "plant m-2",
      dim = list(lon, lat, t, p),
@@ -1319,13 +1557,13 @@ put_E_values <-
 
     # longname of this variable will be parsed by read.output
     # so that read.output has a way of accessing PFT names
-    ncdf4::ncvar_def(
+   "PFT" = ncdf4::ncvar_def(
       "PFT",
       units = "",
       dim = list(p),
       longname = paste(names(pfts), collapse = ",")
     ),
-    ncdf4::ncvar_def(
+   "dtime_bounds" = ncdf4::ncvar_def(
       name = "dtime_bounds",
       units = "",
       longname = "monthly history time interval endpoints",
@@ -1334,6 +1572,7 @@ put_E_values <-
     )
  )
  #TODO: assure that nc_var and var_list are of same length and same order?
+ #TODO: don't append to nc_var and var_list here.  Move that to main model2netcdf.ED2 function
  nc_var <- append(nc_var, evars)
  var_list <- append(var_list, list(dtime_bounds = c(bounds)))
  
@@ -1542,14 +1781,13 @@ extract_pfts <- function(pfts) {
   pfts_out
 }
 
-
 # A version of ncvar_put that returns the varid in warning messages
-var_put <- function(nc, varid, vals, start = NA, count = NA) {
-  output <- utils::capture.output(
-    ncdf4::ncvar_put(nc = nc, varid = varid, vals = vals, start = start, count = count)
+varput <- function(nc, varid, vals, start = NA, count = NA) {
+  output <- capture.output(
+    ncvar_put(nc = nc, varid = varid, vals = vals, start = start, count = count)
   )
   if(length(output)!=0) {
-    cat(paste0("With '", varid$name, "':"), output, "\n")
+    cat(paste0("Message for var '", varid$name, "':"), output)
   }
 }
 
diff --git a/models/ed/tests/testthat/test.model2netcdf.ED2.R b/models/ed/tests/testthat/test.model2netcdf.ED2.R
index 03ccddc7e01..7a7241166bf 100644
--- a/models/ed/tests/testthat/test.model2netcdf.ED2.R
+++ b/models/ed/tests/testthat/test.model2netcdf.ED2.R
@@ -4,23 +4,36 @@ testdir <- tempfile()
 dir.create(testdir)
 withr::defer(unlink(testdir, recursive = TRUE))
 unzip("data/outdir.zip", exdir = testdir)
+#for interactive debugging:
+# unzip("models/ed/tests/testthat/data/outdir.zip", exdir = testdir)
 settings <-
   PEcAn.settings::read.settings(file.path(testdir, "outdir", "pecan_checked.xml"))
 settings$outdir <- file.path(testdir, "outdir")
 #not to be confused with outdir
 outdir <- file.path(settings$outdir, "out", "ENS-00001-76")
 
+# set logger level so test output is more readable
+old_level <- PEcAn.logger::logger.setLevel("ERROR")
+withr::defer(PEcAn.logger::logger.setLevel(old_level))
+
 test_that("model2netcdf.ED2 runs without error", {
-  
   #hacky way to check for errors b/c PEcAn.logger errors are non-standard and
   #not captured by testthat::expect_message() or expect_error()
-  x <- capture.output(
+  msg <- capture.output(
     model2netcdf.ED2(outdir = outdir,
                      settings = settings),
     type = "message"
   )
+  expect_false(any(str_detect(msg, "ERROR")))
+})
 
-  expect_false(any(str_detect(x, "ERROR")))
+test_that("no ncvar_put warnings", {
+  out <- capture.output(
+    model2netcdf.ED2(outdir = outdir,
+                     settings = settings),
+    type = "output"
+  )
+  expect_false(any(str_detect(out, "ncvar_put: warning:")))
 })
 
 #remove .nc file
@@ -111,4 +124,13 @@ test_that("variables have MsTMIP standard units",{
   )
 })
 
-
+test_that("Errors when PFTs in settings aren't in .h5 files", {
+  #weird edge case: https://github.com/PecanProject/pecan/issues/3068
+  #fake this by manually editing settings
+  settings$pfts <-
+    append(settings$pfts, list(pft = list(
+      name = "temperate.Early_Hardwood", ed2_pft_number = "9"
+    )))
+  expect_error(model2netcdf.ED2(outdir = outdir, settings = settings))
+  #TODO: also expect no output .nc files.  Currently this error doesn't stop the function from writing the .nc files for some reason
+})