diff --git a/Biomass_core.R b/Biomass_core.R
index 8f5b259..a65b0fd 100644
--- a/Biomass_core.R
+++ b/Biomass_core.R
@@ -146,7 +146,7 @@ defineModule(sim, list(
                                     "trunk/tests/v6.0-2.0/ecoregions.txt")),
     expectsInput("ecoregionMap", "RasterLayer",
                  desc = paste("ecoregion map that has mapcodes match ecoregion table and speciesEcoregion table.",
-                              "Defaults to a dummy map matching rasterToMatch with two regions")),
+                              "Defaults to a dummy map matching rasterToMatch with two regions. Must be a factor raster.")),
     # expectsInput("initialCommunities", "data.table",
     #              desc = "initial community table",
     #              sourceURL = "https://raw.githubusercontent.com/LANDIS-II-Foundation/Extensions-Succession/master/biomass-succession-archive/trunk/tests/v6.0-2.0/initial-communities.txt"),
@@ -524,9 +524,18 @@ Init <- function(sim, verbose = getOption("LandR.verbose", TRUE)) {
   ##############################################
   ## species
   ##############################################
+  if (!is.data.table(sim$species)) setDT(sim$species)
   if (is.null(sim$species))
     stop("'species' object must be provided")
-  species <- setDT(sim$species)[, speciesCode := as.factor(species)]
+  if ("speciesCode" %in% colnames(sim$species)) {
+    if ("species" %in% colnames(sim$species))
+      sim$species[, species := as.character(speciesCode)]
+  } else {
+    if (!"species" %in% colnames(sim$species))
+      stop("species object must have either species or speciesCode column, or both")
+    sim$species[, speciesCode := as.factor(species)]
+  }
+  species <- sim$species
   LandR::assertColumns(species,
                        c(species = "character", Area = "factor", longevity = "integer",
                          sexualmature = "integer", shadetolerance = "numeric",
@@ -652,6 +661,11 @@ Init <- function(sim, verbose = getOption("LandR.verbose", TRUE)) {
     biomassModel <- quote(lme4::lmer(B ~ logAge * speciesCode + cover * speciesCode +
                                        (logAge + cover + speciesCode | ecoregionGroup)))
 
+    coverModel <- quote(lme4::glmer(cbind(coverPres, coverNum) ~
+                                      (1 | ecoregionGroup), family = binomial))
+    biomassModel <- quote(lme4::lmer(B ~ logAge  + cover +
+                                       (logAge + cover  | ecoregionGroup)))
+
     ## COVER
     message(blue("Estimating Species Establishment Probability from "), red("DUMMY values of ecoregionGroup "),
             blue("using the formula:\n"), magenta(format(coverModel)))
@@ -742,10 +756,10 @@ Init <- function(sim, verbose = getOption("LandR.verbose", TRUE)) {
     LandR::assertUniqueCohortData(sim$cohortData, c("pixelGroup", "ecoregionGroup", "speciesCode"))
   }
 
-  if (!is.null(sim$ecoregionMap) && !is.null(sim$pixelGroupMap) && !is.null(sim$biomassMap)) {
-    compareRaster(sim$biomassMap, sim$ecoregionMap, sim$pixelGroupMap, sim$rasterToMatch, orig = TRUE)
+  if (!is.null(sim$ecoregionMap) && !is.null(sim$pixelGroupMap)) {
+    compareRaster(sim$ecoregionMap, sim$pixelGroupMap, sim$rasterToMatch, orig = TRUE)
   } else {
-    stop("Expecting 3 rasters at this point: sim$biomassMap, sim$ecoregionMap, ",
+    stop("Expecting 2 rasters at this point: sim$ecoregionMap, ",
          "sim$pixelGroupMap and they must match sim$rasterToMatch")
   }
 
@@ -885,6 +899,7 @@ Init <- function(sim, verbose = getOption("LandR.verbose", TRUE)) {
       set(cohortData, NULL, "B", asInteger(cohortData[["B"]]))
   }
 
+
   pixelAll <- cohortData[, .(uniqueSumB = sum(B, na.rm = TRUE)), by = pixelGroup]
   if (!is.integer(pixelAll[["uniqueSumB"]]))
     set(pixelAll, NULL, "uniqueSumB", asInteger(pixelAll[["uniqueSumB"]]))
@@ -935,7 +950,8 @@ Init <- function(sim, verbose = getOption("LandR.verbose", TRUE)) {
   return(invisible(sim))
 }
 
-SummaryBGM <- compiler::cmpfun(function(sim) {
+SummaryBGM <- #compiler::cmpfun(
+  function(sim) {
   pixelGroups <- data.table(pixelGroupIndex = unique(sim$cohortData$pixelGroup),
                             temID = 1:length(unique(sim$cohortData$pixelGroup)))
   cutpoints <- sort(unique(c(seq(1, max(pixelGroups$temID), by = P(sim)$cutpoint),
@@ -1019,7 +1035,8 @@ SummaryBGM <- compiler::cmpfun(function(sim) {
 
   rm(cutpoints, pixelGroups, tempOutput_All, summaryBGMtable) ## TODO: is this needed? on exit, should free the mem used for these
   return(invisible(sim))
-})
+}
+# )
 
 MortalityAndGrowth <- compiler::cmpfun(function(sim) {
   # If cohortData has length 0, don't do this -- this can happen in more theoretical use cases where e.g., end(sim) is longer than longevity
@@ -1223,7 +1240,8 @@ MortalityAndGrowth <- compiler::cmpfun(function(sim) {
     LandR::assertCohortData(sim$cohortData, sim$pixelGroupMap, cohortDefinitionCols = P(sim)$cohortDefinitionCols)
   }
   return(invisible(sim))
-})
+}
+#)
 
 Dispersal <- function(sim) {
 
@@ -1246,7 +1264,8 @@ Dispersal <- function(sim) {
   return(invisible(sim))
 }
 
-NoDispersalSeeding <- compiler::cmpfun(function(sim, tempActivePixel, pixelsFromCurYrBurn) {
+NoDispersalSeeding <- #compiler::cmpfun(
+  function(sim, tempActivePixel, pixelsFromCurYrBurn) {
   # if (sim$lastFireYear == round(time(sim))) { # if current year is both fire year and succession year
   #   # find new active pixel that remove successful postfire regeneration
   #   # since this is on site regeneration, all the burnt pixels can not seeding
@@ -1312,27 +1331,50 @@ NoDispersalSeeding <- compiler::cmpfun(function(sim, tempActivePixel, pixelsFrom
 
   sim$lastReg <- round(time(sim))
   return(invisible(sim))
-})
+}#)
 
-UniversalDispersalSeeding <- compiler::cmpfun(function(sim, tempActivePixel) {
+UniversalDispersalSeeding <- #compiler::cmpfun(
+  function(sim, tempActivePixel, pixelsFromCurYrBurn) {
   # if (sim$lastFireYear == round(time(sim))) { # the current year is both fire year and succession year
   #   tempActivePixel <- sim$activePixelIndex[!(sim$activePixelIndex %in% sim$postFirePixel)]
   # } else {
   #   tempActivePixel <- sim$activePixelIndex
   # }
-  sim$cohortData <- calculateSumB(sim$cohortData, lastReg = sim$lastReg, currentTime = round(time(sim)),
-                                  successionTimestep = P(sim)$successionTimestep)
-  species <- sim$species
+
+    sim$cohortData <- calculateSumB(cohortData = sim$cohortData,
+                                    lastReg = sim$lastReg, currentTime = round(time(sim)),
+                                    successionTimestep = P(sim)$successionTimestep)
+    if (is.null(sim$speciesEcoregion)) stop("WardDispersalSeeding requires a speciesEcoregion")
+    siteShade <- calcSiteShade(currentTime = round(time(sim)), cohortData = sim$cohortData,
+                               sim$speciesEcoregion, sim$minRelativeB)
+    activePixelGroup <- data.table(pixelGroup = unique(getValues(sim$pixelGroupMap)[tempActivePixel])) %>%
+      na.omit()
+    siteShade <- siteShade[activePixelGroup, on = "pixelGroup"]
+    siteShade[is.na(siteShade), siteShade := 0]
+
+    ## Seed source cells:
+    ## 1. Select only sexually mature cohorts, then
+    ## 2. collapse to pixelGroup by species, i.e,. doesn't matter that there is >1 cohort of same species
+    sim$cohortData <- sim$species[, c("speciesCode", "sexualmature")][sim$cohortData,
+                                                                      on = "speciesCode"]
+    # sim$cohortData <- setkey(sim$cohortData, speciesCode)[setkey(sim$species[, .(speciesCode, sexualmature)],
+    #                                                              speciesCode),
+    #                                                       nomatch = 0]
+
+
+  # sim$cohortData <- calculateSumB(sim$cohortData, lastReg = sim$lastReg, currentTime = round(time(sim)),
+  #                                 successionTimestep = P(sim)$successionTimestep)
+  # species <- sim$species
   # all species can provide seed source, i.e. age>=sexualmature
   speciessource <- setkey(sim$species[, .(speciesCode, k = 1)], k)
-  siteShade <- data.table(calcSiteShade(currentTime = round(time(sim)), sim$cohortData,
-                                        sim$speciesEcoregion, sim$minRelativeB))
-  activePixelGroup <- unique(data.table(pixelGroup = getValues(sim$pixelGroupMap)[tempActivePixel],
-                                        ecoregionGroup = factorValues2(sim$ecoregionMap, getValues(sim$ecoregionMap),
-                                                                       att = "ecoregionGroup")[tempActivePixel]),
-                             by = "pixelGroup")
-  siteShade <- dplyr::left_join(activePixelGroup, siteShade, by = "pixelGroup") %>% data.table()
-  siteShade[is.na(siteShade), siteShade := 0]
+  # siteShade <- data.table(calcSiteShade(currentTime = round(time(sim)), sim$cohortData,
+  #                                       sim$speciesEcoregion, sim$minRelativeB))
+  # activePixelGroup <- unique(data.table(pixelGroup = getValues(sim$pixelGroupMap)[tempActivePixel],
+  #                                       ecoregionGroup = factorValues2(sim$ecoregionMap, getValues(sim$ecoregionMap),
+  #                                                                      att = "ecoregionGroup")[tempActivePixel]),
+  #                            by = "pixelGroup")
+  # siteShade <- dplyr::left_join(activePixelGroup, siteShade, by = "pixelGroup") %>% data.table()
+  # siteShade[is.na(siteShade), siteShade := 0]
   setkey(siteShade[, k := 1], k)
   # i believe this is the latest version how the landis guys calculate sufficient light
   # http://landis-extensions.googlecode.com/svn/trunk/succession-library/trunk/src/ReproductionDefaults.cs
@@ -1383,148 +1425,150 @@ UniversalDispersalSeeding <- compiler::cmpfun(function(sim, tempActivePixel) {
   }
   sim$lastReg <- round(time(sim))
   return(invisible(sim))
-})
-
-WardDispersalSeeding <- compiler::cmpfun(function(sim, tempActivePixel, pixelsFromCurYrBurn,
-                                                  verbose = getOption("LandR.verbose", TRUE)) {
-
-  sim$cohortData <- calculateSumB(cohortData = sim$cohortData,
-                                  lastReg = sim$lastReg, currentTime = round(time(sim)),
-                                  successionTimestep = P(sim)$successionTimestep)
-  siteShade <- calcSiteShade(currentTime = round(time(sim)), cohortData = sim$cohortData,
-                             sim$speciesEcoregion, sim$minRelativeB)
-  activePixelGroup <- data.table(pixelGroup = unique(getValues(sim$pixelGroupMap)[tempActivePixel])) %>%
-    na.omit()
-  siteShade <- siteShade[activePixelGroup, on = "pixelGroup"]
-  siteShade[is.na(siteShade), siteShade := 0]
-
-  ## Seed source cells:
-  ## 1. Select only sexually mature cohorts, then
-  ## 2. collapse to pixelGroup by species, i.e,. doesn't matter that there is >1 cohort of same species
-  sim$cohortData <- sim$species[, c("speciesCode", "sexualmature")][sim$cohortData,
-                                                                    on = "speciesCode"]
-  # sim$cohortData <- setkey(sim$cohortData, speciesCode)[setkey(sim$species[, .(speciesCode, sexualmature)],
-  #                                                              speciesCode),
-  #                                                       nomatch = 0]
-  matureCohorts <- sim$cohortData[age >= sexualmature] %>%
-    unique(by = c("pixelGroup", "speciesCode")) %>%
-    setkey(., speciesCode)
-  matureCohorts <- matureCohorts[, .(pixelGroup, speciesCode)]
-  set(sim$cohortData, NULL, "sexualmature", NULL)
-
-  if (NROW(matureCohorts) > 0) {
-    seedSource <- sim$species[, list(speciesCode, seeddistance_eff, seeddistance_max)] %>%
-      setkey(., speciesCode) %>%
-      .[matureCohorts]
-    setkey(seedSource, speciesCode)
-    #  Seed Receiving cells:
-    #  1. Must be sufficient light
-    # seed receive just for the species that are seed source
-    tempspecies1 <- sim$species[speciesCode %in% unique(matureCohorts$speciesCode),][
-      , .(speciesCode, shadetolerance, seeddistance_eff, seeddistance_max)]
-    seedReceive <- setkey(tempspecies1[, c(k = 1, .SD)], k)[setkey(siteShade[
-      , c(k = 1, .SD)], k), allow.cartesian = TRUE][, k := NULL]
-    seedReceive <- assignLightProb(sufficientLight = sim$sufficientLight, seedReceive)
-    set(seedReceive, NULL, "siteShade", NULL)
-
-    # 3. Remove any species from the seedSource that couldn't regeneration anywhere on the map due to insufficient light
-    seedReceive <- seedReceive[lightProb %>>% runif(NROW(seedReceive), 0, 1), ][
-      , .(pixelGroup, speciesCode, seeddistance_eff, seeddistance_max)]
-    setkey(seedReceive, speciesCode)
-
-    # rm ones that had successful serotiny or resprouting
-    seedReceive <- seedReceive[!sim$cohortData[age == 1L], on = c("pixelGroup", "speciesCode")]
-    #    (info contained within seedReceive)
-    # this is should be a inner join, needs to specify the nomatch=0, nomatch = NA is default that sugest the full joint.
-    seedSource <- seedSource[speciesCode %in% unique(seedReceive$speciesCode),]
-
-    # it could be more efficient if sim$pixelGroupMap is reduced map by removing the pixels that have
-    # successful postdisturbance regeneration and the inactive pixels
-    # how to subset the reducedmap
-    # if (sim$lastFireYear == round(time(sim))) { # the current year is both fire year and succession year
-    #   inactivePixelIndex <- c(sim$inactivePixelIndex, sim$treedFirePixelTableSinceLastDisp$pixelIndex)
-    # } else {
-    #   inactivePixelIndex <- sim$inactivePixelIndex
-    # }
-    reducedPixelGroupMap <- sim$pixelGroupMap
-
-    # Calculate the maximum size of the chunks for LANDISDisp
-    if (length(pixelsFromCurYrBurn) > 0) {
-      reducedPixelGroupMap[pixelsFromCurYrBurn] <- NA
-    }
-
-    seedingData <- LANDISDisp(dtRcv = seedReceive, plot.it = FALSE,
-                              dtSrc = seedSource,
-                              speciesTable = sim$species,
-                              pixelGroupMap = reducedPixelGroupMap,
-                              successionTimestep = P(sim)$successionTimestep,
-                              verbose = getOption("LandR.verbose", TRUE) > 0)
+}#)
+
+WardDispersalSeeding <- # compiler::cmpfun(
+  function(sim, tempActivePixel, pixelsFromCurYrBurn,
+           verbose = getOption("LandR.verbose", TRUE)) {
+
+    sim$cohortData <- calculateSumB(cohortData = sim$cohortData,
+                                    lastReg = sim$lastReg, currentTime = round(time(sim)),
+                                    successionTimestep = P(sim)$successionTimestep)
+    if (is.null(sim$speciesEcoregion)) stop("WardDispersalSeeding requires a speciesEcoregion")
+    siteShade <- calcSiteShade(currentTime = round(time(sim)), cohortData = sim$cohortData,
+                               sim$speciesEcoregion, sim$minRelativeB)
+    activePixelGroup <- data.table(pixelGroup = unique(getValues(sim$pixelGroupMap)[tempActivePixel])) %>%
+      na.omit()
+    siteShade <- siteShade[activePixelGroup, on = "pixelGroup"]
+    siteShade[is.na(siteShade), siteShade := 0]
+
+    ## Seed source cells:
+    ## 1. Select only sexually mature cohorts, then
+    ## 2. collapse to pixelGroup by species, i.e,. doesn't matter that there is >1 cohort of same species
+    sim$cohortData <- sim$species[, c("speciesCode", "sexualmature")][sim$cohortData,
+                                                                      on = "speciesCode"]
+    # sim$cohortData <- setkey(sim$cohortData, speciesCode)[setkey(sim$species[, .(speciesCode, sexualmature)],
+    #                                                              speciesCode),
+    #                                                       nomatch = 0]
+    matureCohorts <- sim$cohortData[age >= sexualmature] %>%
+      unique(by = c("pixelGroup", "speciesCode")) %>%
+      setkey(., speciesCode)
+    matureCohorts <- matureCohorts[, .(pixelGroup, speciesCode)]
+    set(sim$cohortData, NULL, "sexualmature", NULL)
+
+    if (NROW(matureCohorts) > 0) {
+      seedSource <- sim$species[, list(speciesCode, seeddistance_eff, seeddistance_max)] %>%
+        setkey(., speciesCode) %>%
+        .[matureCohorts]
+      setkey(seedSource, speciesCode)
+      #  Seed Receiving cells:
+      #  1. Must be sufficient light
+      # seed receive just for the species that are seed source
+      tempspecies1 <- sim$species[speciesCode %in% unique(matureCohorts$speciesCode),][
+        , .(speciesCode, shadetolerance, seeddistance_eff, seeddistance_max)]
+      seedReceive <- setkey(tempspecies1[, c(k = 1, .SD)], k)[setkey(siteShade[
+        , c(k = 1, .SD)], k), allow.cartesian = TRUE][, k := NULL]
+      seedReceive <- assignLightProb(sufficientLight = sim$sufficientLight, seedReceive)
+      set(seedReceive, NULL, "siteShade", NULL)
+
+      # 3. Remove any species from the seedSource that couldn't regeneration anywhere on the map due to insufficient light
+      seedReceive <- seedReceive[lightProb %>>% runif(NROW(seedReceive), 0, 1), ][
+        , .(pixelGroup, speciesCode, seeddistance_eff, seeddistance_max)]
+      setkey(seedReceive, speciesCode)
+
+      # rm ones that had successful serotiny or resprouting
+      seedReceive <- seedReceive[!sim$cohortData[age == 1L], on = c("pixelGroup", "speciesCode")]
+      #    (info contained within seedReceive)
+      # this is should be a inner join, needs to specify the nomatch=0, nomatch = NA is default that sugest the full joint.
+      seedSource <- seedSource[speciesCode %in% unique(seedReceive$speciesCode),]
+
+      # it could be more efficient if sim$pixelGroupMap is reduced map by removing the pixels that have
+      # successful postdisturbance regeneration and the inactive pixels
+      # how to subset the reducedmap
+      # if (sim$lastFireYear == round(time(sim))) { # the current year is both fire year and succession year
+      #   inactivePixelIndex <- c(sim$inactivePixelIndex, sim$treedFirePixelTableSinceLastDisp$pixelIndex)
+      # } else {
+      #   inactivePixelIndex <- sim$inactivePixelIndex
+      # }
+      reducedPixelGroupMap <- sim$pixelGroupMap
+
+      # Calculate the maximum size of the chunks for LANDISDisp
+      if (length(pixelsFromCurYrBurn) > 0) {
+        reducedPixelGroupMap[pixelsFromCurYrBurn] <- NA
+      }
 
-    if (getOption("LandR.verbose", TRUE) > 0) {
-      emptyForestPixels <- sim$treedFirePixelTableSinceLastDisp[burnTime < time(sim)]
-      seedsArrivedPixels <- unique(seedingData[emptyForestPixels, on = "pixelIndex", nomatch = 0], by = "pixelIndex")
-      message(blue("Of", NROW(emptyForestPixels),
-                   "burned and empty pixels: Num pixels where seeds arrived:",
-                   NROW(seedsArrivedPixels)))
-    }
+      seedingData <- LANDISDisp(dtRcv = seedReceive, plot.it = FALSE,
+                                dtSrc = seedSource,
+                                speciesTable = sim$species,
+                                pixelGroupMap = reducedPixelGroupMap,
+                                successionTimestep = P(sim)$successionTimestep,
+                                verbose = getOption("LandR.verbose", TRUE) > 0)
 
-    rm(seedReceive, seedSource)
-    if (NROW(seedingData) > 0) {
-      seedingData[, ecoregionGroup := factorValues2(sim$ecoregionMap, getValues(sim$ecoregionMap),
-                                                    att = "ecoregionGroup")[seedingData$pixelIndex]]
-      seedingData <- setkey(seedingData, ecoregionGroup, speciesCode)
-
-      specieseco_current <- speciesEcoregionLatestYear(
-        sim$speciesEcoregion[, .(year, speciesCode, establishprob, ecoregionGroup)],
-        round(time(sim)))
-      specieseco_current <- setkeyv(specieseco_current, c("ecoregionGroup", "speciesCode"))
-
-      # specieseco_current <- sim$speciesEcoregion[year <= round(time(sim))]
-      # specieseco_current <- setkey(specieseco_current[year == max(specieseco_current$year),
-      #                                                 .(speciesCode, establishprob, ecoregionGroup)],
-      #                              ecoregionGroup, speciesCode)
-      seedingData <- seedingData[specieseco_current, nomatch = 0]
-
-      ##############################################
-      # Run probability of establishment
-      ##############################################
-      LandR::assertCohortData(sim$cohortData, sim$pixelGroupMap, cohortDefinitionCols = P(sim)$cohortDefinitionCols)
-
-      seedingData <- seedingData[runif(nrow(seedingData)) <= establishprob, ]
       if (getOption("LandR.verbose", TRUE) > 0) {
-        seedsArrivedPixels <- unique(seedingData[emptyForestPixels, on = "pixelIndex", nomatch = 0],
-                                     by = "pixelIndex")
+        emptyForestPixels <- sim$treedFirePixelTableSinceLastDisp[burnTime <= time(sim)]
+        seedsArrivedPixels <- unique(seedingData[emptyForestPixels, on = "pixelIndex", nomatch = 0], by = "pixelIndex")
         message(blue("Of", NROW(emptyForestPixels),
-                     "burned and empty pixels: Num pixels where seedlings established:",
+                     "burned and empty pixels: Num pixels where seeds arrived:",
                      NROW(seedsArrivedPixels)))
       }
 
-      set(seedingData, NULL, "establishprob", NULL)
-      if (P(sim)$calibrate == TRUE) {
-        seedingData_summ <- seedingData[
-          , .(seedingAlgorithm = P(sim)$seedingAlgorithm, Year = round(time(sim)),
-              numberOfReg = length(pixelIndex)),
-          by = speciesCode]
-        seedingData_summ <- setkey(seedingData_summ, speciesCode)[
-          setkey(sim$species[, .(species,speciesCode)], speciesCode), nomatch = 0][
-            , .(species, seedingAlgorithm, Year, numberOfReg)]
-        sim$regenerationOutput <- rbindlist(list(sim$regenerationOutput, seedingData_summ))
-      }
-      if (nrow(seedingData) > 0) {
-        outs <- updateCohortData(seedingData, cohortData = sim$cohortData,
-                                 pixelGroupMap = sim$pixelGroupMap,
-                                 currentTime = round(time(sim)), speciesEcoregion = sim$speciesEcoregion,
-                                 treedFirePixelTableSinceLastDisp = NULL,
-                                 successionTimestep = P(sim)$successionTimestep)
-
-        sim$cohortData <- outs$cohortData
-        sim$pixelGroupMap <- outs$pixelGroupMap
+      rm(seedReceive, seedSource)
+      if (NROW(seedingData) > 0) {
+        seedingData[, ecoregionGroup := factorValues2(sim$ecoregionMap, getValues(sim$ecoregionMap),
+                                                      att = "ecoregionGroup")[seedingData$pixelIndex]]
+        seedingData <- setkey(seedingData, ecoregionGroup, speciesCode)
+
+        specieseco_current <- speciesEcoregionLatestYear(
+          sim$speciesEcoregion[, .(year, speciesCode, establishprob, ecoregionGroup)],
+          round(time(sim)))
+        specieseco_current <- setkeyv(specieseco_current, c("ecoregionGroup", "speciesCode"))
+
+        # specieseco_current <- sim$speciesEcoregion[year <= round(time(sim))]
+        # specieseco_current <- setkey(specieseco_current[year == max(specieseco_current$year),
+        #                                                 .(speciesCode, establishprob, ecoregionGroup)],
+        #                              ecoregionGroup, speciesCode)
+        seedingData <- seedingData[specieseco_current, nomatch = 0]
+
+        ##############################################
+        # Run probability of establishment
+        ##############################################
+        LandR::assertCohortData(sim$cohortData, sim$pixelGroupMap, cohortDefinitionCols = P(sim)$cohortDefinitionCols)
+
+        seedingData <- seedingData[runif(nrow(seedingData)) <= establishprob, ]
+        if (getOption("LandR.verbose", TRUE) > 0) {
+          seedsArrivedPixels <- unique(seedingData[emptyForestPixels, on = "pixelIndex", nomatch = 0],
+                                       by = "pixelIndex")
+          message(blue("Of", NROW(emptyForestPixels),
+                       "burned and empty pixels: Num pixels where seedlings established:",
+                       NROW(seedsArrivedPixels)))
+        }
+
+        set(seedingData, NULL, "establishprob", NULL)
+        if (P(sim)$calibrate == TRUE) {
+          seedingData_summ <- seedingData[
+            , .(seedingAlgorithm = P(sim)$seedingAlgorithm, Year = round(time(sim)),
+                numberOfReg = length(pixelIndex)),
+            by = speciesCode]
+          seedingData_summ <- setkey(seedingData_summ, speciesCode)[
+            setkey(sim$species[, .(species,speciesCode)], speciesCode), nomatch = 0][
+              , .(species, seedingAlgorithm, Year, numberOfReg)]
+          sim$regenerationOutput <- rbindlist(list(sim$regenerationOutput, seedingData_summ))
+        }
+        if (nrow(seedingData) > 0) {
+          outs <- updateCohortData(seedingData, cohortData = sim$cohortData,
+                                   pixelGroupMap = sim$pixelGroupMap,
+                                   currentTime = round(time(sim)), speciesEcoregion = sim$speciesEcoregion,
+                                   treedFirePixelTableSinceLastDisp = NULL,
+                                   successionTimestep = P(sim)$successionTimestep)
+
+          sim$cohortData <- outs$cohortData
+          sim$pixelGroupMap <- outs$pixelGroupMap
+        }
       }
     }
-  }
-  sim$lastReg <- round(time(sim))
-  return(invisible(sim))
-})
+    sim$lastReg <- round(time(sim))
+    return(invisible(sim))
+  }# )
 
 summaryRegen <- compiler::cmpfun(function(sim) {
   #cohortData <- sim$cohortData
@@ -1538,49 +1582,50 @@ summaryRegen <- compiler::cmpfun(function(sim) {
     if (!is.integer(pixelAll[["uniqueSumReproduction"]]))
       set(pixelAll, NULL, 'uniqueSumReproduction', asInteger(pixelAll[["uniqueSumReproduction"]]))
 
-    if (NROW(pixelAll) > 0) {
-      reproductionMap <- rasterizeReduced(pixelAll, pixelGroupMap, "uniqueSumReproduction")
-      setColors(reproductionMap) <- c("light green", "dark green")
-    } else {
-      reproductionMap <- setValues(pixelGroupMap, 0L)
+      if (NROW(pixelAll) > 0) {
+        reproductionMap <- rasterizeReduced(pixelAll, pixelGroupMap, "uniqueSumReproduction")
+        setColors(reproductionMap) <- c("light green", "dark green")
+      } else {
+        reproductionMap <- setValues(pixelGroupMap, 0L)
+      }
+      rm(pixelAll)
+      sim$reproductionMap <- reproductionMap
+      rm(pixelGroupMap)
     }
-    rm(pixelAll)
-    sim$reproductionMap <- reproductionMap
-    rm(pixelGroupMap)
+    return(invisible(sim))
   }
-  return(invisible(sim))
-})
+#)
 
 plotSummaryBySpecies <- compiler::cmpfun(function(sim) {
   LandR::assertSpeciesPlotLabels(sim$species$species, sim$sppEquiv)
 
-  checkPath(file.path(outputPath(sim), "figures"), create = TRUE)
+    checkPath(file.path(outputPath(sim), "figures"), create = TRUE)
 
-  ## BIOMASS, WEIGHTED AVERAGE AGE, AVERAGE ANPP
-  ## AND AGE OF OLDEST COHORT PER SPECIES
+    ## BIOMASS, WEIGHTED AVERAGE AGE, AVERAGE ANPP
+    ## AND AGE OF OLDEST COHORT PER SPECIES
 
-  ## Averages are calculated across pixels
-  ## don't expand table, multiply by no. pixels - faster
-  thisPeriod <- addNoPixel2CohortData(sim$cohortData, sim$pixelGroupMap, cohortDefinitionCols = P(sim)$cohortDefinitionCols)
+    ## Averages are calculated across pixels
+    ## don't expand table, multiply by no. pixels - faster
+    thisPeriod <- addNoPixel2CohortData(sim$cohortData, sim$pixelGroupMap, cohortDefinitionCols = P(sim)$cohortDefinitionCols)
 
-  for (column in names(thisPeriod)) if (is.integer(thisPeriod[[column]]))
-    set(thisPeriod, NULL, column, as.numeric(thisPeriod[[column]]))
+    for (column in names(thisPeriod)) if (is.integer(thisPeriod[[column]]))
+      set(thisPeriod, NULL, column, as.numeric(thisPeriod[[column]]))
 
-  thisPeriod <- thisPeriod[, list(year = time(sim),
-                                  BiomassBySpecies = sum(B * noPixels, na.rm = TRUE),
-                                  AgeBySppWeighted = sum(age * B * noPixels, na.rm = TRUE) /
-                                    sum(B * noPixels, na.rm = TRUE),
-                                  aNPPBySpecies = sum(aNPPAct * noPixels, na.rm = TRUE),
-                                  OldestCohortBySpp = max(age, na.rm = TRUE)),
-                           by = .(speciesCode)]
+    thisPeriod <- thisPeriod[, list(year = time(sim),
+                                    BiomassBySpecies = sum(B * noPixels, na.rm = TRUE),
+                                    AgeBySppWeighted = sum(age * B * noPixels, na.rm = TRUE) /
+                                      sum(B * noPixels, na.rm = TRUE),
+                                    aNPPBySpecies = sum(aNPPAct * noPixels, na.rm = TRUE),
+                                    OldestCohortBySpp = max(age, na.rm = TRUE)),
+                             by = .(speciesCode)]
 
-  #overstory
-  cohortData <-  addNoPixel2CohortData(sim$cohortData, sim$pixelGroupMap, cohortDefinitionCols = P(sim)$cohortDefinitionCols)
-  cohortData[, bWeightedAge := floor(sum(age * B) / sum(B) / 10) * 10, .(pixelGroup)]
-  # B was set as numeric to avoid problems with big numbers being integers
-  overstory <- cohortData[age >= bWeightedAge, .(overstoryBiomass = sum(as.numeric(B) * noPixels)),
-                          .(speciesCode)]
-  thisPeriod <- thisPeriod[overstory, on = 'speciesCode']
+    #overstory
+    cohortData <-  addNoPixel2CohortData(sim$cohortData, sim$pixelGroupMap, cohortDefinitionCols = P(sim)$cohortDefinitionCols)
+    cohortData[, bWeightedAge := floor(sum(age * B) / sum(B) / 10) * 10, .(pixelGroup)]
+    # B was set as numeric to avoid problems with big numbers being integers
+    overstory <- cohortData[age >= bWeightedAge, .(overstoryBiomass = sum(as.numeric(B) * noPixels)),
+                            .(speciesCode)]
+    thisPeriod <- thisPeriod[overstory, on = 'speciesCode']
 
   if (is.null(sim$summaryBySpecies)) {
     summaryBySpecies <- thisPeriod
@@ -1588,27 +1633,27 @@ plotSummaryBySpecies <- compiler::cmpfun(function(sim) {
     summaryBySpecies <- rbindlist(list(sim$summaryBySpecies, thisPeriod))
   }
 
-  ## MEAN NO. PIXELS PER LEADING SPECIES
-  vtm <- raster::mask(sim$vegTypeMap, sim$studyAreaReporting)
-  freqs <- table(na.omit(factorValues2(vtm, vtm[], att = 2)))
-  tabl <- as.vector(freqs)
-  summaryBySpecies1 <- data.frame(year = rep(floor(time(sim)), length(freqs)),
-                                  leadingType = names(freqs),
-                                  #freqs = freqs,
-                                  counts = tabl,
-                                  stringsAsFactors = FALSE)
-
-  whMixedLeading <- which(summaryBySpecies1$leadingType == "Mixed")
-  summaryBySpecies1$leadingType <- equivalentName(summaryBySpecies1$leadingType, sim$sppEquiv,
-                                                  "EN_generic_short")
-  summaryBySpecies1$leadingType[whMixedLeading] <- "Mixed"
-
-  colours <- equivalentName(names(sim$sppColorVect), sim$sppEquiv, "EN_generic_short")
-  whMixedSppColors <- which(names(sim$sppColorVect) == "Mixed")
-  colours[whMixedSppColors] <- "Mixed"
+    ## MEAN NO. PIXELS PER LEADING SPECIES
+    vtm <- raster::mask(sim$vegTypeMap, sim$studyAreaReporting)
+    freqs <- table(na.omit(factorValues2(vtm, vtm[], att = 2)))
+    tabl <- as.vector(freqs)
+    summaryBySpecies1 <- data.frame(year = rep(floor(time(sim)), length(freqs)),
+                                    leadingType = names(freqs),
+                                    #freqs = freqs,
+                                    counts = tabl,
+                                    stringsAsFactors = FALSE)
+
+    whMixedLeading <- which(summaryBySpecies1$leadingType == "Mixed")
+    summaryBySpecies1$leadingType <- equivalentName(summaryBySpecies1$leadingType, sim$sppEquiv,
+                                                    "EN_generic_short")
+    summaryBySpecies1$leadingType[whMixedLeading] <- "Mixed"
 
-  colorIDs <- match(summaryBySpecies1$leadingType, colours)
-  summaryBySpecies1$cols <- sim$sppColorVect[colorIDs]
+    colours <- equivalentName(names(sim$sppColorVect), sim$sppEquiv, "EN_generic_short")
+    whMixedSppColors <- which(names(sim$sppColorVect) == "Mixed")
+    colours[whMixedSppColors] <- "Mixed"
+
+    colorIDs <- match(summaryBySpecies1$leadingType, colours)
+    summaryBySpecies1$cols <- sim$sppColorVect[colorIDs]
 
   if (!is.null(sim$summaryBySpecies1)) {
     summaryBySpecies1 <- rbindlist(list(sim$summaryBySpecies1, summaryBySpecies1))
@@ -1800,36 +1845,50 @@ plotVegAttributesMaps <- compiler::cmpfun(function(sim) {
     grid.text(label = paste0("Year = ", round(time(sim))), x = 0.93, y = 0.97)
   }
 
-  return(invisible(sim))
-})
-
-plotAvgVegAttributes <- compiler::cmpfun(function(sim) {
-  LandR::assertSpeciesPlotLabels(sim$species$species, sim$sppEquiv)
-
-  checkPath(file.path(outputPath(sim), "figures"), create = TRUE)
-
-  ## AVERAGE STAND BIOMASS/AGE/ANPP
-  ## calculate acrosS pixels
-  ## don't expand table, multiply by no. pixels - faster
-  pixelCohortData <- addNoPixel2CohortData(sim$cohortData, sim$pixelGroupMap, cohortDefinitionCols = P(sim)$cohortDefinitionCols)
-  thisPeriod <- pixelCohortData[, list(year = time(sim),
-                                       sumB = sum(B*noPixels, na.rm = TRUE),
-                                       maxAge = asInteger(max(age, na.rm = TRUE)),
-                                       sumANPP = asInteger(sum(aNPPAct*noPixels, na.rm = TRUE)))]
-  denominator <- length(sim$pixelGroupMap[!is.na(sim$pixelGroupMap)]) * 100 #to get tonnes/ha
-  thisPeriod[, sumB := asInteger(sumB/denominator)]
-  thisPeriod[, sumANPP := asInteger(sumANPP/denominator)]
+      shortNames <- equivalentName(levsLeading, sppEquiv, "EN_generic_short")
+      shortNames[whMixedLevs] <- "Mixed"
+      levs[[levelsName]] <- shortNames
+      levels(sim$vegTypeMap) <- levs
+
+      colsLeading <- equivalentName(names(sim$sppColorVect), sppEquiv, "Leading")
+      colsLeading[whMixedSppColors] <- "Mixed"
+      sppColorVect <- sim$sppColorVect
+      names(sppColorVect) <- colsLeading
+      colours <- sppColorVect[na.omit(match(levsLeading, colsLeading))]
+      setColors(sim$vegTypeMap, levs$ID) <- colours
+
+      # Mask out NAs based on rasterToMatch (for plotting only!)
+      vegTypeMapForPlot <- raster::mask(sim$vegTypeMap, sim$studyAreaReporting)
+
+      ## Plot
+      tryCatch({
+        dev(mod$mapWindow) # Protecting from error of headless/terminal run
+        if (!is.null(biomassMapForPlot))
+          Plot(biomassMapForPlot, title = "Biomass", new = TRUE)
+        if (!is.null(ANPPMapForPlot))
+          Plot(ANPPMapForPlot, title = "ANPP", new = TRUE)
+        if (!is.null(mortalityMapForPlot))
+          Plot(mortalityMapForPlot, title = "Mortality", new = TRUE)
+        Plot(vegTypeMapForPlot, new = TRUE, title = "Leading vegetation")
+        grid.rect(0.93, 0.97, width = 0.2, height = 0.06, gp = gpar(fill = "white", col = "white"))
+        grid.text(label = paste0("Year = ", round(time(sim))), x = 0.93, y = 0.97)
+        #if (!is.null(reproductionMapForPlot))
+        #  Plot(reproductionMapForPlot, title = "Reproduction", new = TRUE)
+      }, error = function(e)
+        message("Can't open the device for plotting. Plotting will be disabled to avoid errors"))
+    }
 
   if (is.null(sim$summaryLandscape)) {
     summaryLandscape <- thisPeriod
   } else {
     summaryLandscape <- rbindlist(list(sim$summaryLandscape, thisPeriod))
   }
+#)
 
   if (length(unique(summaryLandscape$year)) > 1) {
     df2 <- melt(summaryLandscape, id.vars = "year")
 
-    varLabels <- c(sumB = "Biomass", maxAge = "Age", sumANPP = "aNPP")
+    checkPath(file.path(outputPath(sim), "figures"), create = TRUE)
 
     if (any(P(sim)$.plots == "screen")) {
       dev(mod$statsWindow)
@@ -1848,30 +1907,32 @@ plotAvgVegAttributes <- compiler::cmpfun(function(sim) {
   return(invisible(sim))
 })
 
-Save <- compiler::cmpfun(function(sim) {
-  raster::projection(sim$simulatedBiomassMap) <- raster::projection(sim$ecoregionMap)
-  raster::projection(sim$ANPPMap) <- raster::projection(sim$ecoregionMap)
-  raster::projection(sim$mortalityMap) <- raster::projection(sim$ecoregionMap)
-  raster::projection(sim$reproductionMap) <- raster::projection(sim$ecoregionMap)
-  writeRaster(sim$simulatedBiomassMap,
-              file.path(outputPath(sim), "figures",
-                        paste0("simulatedBiomassMap_Year", round(time(sim)), ".tif")),
-              datatype = 'INT4S', overwrite = TRUE)
-  writeRaster(sim$ANPPMap,
-              file.path(outputPath(sim), "figures",
-                        paste0("ANPP_Year", round(time(sim)), ".tif")),
-              datatype = 'INT4S', overwrite = TRUE)
-  writeRaster(sim$mortalityMap,
-              file.path(outputPath(sim), "figures",
-                        paste0("mortalityMap_Year", round(time(sim)), ".tif")),
-              datatype = 'INT4S', overwrite = TRUE)
-  writeRaster(sim$reproductionMap,
-              file.path(outputPath(sim), "figures",
-                        paste0("reproductionMap_Year", round(time(sim)), ".tif")),
-              datatype = 'INT4S', overwrite = TRUE)
+Save <- #compiler::cmpfun(
+  function(sim) {
+    raster::projection(sim$simulatedBiomassMap) <- raster::projection(sim$ecoregionMap)
+    raster::projection(sim$ANPPMap) <- raster::projection(sim$ecoregionMap)
+    raster::projection(sim$mortalityMap) <- raster::projection(sim$ecoregionMap)
+    raster::projection(sim$reproductionMap) <- raster::projection(sim$ecoregionMap)
+    writeRaster(sim$simulatedBiomassMap,
+                file.path(outputPath(sim), "figures",
+                          paste0("simulatedBiomassMap_Year", round(time(sim)), ".tif")),
+                datatype = 'INT4S', overwrite = TRUE)
+    writeRaster(sim$ANPPMap,
+                file.path(outputPath(sim), "figures",
+                          paste0("ANPP_Year", round(time(sim)), ".tif")),
+                datatype = 'INT4S', overwrite = TRUE)
+    writeRaster(sim$mortalityMap,
+                file.path(outputPath(sim), "figures",
+                          paste0("mortalityMap_Year", round(time(sim)), ".tif")),
+                datatype = 'INT4S', overwrite = TRUE)
+    writeRaster(sim$reproductionMap,
+                file.path(outputPath(sim), "figures",
+                          paste0("reproductionMap_Year", round(time(sim)), ".tif")),
+                datatype = 'INT4S', overwrite = TRUE)
 
-  return(invisible(sim))
-})
+    return(invisible(sim))
+  }
+#)
 
 CohortAgeReclassification <- function(sim) {
   if (time(sim) != start(sim)) {
@@ -1885,38 +1946,40 @@ CohortAgeReclassification <- function(sim) {
   }
 }
 
-.inputObjects <- compiler::cmpfun(function(sim) {
-  cacheTags <- c(currentModule(sim), "function:.inputObjects")
-  dPath <- asPath(getOption("reproducible.destinationPath", dataPath(sim)), 1)
-  if (getOption("LandR.verbose", TRUE) > 0)
-    message(currentModule(sim), ": using dataPath '", dPath, "'.")
+.inputObjects <- #compiler::cmpfun(
+  function(sim) {
+    cacheTags <- c(currentModule(sim), "function:.inputObjects")
+    dPath <- asPath(getOption("reproducible.destinationPath", dataPath(sim)), 1)
+    if (getOption("LandR.verbose", TRUE) > 0)
+      message(currentModule(sim), ": using dataPath '", dPath, "'.")
 
-  #######################################################
+    #######################################################
 
-  if (!suppliedElsewhere("studyArea", sim)) {
-    stop("Please provide a 'studyArea' polygon")
-    # message("'studyArea' was not provided by user. Using a polygon (6250000 m^2) in southwestern Alberta, Canada")
-    # sim$studyArea <- randomStudyArea(seed = 1234, size = (250^2)*100)  # Jan 2021 we agreed to force user to provide a SA/SAL
-  }
+    if (!suppliedElsewhere("studyArea", sim)) {
+      sim$studyArea <- randomStudyArea(seed = 1234, size = (250^2)*100)  # Jan 2021 we agreed to force user to provide a SA/SAL
+      # stop("Please provide a 'studyArea' polygon")
+      # message("'studyArea' was not provided by user. Using a polygon (6250000 m^2) in southwestern Alberta, Canada")
+      # sim$studyArea <- randomStudyArea(seed = 1234, size = (250^2)*100)  # Jan 2021 we agreed to force user to provide a SA/SAL
+    }
 
-  if (is.na(P(sim)$.studyAreaName)) {
-    params(sim)[[currentModule(sim)]][[".studyAreaName"]] <- reproducible::studyAreaName(sim$studyArea)
-    message("The .studyAreaName is not supplied; derived name from sim$studyArea: ",
-            params(sim)[[currentModule(sim)]][[".studyAreaName"]])
-  }
+    if (is.na(P(sim)$.studyAreaName)) {
+      params(sim)[[currentModule(sim)]][[".studyAreaName"]] <- reproducible::studyAreaName(sim$studyArea)
+      message("The .studyAreaName is not supplied; derived name from sim$studyArea: ",
+              params(sim)[[currentModule(sim)]][[".studyAreaName"]])
+    }
 
-  needRTM <- FALSE
-  if (is.null(sim$rasterToMatch)) {
-    if (!suppliedElsewhere("rasterToMatch", sim)) {
-      needRTM <- TRUE
-      message("There is no rasterToMatch supplied; will attempt to use rawBiomassMap")
-    } else {
-      stop("rasterToMatch is going to be supplied, but ", currentModule(sim), " requires it ",
-           "as part of its .inputObjects. Please make it accessible to ", currentModule(sim),
-           " in the .inputObjects by passing it in as an object in simInit(objects = list(rasterToMatch = aRaster)",
-           " or in a module that gets loaded prior to ", currentModule(sim))
+    needRTM <- FALSE
+    if (is.null(sim$rasterToMatch)) {
+      if (!suppliedElsewhere("rasterToMatch", sim)) {
+        needRTM <- TRUE
+        message("There is no rasterToMatch supplied; will attempt to use rawBiomassMap")
+      } else {
+        stop("rasterToMatch is going to be supplied, but ", currentModule(sim), " requires it ",
+             "as part of its .inputObjects. Please make it accessible to ", currentModule(sim),
+             " in the .inputObjects by passing it in as an object in simInit(objects = list(rasterToMatch = aRaster)",
+             " or in a module that gets loaded prior to ", currentModule(sim))
+      }
     }
-  }
 
   if (needRTM) {
     if (!suppliedElsewhere("rawBiomassMap", sim) ||
@@ -1956,154 +2019,154 @@ CohortAgeReclassification <- function(sim) {
                              omitArgs = c("destinationPath", "targetFile", "userTags", "stable"))
     }
 
-    ## if we need rasterToMatchLarge, that means a) we don't have it, but b) we will have rawBiomassMap
-
-    warning("rasterToMatch is missing and will be created",
-            " from rawBiomassMap and studyAreaLarge.",
-            " If this is wrong, provide raster.")
-
-    sim$rasterToMatch <- rawBiomassMap
-    RTMvals <- getValues(sim$rasterToMatch)
-    sim$rasterToMatch[!is.na(RTMvals)] <- 1
-    sim$rasterToMatch <- Cache(writeOutputs, sim$rasterToMatch,
-                               filename2 = .suffix(file.path(dPath, "rasterToMatch.tif"),
-                                                   paste0("_", P(sim)$.studyAreaName)),
-                               datatype = "INT2U", overwrite = TRUE,
-                               userTags = c(cacheTags, "rasterToMatch"),
-                               omitArgs = c("userTags"))
-  }
+      ## if we need rasterToMatchLarge, that means a) we don't have it, but b) we will have rawBiomassMap
+
+      warning("rasterToMatch is missing and will be created",
+              " from rawBiomassMap and studyAreaLarge.",
+              " If this is wrong, provide raster.")
+
+      sim$rasterToMatch <- rawBiomassMap
+      RTMvals <- getValues(sim$rasterToMatch)
+      sim$rasterToMatch[!is.na(RTMvals)] <- 1
+      sim$rasterToMatch <- Cache(writeOutputs, sim$rasterToMatch,
+                                 filename2 = .suffix(file.path(dPath, "rasterToMatch.tif"),
+                                                     paste0("_", P(sim)$.studyAreaName)),
+                                 datatype = "INT2U", overwrite = TRUE,
+                                 userTags = c(cacheTags, "rasterToMatch"),
+                                 omitArgs = c("userTags"))
+    }
 
-  if (!compareCRS(sim$studyArea, sim$rasterToMatch)) {
-    warning(paste0("studyArea and rasterToMatch projections differ.\n",
-                   "studyArea will be projected to match rasterToMatch"))
-    sim$studyArea <- spTransform(sim$studyArea, crs(sim$rasterToMatch))
-    sim$studyArea <- fixErrors(sim$studyArea)
-  }
+    if (!compareCRS(sim$studyArea, sim$rasterToMatch)) {
+      warning(paste0("studyArea and rasterToMatch projections differ.\n",
+                     "studyArea will be projected to match rasterToMatch"))
+      sim$studyArea <- spTransform(sim$studyArea, crs(sim$rasterToMatch))
+      sim$studyArea <- fixErrors(sim$studyArea)
+    }
 
-  if (!suppliedElsewhere("studyAreaReporting", sim)) {
-    if (getOption("LandR.verbose", TRUE) > 0)
-      message("'studyAreaReporting' was not provided by user. Using the same as 'studyArea'.")
-    sim$studyAreaReporting <- sim$studyArea
-  }
+    if (!suppliedElsewhere("studyAreaReporting", sim)) {
+      if (getOption("LandR.verbose", TRUE) > 0)
+        message("'studyAreaReporting' was not provided by user. Using the same as 'studyArea'.")
+      sim$studyAreaReporting <- sim$studyArea
+    }
 
-  ## make light requirements table
-  if (!suppliedElsewhere("sufficientLight", sim)) {
-    ## load the biomass_succession.txt to get shade tolerance parameters
-    mainInput <- prepInputsMainInput(url = extractURL("sufficientLight"),
-                                     dPath,
-                                     cacheTags = c(cacheTags, "mainInput")) ## uses default URL
-
-    sufficientLight <- data.frame(mainInput, stringsAsFactors = FALSE)
-    startRow <- which(sufficientLight$col1 == "SufficientLight")
-    sufficientLight <- sufficientLight[(startRow + 1):(startRow + 5), 1:7]
-    sufficientLight <- data.table(sufficientLight)
-    sufficientLight <- sufficientLight[, lapply(.SD, function(x) as.numeric(x))]
-
-    names(sufficientLight) <- c("speciesshadetolerance",
-                                "X0", "X1", "X2", "X3", "X4", "X5")
-    sim$sufficientLight <- data.frame(sufficientLight, stringsAsFactors = FALSE)
-  }
+    ## make light requirements table
+    if (!suppliedElsewhere("sufficientLight", sim)) {
+      ## load the biomass_succession.txt to get shade tolerance parameters
+      mainInput <- prepInputsMainInput(url = extractURL("sufficientLight"),
+                                       dPath,
+                                       cacheTags = c(cacheTags, "mainInput")) ## uses default URL
+
+      sufficientLight <- data.frame(mainInput, stringsAsFactors = FALSE)
+      startRow <- which(sufficientLight$col1 == "SufficientLight")
+      sufficientLight <- sufficientLight[(startRow + 1):(startRow + 5), 1:7]
+      sufficientLight <- data.table(sufficientLight)
+      sufficientLight <- sufficientLight[, lapply(.SD, function(x) as.numeric(x))]
+
+      names(sufficientLight) <- c("speciesshadetolerance",
+                                  "X0", "X1", "X2", "X3", "X4", "X5")
+      sim$sufficientLight <- data.frame(sufficientLight, stringsAsFactors = FALSE)
+    }
 
-  if (!suppliedElsewhere("sppEquiv", sim)) {
-    if (!is.null(sim$sppColorVect))
-      stop("If you provide sppColorVect, you MUST also provide sppEquiv")
+    if (!suppliedElsewhere("sppEquiv", sim)) {
+      if (!is.null(sim$sppColorVect))
+        stop("If you provide sppColorVect, you MUST also provide sppEquiv")
 
-    data("sppEquivalencies_CA", package = "LandR", envir = environment())
-    sim$sppEquiv <- as.data.table(sppEquivalencies_CA)
-    ## By default, Abies_las is renamed to Abies_sp
-    sim$sppEquiv[KNN == "Abie_Las", LandR := "Abie_sp"]
+      data("sppEquivalencies_CA", package = "LandR", envir = environment())
+      sim$sppEquiv <- as.data.table(sppEquivalencies_CA)
+      ## By default, Abies_las is renamed to Abies_sp
+      sim$sppEquiv[KNN == "Abie_Las", LandR := "Abie_sp"]
 
-    ## check spp column to use
-    if (P(sim)$sppEquivCol == "Boreal") {
-      message(paste("There is no 'sppEquiv' table supplied;",
-                    "will attempt to use species listed under 'Boreal'",
-                    "in the 'LandR::sppEquivalencies_CA' table"))
-    } else {
-      if (grepl(P(sim)$sppEquivCol, names(sim$sppEquiv))) {
-        message(paste("There is no 'sppEquiv' table supplied,",
-                      "will attempt to use species listed under", P(sim)$sppEquivCol,
+      ## check spp column to use
+      if (P(sim)$sppEquivCol == "Boreal") {
+        message(paste("There is no 'sppEquiv' table supplied;",
+                      "will attempt to use species listed under 'Boreal'",
                       "in the 'LandR::sppEquivalencies_CA' table"))
       } else {
-        stop("You changed 'sppEquivCol' without providing 'sppEquiv',",
-             "and the column name can't be found in the default table ('LandR::sppEquivalencies_CA').",
-             "Please provide conforming 'sppEquivCol', 'sppEquiv' and 'sppColorVect'")
+        if (grepl(P(sim)$sppEquivCol, names(sim$sppEquiv))) {
+          message(paste("There is no 'sppEquiv' table supplied,",
+                        "will attempt to use species listed under", P(sim)$sppEquivCol,
+                        "in the 'LandR::sppEquivalencies_CA' table"))
+        } else {
+          stop("You changed 'sppEquivCol' without providing 'sppEquiv',",
+               "and the column name can't be found in the default table ('LandR::sppEquivalencies_CA').",
+               "Please provide conforming 'sppEquivCol', 'sppEquiv' and 'sppColorVect'")
+        }
       }
-    }
 
-    ## remove empty lines/NAs
-    sim$sppEquiv <- sim$sppEquiv[!"", on = P(sim)$sppEquivCol]
-    sim$sppEquiv <- na.omit(sim$sppEquiv, P(sim)$sppEquivCol)
+      ## remove empty lines/NAs
+      sim$sppEquiv <- sim$sppEquiv[!"", on = P(sim)$sppEquivCol]
+      sim$sppEquiv <- na.omit(sim$sppEquiv, P(sim)$sppEquivCol)
 
-    ## add default colors for species used in model
-    sim$sppColorVect <- sppColors(sim$sppEquiv, P(sim)$sppEquivCol,
-                                  newVals = "Mixed", palette = "Accent")
-  } else {
-    if (is.null(sim$sppColorVect)) {
-      message("'sppEquiv' is provided without a 'sppColorVect'. Running:
-              LandR::sppColors with column ", P(sim)$sppEquivCol)
+      ## add default colors for species used in model
       sim$sppColorVect <- sppColors(sim$sppEquiv, P(sim)$sppEquivCol,
                                     newVals = "Mixed", palette = "Accent")
+    } else {
+      if (is.null(sim$sppColorVect)) {
+        message("'sppEquiv' is provided without a 'sppColorVect'. Running:
+              LandR::sppColors with column ", P(sim)$sppEquivCol)
+        sim$sppColorVect <- sppColors(sim$sppEquiv, P(sim)$sppEquivCol,
+                                      newVals = "Mixed", palette = "Accent")
+      }
     }
-  }
-
-  if (P(sim)$vegLeadingProportion > 0 & is.na(sim$sppColorVect['Mixed'])) {
-    stop("vegLeadingProportion  is > 0 but there is no 'Mixed' color in sim$sppColorVect. ",
-         "Please supply sim$sppColorVect with a 'Mixed' color or set vegLeadingProportion to zero.")
-  }
 
+    if (P(sim)$vegLeadingProportion > 0 & is.na(sim$sppColorVect['Mixed'])) {
+      stop("vegLeadingProportion  is > 0 but there is no 'Mixed' color in sim$sppColorVect. ",
+           "Please supply sim$sppColorVect with a 'Mixed' color or set vegLeadingProportion to zero.")
+    }
 
-  if (!suppliedElsewhere("treedFirePixelTableSinceLastDisp", sim)) {
-    sim$treedFirePixelTableSinceLastDisp <- data.table(pixelIndex = integer(),
-                                                       pixelGroup = integer(),
-                                                       burnTime = numeric())
-  }
 
-  if (!suppliedElsewhere("speciesLayers", sim)) {
-    message("No RasterStack map of biomass X species is provided; using KNN")
-    url <- paste0("http://ftp.maps.canada.ca/pub/nrcan_rncan/Forests_Foret/",
-                  "canada-forests-attributes_attributs-forests-canada/2001-attributes_attributs-2001/")
-    sim$speciesLayers <- Cache(loadkNNSpeciesLayers,
-                               dPath = dPath,
-                               rasterToMatch = sim$rasterToMatch,
-                               studyArea = sim$studyArea,
-                               sppEquiv = sim$sppEquiv,
-                               knnNamesCol = "KNN",
-                               sppEquivCol = P(sim)$sppEquivCol,
-                               thresh = 10,
-                               url = url,
-                               userTags = c(cacheTags, "speciesLayers"),
-                               omitArgs = c("userTags"))
-  }
+    if (!suppliedElsewhere("treedFirePixelTableSinceLastDisp", sim)) {
+      sim$treedFirePixelTableSinceLastDisp <- data.table(pixelIndex = integer(),
+                                                         pixelGroup = integer(),
+                                                         burnTime = numeric())
+    }
 
-  ## additional species traits
-  if (!suppliedElsewhere("species", sim)) {
-    speciesTable <- getSpeciesTable(dPath = dPath, url = extractURL("species"),
-                                    cacheTags = c(cacheTags, "speciesTable"))
-    sim$species <- prepSpeciesTable(speciesTable = speciesTable,
-                                    # speciesLayers = sim$speciesLayers,
-                                    sppEquiv = sim$sppEquiv[get(P(sim)$sppEquivCol) %in%
-                                                              names(sim$speciesLayers)],
-                                    sppEquivCol = P(sim)$sppEquivCol)
-  }
+    if (!suppliedElsewhere("speciesLayers", sim)) {
+      message("No RasterStack map of biomass X species is provided; using KNN")
+      url <- paste0("http://ftp.maps.canada.ca/pub/nrcan_rncan/Forests_Foret/",
+                    "canada-forests-attributes_attributs-forests-canada/2001-attributes_attributs-2001/")
+      sim$speciesLayers <- Cache(loadkNNSpeciesLayers,
+                                 dPath = dPath,
+                                 rasterToMatch = sim$rasterToMatch,
+                                 studyArea = sim$studyArea,
+                                 sppEquiv = sim$sppEquiv,
+                                 knnNamesCol = "KNN",
+                                 sppEquivCol = P(sim)$sppEquivCol,
+                                 thresh = 10,
+                                 url = url,
+                                 userTags = c(cacheTags, "speciesLayers"),
+                                 omitArgs = c("userTags"))
+    }
 
-  if (P(sim)$growthAndMortalityDrivers != 'LandR') {
-    if (!suppliedElsewhere("cceArgs", sim)) {
-      sim$cceArgs <- list(quote(CMI),
-                          quote(ATA),
-                          quote(CMInormal),
-                          quote(mcsModel),
-                          quote(gcsModel))
-      names(sim$cceArgs) <- paste(sim$cceArgs)
+    ## additional species traits
+    if (!suppliedElsewhere("species", sim)) {
+      speciesTable <- getSpeciesTable(dPath = dPath, url = extractURL("species"),
+                                      cacheTags = c(cacheTags, "speciesTable"))
+      sim$species <- prepSpeciesTable(speciesTable = speciesTable,
+                                      # speciesLayers = sim$speciesLayers,
+                                      sppEquiv = sim$sppEquiv[get(P(sim)$sppEquivCol) %in%
+                                                                names(sim$speciesLayers)],
+                                      sppEquivCol = P(sim)$sppEquivCol)
     }
 
-    #check for climate args
-    # if (!all(unlist(lapply(names(sim$cceArgs), suppliedElsewhere, sim = sim)))) {
-    #   stop("Some or all of sim$cceArgs are not supplied")
-    # }
-  }
+    if (P(sim)$growthAndMortalityDrivers != 'LandR') {
+      if (!suppliedElsewhere("cceArgs", sim)) {
+        sim$cceArgs <- list(quote(CMI),
+                            quote(ATA),
+                            quote(CMInormal),
+                            quote(mcsModel),
+                            quote(gcsModel))
+        names(sim$cceArgs) <- paste(sim$cceArgs)
+      }
 
+      #check for climate args
+      # if (!all(unlist(lapply(names(sim$cceArgs), suppliedElsewhere, sim = sim)))) {
+      #   stop("Some or all of sim$cceArgs are not supplied")
+      # }
+    }
 
-  gc() ## AMC added this 2019-08-20
+    # gc() ## AMC added this 2019-08-20
 
-  return(invisible(sim))
-})
+    return(invisible(sim))
+  }
+#)
diff --git a/Biomass_core.Rmd b/Biomass_core.Rmd
index 7b1e277..df5c8d8 100644
--- a/Biomass_core.Rmd
+++ b/Biomass_core.Rmd
@@ -160,11 +160,13 @@ setPaths(inputPath = file.path(tempDir, "inputs"),
 
 times <- list(start = 0, end = 10)
 
-studyArea <- Cache(randomStudyArea, size = 1e7) # cache this so it creates a random one only once on a machine
+studyArea <- Cache(randomStudyArea, size = 1e10) # cache this so it creates a random one only once on a machine
 
+set.seed(123)
 # Pick the species you want to work with -- using the naming convention in "Boreal" column of LandR::sppEquivalencies_CA
 speciesNameConvention <- "Boreal"
 speciesToUse <- c("Pice_Gla", "Popu_Tre", "Pinu_Con")
+#speciesToUse <- c("Pice_Gla")
 
 sppEquiv <- LandR::sppEquivalencies_CA[get(speciesNameConvention) %in% speciesToUse]
 # Assign a colour convention for graphics for each species
diff --git a/R/age-cohorts.R b/R/age-cohorts.R
index 003b625..96e994b 100644
--- a/R/age-cohorts.R
+++ b/R/age-cohorts.R
@@ -1,4 +1,5 @@
-ageReclassification <- compiler::cmpfun(function(cohortData, successionTimestep, stage, byGroups = c('pixelGroup', 'speciesCode', 'age')) {
+ageReclassification <- # ompiler::cmpfun(
+  function(cohortData, successionTimestep, stage, byGroups = c('pixelGroup', 'speciesCode', 'age')) {
 
   byGroupsNoAge <- byGroups[!byGroups %in% 'age'] #age is what will be lumped
   #byGroups default added for backwards compatibility
@@ -23,9 +24,8 @@ ageReclassification <- compiler::cmpfun(function(cohortData, successionTimestep,
     cdColNames <- colnames(cohortData)
     message("  Setting all ages <= ", successionTimestep, " to ", successionTimestepPlusOne)
     if (any(anyDuplicates)) {
-
       # pull out only duplicated types -- note "which = TRUE" gives only the indices of the joined rows -- will use the inverse below
-      tdDuplicates <- targetData[targetData[anyDuplicates], nomatch = NULL,
+      tdDuplicates <- targetData[unique(targetData[anyDuplicates], by = byGroupsNoAge), nomatch = NULL,
                                  on = byGroupsNoAge, which = TRUE]
 
       td <- targetData[tdDuplicates]
@@ -57,4 +57,4 @@ ageReclassification <- compiler::cmpfun(function(cohortData, successionTimestep,
     }
   }
   return(cohortData)
-})
+}#)
diff --git a/R/helpers.R b/R/helpers.R
index b622b92..7a26b58 100644
--- a/R/helpers.R
+++ b/R/helpers.R
@@ -38,6 +38,11 @@ updateSpeciesEcoregionAttributes <- function(speciesEcoregion, currentTime, coho
 #' @importFrom data.table setkey
 updateSpeciesAttributes <- function(species, cohortData) {
   # to assign longevity, mortalityshape, growthcurve to cohortData
+  if (!("speciesCode" %in% colnames(species)) | !("speciesCode" %in% colnames(cohortData))) {
+    browser()
+    species <- species[, speciesCode := as.factor(species)]
+  }
+
   species_temp <- setkey(species[, .(speciesCode, longevity, mortalityshape, growthcurve)], speciesCode)
   setkey(cohortData, speciesCode)
   cohortData <- cohortData[species_temp, nomatch = 0]
@@ -60,7 +65,8 @@ updateSpeciesAttributes <- function(species, cohortData) {
 #'
 #' @export
 #' @importFrom data.table copy rbindlist setkey
-calculateSumB <- compiler::cmpfun(function(cohortData, lastReg, currentTime, successionTimestep,
+calculateSumB2 <- #compiler::cmpfun(
+  function(cohortData, lastReg, currentTime, successionTimestep,
                                            doAssertion = getOption("LandR.assertions", TRUE)) {
   nrowCohortData <- NROW(cohortData)
 
@@ -123,7 +129,8 @@ calculateSumB <- compiler::cmpfun(function(cohortData, lastReg, currentTime, suc
   if (algo == 1 || isTRUE(doAssertion)) {
     ## this older version is typically much slower than the newer one below (Eliot June 2, 2019)
     cohortData1 <- copy(cohortData)
-    old1 <- Sys.time()
+    if (isTRUE(doAssertion))
+      old1 <- Sys.time()
     oldKey <- checkAndChangeKey(cohortData1, "pixelGroup")
     cohortData1[age >= successionTimestep, sumB := sum(B, na.rm = TRUE), by = "pixelGroup"]
     setorderv(cohortData1, c("sumB"), na.last = TRUE)
@@ -136,7 +143,8 @@ calculateSumB <- compiler::cmpfun(function(cohortData, lastReg, currentTime, suc
     set(cohortData1, NULL, "sumB2", NULL)
     if (!is.null(oldKey))
       setkeyv(cohortData1, oldKey)
-    old2 <- Sys.time()
+    if (isTRUE(doAssertion))
+      old2 <- Sys.time()
     if (!is.integer(cohortData1[["sumB"]]))
       set(cohortData1, NULL, "sumB", asInteger(cohortData1[["sumB"]]))
   }
@@ -177,15 +185,14 @@ calculateSumB <- compiler::cmpfun(function(cohortData, lastReg, currentTime, suc
   }
 
   cohortData <- if (algo == 1) copy(cohortData1) else copy(cohortData2)
-
   if (isTRUE(doAssertion)) {
-    if (!exists("oldAlgoSumB")) mod$oldAlgoSumB <- 0
-    if (!exists("newAlgoSumB")) mod$newAlgoSumB <- 0
+    if (!exists("oldAlgoSumB", envir = mod)) mod$oldAlgoSumB <- 0
+    if (!exists("newAlgoSumB", envir = mod)) mod$newAlgoSumB <- 0
     mod$oldAlgoSumB <- mod$oldAlgoSumB + (old2 - old1)
     mod$newAlgoSumB <- mod$newAlgoSumB + (new2 - new1)
 
-    print(paste("Biomass_core:calculateSumB: new algo", mod$newAlgoSumB))
-    print(paste("Biomass_core:calculateSumB: old algo", mod$oldAlgoSumB))
+    print(paste("Biomass_core:calculateSumB: new algo (cumulative)", mod$newAlgoSumB))
+    print(paste("Biomass_core:calculateSumB: old algo (cumulative)", mod$oldAlgoSumB))
 
     setkeyv(cohortData, c("pixelGroup", "speciesCode", "age"))
     setkeyv(cohortData1, c("pixelGroup", "speciesCode", "age"))
@@ -202,11 +209,11 @@ calculateSumB <- compiler::cmpfun(function(cohortData, lastReg, currentTime, suc
     #}
   }
   return(cohortData)
-})
+}#)
 
 #' calculateAgeMortality
 #'
-#' TODO: description and title needed
+#' Follows https://github.com/LANDIS-II-Foundation/Extension-Biomass-Succession/blob/c07f044c4ffaeb075c714b3c0a23e3c22761e7ab/src/CohortBiomass.cs#L141
 #'
 #' @param cohortData \code{data.table} TODO: description needed
 #' @param stage TODO: description needed
@@ -234,7 +241,7 @@ calculateAgeMortality <- function(cohortData, stage = "nonSpinup", spinupMortali
 
 #' calculateANPP
 #'
-#' TODO: description and title needed
+#' Follows https://github.com/LANDIS-II-Foundation/Extension-Biomass-Succession/blob/c07f044c4ffaeb075c714b3c0a23e3c22761e7ab/src/CohortBiomass.cs#L163
 #'
 #' @param cohortData \code{data.table} TODO: description needed
 #' @param stage TODO: description needed
@@ -243,7 +250,8 @@ calculateAgeMortality <- function(cohortData, stage = "nonSpinup", spinupMortali
 #'
 #' @export
 #' @importFrom data.table set
-calculateANPP <- compiler::cmpfun(function(cohortData, stage = "nonSpinup") {
+calculateANPP <- # compiler::cmpfun(
+  function(cohortData, stage = "nonSpinup") {
   if (stage == "spinup") {
     cohortData[age > 0, aNPPAct := maxANPP * exp(1) * (bAP^growthcurve) *
                  exp(-(bAP^growthcurve)) * bPM]
@@ -261,11 +269,11 @@ calculateANPP <- compiler::cmpfun(function(cohortData, stage = "nonSpinup") {
         pmin(cohortData$maxANPP*cohortData$bPM, aNPPAct))
   }
   return(cohortData)
-})
+}# )
 
 #' calculateGrowthMortality
 #'
-#' TODO: description and title needed
+#' Follows https://github.com/LANDIS-II-Foundation/Extension-Biomass-Succession/blob/c07f044c4ffaeb075c714b3c0a23e3c22761e7ab/src/CohortBiomass.cs#L243
 #'
 #' @param cohortData \code{data.table} TODO: description needed
 #' @param stage TODO: description needed
@@ -275,7 +283,8 @@ calculateANPP <- compiler::cmpfun(function(cohortData, stage = "nonSpinup") {
 #' @export
 #' @importFrom data.table set
 #' @importFrom fpCompare %>>% %<=%
-calculateGrowthMortality <- compiler::cmpfun(function(cohortData, stage = "nonSpinup") {
+calculateGrowthMortality <- #compiler::cmpfun(
+  function(cohortData, stage = "nonSpinup") {
   if (stage == "spinup") {
     cohortData[age > 0 & bAP %>>% 1.0, mBio := maxANPP*bPM]
     cohortData[age > 0 & bAP %<=% 1.0, mBio := maxANPP*(2*bAP) / (1 + bAP)*bPM]
@@ -290,10 +299,11 @@ calculateGrowthMortality <- compiler::cmpfun(function(cohortData, stage = "nonSp
         pmin(cohortData$maxANPP*cohortData$bPM, cohortData$mBio))
   }
   return(cohortData)
-})
+}#)
 
 #' calculateCompetition
 #'
+#' Follows https://github.com/LANDIS-II-Foundation/Extension-Biomass-Succession/blob/c07f044c4ffaeb075c714b3c0a23e3c22761e7ab/src/CohortBiomass.cs#L437
 #' TODO: description and title needed
 #'
 #' @param cohortData \code{data.table} TODO: description needed
@@ -303,7 +313,8 @@ calculateGrowthMortality <- compiler::cmpfun(function(cohortData, stage = "nonSp
 #'
 #' @export
 #' @importFrom data.table key setkeyv
-calculateCompetition <- compiler::cmpfun(function(cohortData, stage = "nonSpinup") {
+calculateCompetition <- #compiler::cmpfun(
+  function(cohortData, stage = "nonSpinup") {
   # two competition indics are calculated bAP and bPM
   if (stage == "spinup") {
     cohortData[age > 0, bPot := pmax(1, maxB - sumB + B)]
@@ -339,7 +350,7 @@ calculateCompetition <- compiler::cmpfun(function(cohortData, stage = "nonSpinup
     set(cohortData, NULL, c("cMultiplier"), NULL)
   }
   return(cohortData)
-})
+}#)
 
 checkAndChangeKey <- function(obj, key) {
   oldKey <- key(obj)
diff --git a/R/spinup.R b/R/spinup.R
index 969eb04..7f7fce9 100644
--- a/R/spinup.R
+++ b/R/spinup.R
@@ -12,7 +12,8 @@
 #' @export
 #' @importFrom data.table data.table set setkey
 #'
-spinUp <- compiler::cmpfun(function(cohortData, calibrate, successionTimestep,
+spinUp <- #compiler::cmpfun(
+  function(cohortData, calibrate, successionTimestep,
                                     spinupMortalityfraction, species) {
   maxAge <- max(cohortData$age, na.rm = TRUE) # determine the pre-simulation length
   set(cohortData, NULL, "origAge", cohortData$age)
@@ -113,7 +114,7 @@ spinUp <- compiler::cmpfun(function(cohortData, calibrate, successionTimestep,
     all <- list(cohortData = cohortData)
   }
   return(all)
-})
+}#)
 
 # cacheSpinUpFunction <- function(sim, cachePath) {
 #   # for slow functions, add cached versions. Then use sim$xxx() throughout module instead of xxx()
diff --git a/tests/testthat/test-Biomass_coreFireDisturbance.R b/tests/testthat/test-Biomass_coreFireDisturbance.R
index 296dae5..b0ffd5e 100644
--- a/tests/testthat/test-Biomass_coreFireDisturbance.R
+++ b/tests/testthat/test-Biomass_coreFireDisturbance.R
@@ -1,14 +1,21 @@
 test_that("test regeneration after fire. ",{
-  library(SpaDES)
-  library(data.table)
-  library(raster)
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
+
   pixelGroupMap <- raster(xmn = 50,xmx = 50 + 50*100,
                           ymn = 50,ymx = 50 + 50*100,
                           res = c(100, 100),
diff --git a/tests/testthat/test-Biomass_coreInit.R b/tests/testthat/test-Biomass_coreInit.R
index 939ba69..3ad5280 100644
--- a/tests/testthat/test-Biomass_coreInit.R
+++ b/tests/testthat/test-Biomass_coreInit.R
@@ -1,17 +1,21 @@
-test_that("test Biomass_coreInit", {
-  # define the module and path
-  library(raster)
-  library(data.table)
-  module <- list("Biomass_core")
-  path <- list(
-    modulePath = "..",
-    outputPath = "~/output"
-  )
-  parameters <- list(
-    .progress = list(type = "graphical", interval = 1),
-    .globals = list(verbose = FALSE),
-    Biomass_core = list(.saveInitialTime = NA)
-  )
+test_that("test Init", {
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
+
   initialCommunitiesMap <- raster(
     xmn = 50, xmx = 50 + 3 * 100,
     ymn = 50, ymx = 50 + 3 * 100,
@@ -77,10 +81,10 @@ test_that("test Biomass_coreInit", {
     objects = objects,
     paths = path
   )
-  if (exists("Biomass_coreInit")) {
-    simOutput <- Biomass_coreInit(mySim)
+  if (exists("Init")) {
+    simOutput <- Init(mySim)
   } else {
-    simOutput <- mySim$.mods$Biomass_core$Biomass_coreInit(mySim)
+    simOutput <- mySim$.mods$Biomass_core$Init(mySim)
   }
   # check the cohortData table
   cohortData <- simOutput$cohortData
@@ -179,10 +183,10 @@ test_that("test Biomass_coreInit", {
     objects = objects,
     paths = path
   )
-  if (exists("Biomass_coreInit")) {
-    simOutput <- Biomass_coreInit(mySim)
+  if (exists("Init")) {
+    simOutput <- Init(mySim)
   } else {
-    simOutput <- mySim$.mods$Biomass_core$Biomass_coreInit(mySim)
+    simOutput <- mySim$.mods$Biomass_core$Init(mySim)
   }
 
   # check the calibration mode
diff --git a/tests/testthat/test-Biomass_coreMortalityAndGrowth.R b/tests/testthat/test-Biomass_coreMortalityAndGrowth.R
index 5d100af..6bdfbb6 100644
--- a/tests/testthat/test-Biomass_coreMortalityAndGrowth.R
+++ b/tests/testthat/test-Biomass_coreMortalityAndGrowth.R
@@ -6,14 +6,15 @@ test_that("test growth and mortality at main simulation stage",{
                   LandR.assertions = FALSE,
                   spades.recoveryMode = FALSE)
   on.exit(options(opts))
-  Require::Require(c("SpaDES.core", "reproducible", "data.table", "raster", "LandR"))
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA, .useCache = FALSE))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
 
   endYr <- 3
   startAge <- 1L
diff --git a/tests/testthat/test-Biomass_coreNoDispersalSeeding.R b/tests/testthat/test-Biomass_coreNoDispersalSeeding.R
index faf5b97..bda89d6 100644
--- a/tests/testthat/test-Biomass_coreNoDispersalSeeding.R
+++ b/tests/testthat/test-Biomass_coreNoDispersalSeeding.R
@@ -1,11 +1,20 @@
 test_that("test no dispersal seeding algorithm",{
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath = "..",
-               outputPath = "~/output")
-  parameters <- list(.progress = list(type = "graphical", interval = 1),
-                     .globals = list(verbose = FALSE),
-                     Biomass_core = list( .saveInitialTime = NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
 
   pixelGroupMap <- raster(xmn = 50, xmx = 50 + 100 * 100,
                         ymn = 50, ymx = 50 + 100 * 100,
diff --git a/tests/testthat/test-Biomass_coreSummaryBGM.R b/tests/testthat/test-Biomass_coreSummaryBGM.R
index 5344b13..17463cf 100644
--- a/tests/testthat/test-Biomass_coreSummaryBGM.R
+++ b/tests/testthat/test-Biomass_coreSummaryBGM.R
@@ -1,12 +1,21 @@
 test_that("test summary aboveground biomass, growth, mortality. ",{
-  library(raster)
-  library(data.table)
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
+
   cohortData <- data.table(expand.grid(speciesCode = 1:3,
                                        pixelGroup = 1:5))
   cohortData[,':='(B = seq(700, by = 50, length = 15),
diff --git a/tests/testthat/test-Biomass_coreSummaryRegen.R b/tests/testthat/test-Biomass_coreSummaryRegen.R
index f3750e2..06b47a5 100644
--- a/tests/testthat/test-Biomass_coreSummaryRegen.R
+++ b/tests/testthat/test-Biomass_coreSummaryRegen.R
@@ -1,10 +1,21 @@
 test_that("test summary regeneration. ",{
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
+
   cohortData <- data.table(expand.grid(age = c(1, 10, 15),
                                        pixelGroup = 1:5))[,B:=10000]
   cohortData[age == 1, B := seq(100, by = 50, length = 5)]
diff --git a/tests/testthat/test-Biomass_coreUniversalDispersalSeeding.R b/tests/testthat/test-Biomass_coreUniversalDispersalSeeding.R
index cec6b11..bc1b431 100644
--- a/tests/testthat/test-Biomass_coreUniversalDispersalSeeding.R
+++ b/tests/testthat/test-Biomass_coreUniversalDispersalSeeding.R
@@ -1,11 +1,20 @@
 test_that("test universal dispersal seeding algorithm",{
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
 
   pixelGroupMap <- raster(xmn=50,xmx=50+100*100,
                         ymn=50,ymx=50+100*100,
diff --git a/tests/testthat/test-Biomass_coreWardDispersalSeeding.R b/tests/testthat/test-Biomass_coreWardDispersalSeeding.R
index 20e863b..76dbd87 100644
--- a/tests/testthat/test-Biomass_coreWardDispersalSeeding.R
+++ b/tests/testthat/test-Biomass_coreWardDispersalSeeding.R
@@ -1,26 +1,65 @@
 test_that("test Ward dispersal seeding algorithm",{
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath = "..",
-               outputPath = "~/output")
-  parameters <- list(.progress = list(type = "graphical", interval = 1),
-                     .globals = list(verbose = FALSE),
-                     Biomass_core = list(.saveInitialTime = NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  require("LandR")
+  require("SpaDES.core")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
 
 # 1. testing how the species seeds spread into the neighbor empty cells
 # the cohort data is set to not allow on site regeneration
-  ecoregionMap <- raster(xmn = 50, xmx = 50 + 99 * 100,
-                        ymn = 50, ymx = 50 + 99 * 100,
+
+
+  sppEquiv <- LandR::sppEquivalencies_CA[nzchar(LANDIS_test)]
+  speciesURL <- moduleInputs(module, modulePath)[objectName == "species"]$sourceURL
+  species <- getSpeciesTable(dPath = outputPath, url = speciesURL)[Area == "BSW"]
+  species <- species[sppEquiv, on = c("LandisCode" = "LANDIS_traits"), nomatch = 0]
+  setnames(species, c("Maturity", "SeedEffDist", "SeedMaxDist", "Shade", "Longevity", "MortalityCurve", "GrowthCurve"),
+           c("sexualmature", "seeddistance_eff", "seeddistance_max", "shadetolerance", "longevity", "mortalityshape", "growthcurve"))
+  spFac <- factor(species$LandR)
+  species[, speciesCode := spFac]
+
+  cohortData <- data.table(pixelGroup = 1, ecoregionGroup = factor(1L),
+                           speciesCode = spFac[7:10], age = 41, B = 8000L,
+                           mortality = 50, aNPPAct = 1079.75965551773)
+
+
+  # sim$studyArea <- randomStudyArea(seed = 1234, size = (250^2)*100)  # Jan 2021 we agreed to force user to provide a SA/SAL
+  rasterToMatch <- raster(xmn = 50 - 1.3e6, xmx = 50 + 99 * 100 - 1.3e6,
+                        ymn = 50  + 6.9e6, ymx = 50 + 99 * 100 + 6.9e6,
                         res = c(100, 100), val = 1)
-  pixelGroupMap <- setValues(ecoregionMap,2)
+  crs(rasterToMatch) <- CRS("+proj=lcc +lat_0=0 +lon_0=-95 +lat_1=49 +lat_2=77 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")
+  studyArea <- rasterToPolygons(rasterToMatch)
+
+  ecoregionFiles <- makeDummyEcoregionFiles(rasterToMatch, rasterToMatch, rasterToMatch)
+  ecoregionMap <- ecoregionFiles$ecoregionMap
+  levs <- raster::levels(ecoregionMap)[[1]]
+  colnames(levs) <- gsub("^ecoregion$", "ecoregionGroup", colnames(levs))
+  levels(ecoregionMap) <- levs
+
+  pixelGroupMap <- setValues(raster(ecoregionMap), 2)
+  nPixelsBurned <- 3000
+  pixelGroupMap[1:nPixelsBurned] <- 0
   c <- expand.grid(data.frame(a = seq(5, 99, by = 9), b = seq(5, 99, by = 9)))
   pixelindex <- (c$a - 1) * 99 + c$b #121
   pixelGroupMap[pixelindex] <- 1
   minRelativeB <- data.table(ecoregion = "eco1", X1 = 0.15,
                              X2 = 0.25, X3 = 0.5, X4 = 0.8,
-                             X5 = 0.95, ecoregionGroup = 1)
-  sufficientLight <- read.csv("~/GitHub/nrv-succession/code blitz succession/modeltesting-data/sufficientLight.csv",
-                              header = TRUE, stringsAsFactor = FALSE)
+                             X5 = 0.95, ecoregionGroup = factor(1))
+  # sufficientLight <- read.csv("~/GitHub/nrv-succession/code blitz succession/modeltesting-data/sufficientLight.csv",
+  #                             header = TRUE, stringsAsFactor = FALSE)
   lastFireYear <- "noFire"
   activePixelIndex <- 1:9801
   lastReg <- 0
@@ -28,143 +67,63 @@ test_that("test Ward dispersal seeding algorithm",{
   calibrate <- TRUE
   regenerationOutput <- data.table(seedingAlgorithm = character(), species = character(),
                                    Year = numeric(), numberOfReg = numeric())
-  species <- read.csv("~/GitHub/nrv-succession/code blitz succession/modeltesting-data/species.csv",
-                      header = TRUE, stringsAsFactors = FALSE)
-  species <- data.table(species)[, speciesCode := 1:16]
-  speciesEcoregion <- read.csv("~/GitHub/nrv-succession/code blitz succession/modeltesting-data/speciesEcoregion.csv",
-                               header = TRUE, stringsAsFactors = FALSE)
-  speciesEcoregion <- data.table(speciesEcoregion)[, ecoregionGroup := as.numeric(as.factor(ecoregion))]
-  tempsp <- setkey(species[,.(species,speciesCode)], species)
-  speciesEcoregion <- setkey(speciesEcoregion,species)[tempsp]
-  sufficientLight <- read.csv("~/GitHub/nrv-succession/code blitz succession/modeltesting-data/sufficientLight.csv",
-                              header = TRUE, stringsAsFactors = FALSE)
+
+
+  speciesEcoregion <- data.table(year = 0,
+                                 ecoregionGroup = factor(1),
+                                 speciesCode = spFac,
+                                 establishprob = 0.22,
+                                 maxANPP = 1096,
+                                 maxB = 32880L)
+  # #species <- read.csv("~/GitHub/nrv-succession/code blitz succession/modeltesting-data/species.csv",
+  # #                    header = TRUE, stringsAsFactors = FALSE)
+  # species <- data.table(species)[, speciesCode := 1:16]
+  # speciesEcoregion <- read.csv("~/GitHub/nrv-succession/code blitz succession/modeltesting-data/speciesEcoregion.csv",
+  #                              header = TRUE, stringsAsFactors = FALSE)
+  # speciesEcoregion <- data.table(speciesEcoregion)[, ecoregionGroup := as.numeric(as.factor(ecoregion))]
+  # tempsp <- setkey(species[,.(species,speciesCode)], species)
+  # speciesEcoregion <- setkey(speciesEcoregion,species)[tempsp]
+  # sufficientLight <- read.csv("~/GitHub/nrv-succession/code blitz succession/modeltesting-data/sufficientLight.csv",
+  #                             header = TRUE, stringsAsFactors = FALSE)
   seedingAlgorithm <- "wardDispersal"
-  cohortData <- data.table(pixelGroup = 1, ecoregionGroup = 1L,
-                           speciesCode = 7:10, age = 41, B = 8000L,
-                           mortality = 50, aNPPAct = 1079.75965551773)
-  objects <- list("pixelGroupMap" = pixelGroupMap,
-                  "speciesEcoregion" = speciesEcoregion,
-                  "species" = species,
-                  "successionTimestep" = successionTimestep,
-                  "calibrate" = calibrate,
-                  "seedingAlgorithm" = seedingAlgorithm,
-                  "minRelativeB" = minRelativeB,
-                  "sufficientLight" = sufficientLight,
-                  "lastFireYear" = lastFireYear,
-                  "activePixelIndex" = activePixelIndex,
-                  "lastReg" = lastReg,
-                  "regenerationOutput" = regenerationOutput,
-                  "cohortData" = cohortData,
-                  "ecoregionMap" = ecoregionMap)
-  mySim <- simInit(times = list(start = 0, end = 2),
-                   params = parameters,
-                   modules = module,
-                   objects = objects,
-                   paths = path)
-  set.seed(1)
-
-  source(file.path(modulePath(mySim), "Biomass_core", "R", "seedDispersalLANDIS.R"))
-  if (exists("Biomass_coreWardDispersalSeeding")) {
-    output <- Biomass_coreWardDispersalSeeding(mySim)
-  } else {
-    output <- mySim$.mods$Biomass_core$Biomass_coreWardDispersalSeeding(mySim)
-  }
-  output <- output$regenerationOutput$numberOfReg
-  expect_equal(output,c(283,483,288,358))
-  rm(pixelGroupMap, output,activePixelIndex,cohortData,objects,mySim)
-
-  pixelGroupMap <- setValues(ecoregionMap,0)
-  pixelGroupMap[pixelindex] <- 1
-  pixelGroupMap[1:297] <- -1
-  activePixelIndex <- 298:9801
-  pixelGroupMap[sort(pixelindex)[1:11]] <- 3
-  lastFireYear <- 0
-  cohortData <- rbind(data.table(pixelGroup = 1, ecoregionGroup = 1L,
-                                 speciesCode = 7:10, age = 41, B = 8000L,
-                                 mortality = 50, aNPPAct = 1080),
-                      data.table(pixelGroup = 3, ecoregionGroup = 1L,
-                                 speciesCode = 7:10, age = 1, B = 8000L,
-                                 mortality = 50, aNPPAct = 1090))
-  inactivePixelIndex <- 1:297
-  postFirePixel <- sort(pixelindex)[1:11]
-  objects <- list("pixelGroupMap" = pixelGroupMap,
-                  "speciesEcoregion" = speciesEcoregion,
-                  "species" = species,
-                  "successionTimestep" = successionTimestep,
-                  "calibrate" = calibrate,
-                  "seedingAlgorithm" = seedingAlgorithm,
-                  "minRelativeB" = minRelativeB,
-                  "sufficientLight" = sufficientLight,
-                  "lastFireYear" = lastFireYear,
-                  "activePixelIndex" = activePixelIndex,
-                  "lastReg" = lastReg,
-                  "regenerationOutput" = regenerationOutput,
-                  "cohortData" = cohortData,
-                  "ecoregionMap" = ecoregionMap,
-                  "inactivePixelIndex" = inactivePixelIndex,
-                  "postFirePixel" = postFirePixel)
-  mySim <- simInit(times = list(start = 0, end = 2),
-                   params = parameters,
-                   modules = module,
-                   objects = objects,
-                   paths = path)
-  set.seed(1)
-
-  source(file.path(modulePath(mySim), "Biomass_core", "R", "seedDispersalLANDIS.R"))
-  if(exists("Biomass_coreWardDispersalSeeding")){
-    output <- Biomass_coreWardDispersalSeeding(mySim)
-  } else {
-    output <- mySim$.mods$Biomass_core$Biomass_coreWardDispersalSeeding(mySim)
-  }
-  expect_equal(output$regenerationOutput$numberOfReg,
-               c(239, 439, 254, 337))
-  expect_equal(unique(output$pixelGroupMap[1:297]), -1)
-  rm(pixelGroupMap, ecoregionMap, objects, cohortData, mySim, activePixelIndex)
-
-  # 2. testing on site regeneration in ward dispersal algorithm
-  # the neiboring cells are inactive
-  ecoregionMap <- raster(xmn = 50, xmx = 50 + 99 * 100,
-                         ymn = 50, ymx = 50 + 99 * 100,
-                         res = c(100, 100), val = 2)
-  pixelGroupMap <- ecoregionMap
-  pixelGroupMap[pixelindex] <- 1
-  ecoregionMap[pixelindex] <- 1
-  activePixelIndex <- pixelindex
-  cohortData <- data.table(pixelGroup = 1, ecoregionGroup = 1L,
-                           speciesCode = 7:10, age = 41, B = 80L,
-                           mortality = 50, aNPPAct = 1079.75965551773)
-  objects <- list("pixelGroupMap" = pixelGroupMap,
-                  "speciesEcoregion" = speciesEcoregion,
-                  "species" = species,
-                  "successionTimestep" = successionTimestep,
-                  "calibrate" = calibrate,
-                  "seedingAlgorithm" = seedingAlgorithm,
-                  "minRelativeB" = minRelativeB,
-                  "sufficientLight" = sufficientLight,
-                  "lastFireYear" = lastFireYear,
-                  "activePixelIndex" = activePixelIndex,
-                  "lastReg" = lastReg,
-                  "regenerationOutput" = regenerationOutput,
-                  "cohortData" = cohortData,
-                  "ecoregionMap" = ecoregionMap)
+
+  sppEquivCol <- "LANDIS_test"
+  sppEquiv <- LandR::sppEquivalencies_CA[LANDIS_test == "poputrem"]
+  sppColorVect <- LandR::sppColors(sppEquiv, sppEquivCol, palette = "Accent")
+  initialBiomassSource <- "cohortData"
+
+  pixelCohortData <- LandR::addPixels2CohortData(cohortData, pixelGroupMap)
+  ecoregion <- LandR::makeEcoregionDT(pixelCohortData, speciesEcoregion)
+
+  objects <- mget(c("pixelGroupMap", "studyArea", "ecoregion",
+                    "rasterToMatch", "speciesEcoregion",
+                    "species", "sppEquiv", "successionTimestep", "calibrate",
+                    "seedingAlgorithm", "minRelativeB", "lastFireYear",
+                    "activePixelIndex", "lastReg", "regenerationOutput", "cohortData",
+                    "ecoregionMap"))
+
+  parameters$Biomass_core$initialBiomassSource = initialBiomassSource
+
   mySim <- simInit(times = list(start = 0, end = 2),
                    params = parameters,
                    modules = module,
                    objects = objects,
                    paths = path)
-  set.seed(1)
-  if(exists("Biomass_coreWardDispersalSeeding")){
-    output <- Biomass_coreWardDispersalSeeding(mySim)
-  } else {
-    output <- mySim$.mods$Biomass_core$Biomass_coreWardDispersalSeeding(mySim)
-  }
-  output <- output$regenerationOutput$numberOfReg
-  expect_equal(output,c(64,121,69,84))
-
+  # set.seed(1)
 
-# test ward dispersal seeding after fire disturbance
+  # source(file.path(modulePath(mySim), "Biomass_core", "R", "seedDispersalLANDIS.R"))
+  mySim$treedFirePixelTableSinceLastDisp <-
+    rbindlist(list(mySim$treedFirePixelTableSinceLastDisp,
+                   data.table(pixelIndex = 1:nPixelsBurned, pixelGroup = pixelGroupMap[][1:nPixelsBurned], burnTime = time(mySim))))
 
+  output <- spades(mySim)
 
+  expect_is(output$cohortData, "data.table")
+  news <- output$cohortData[!cohortData, on = c("speciesCode", "age")]
+  expect_true(NROW(new) > 0)
+  # Each pixel group is unique
+  expect_true(all(!duplicated(news[, fastdigest::fastdigest(speciesCode), by = "pixelGroup"]$V1)))
 
+  # remainder of seed dispersal testing is in LandR package
 
 })
diff --git a/tests/testthat/test-addNewCohorts.R b/tests/testthat/test-addNewCohorts.R
index 178e4d5..b572252 100644
--- a/tests/testthat/test-addNewCohorts.R
+++ b/tests/testthat/test-addNewCohorts.R
@@ -1,20 +1,27 @@
 test_that("test add new cohort function",{
-  # define the module and path
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
   require("raster")
   require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
 
-  module <- list("Biomass_core")
-  path <- list(modulePath = "..", # TODO: use general path
-               outputPath = "~/output") # TODO: use general path
-  parameters <- list(.progress = list(type = "graphical", interval = 1),
-                     .globals = list(verbose = FALSE),
-                     Biomass_core = list(.saveInitialTime = NA))
-
+  successionTimestep <- 10
   pixelGroupMap <- raster(xmn = 50, xmx = 50 + 3 * 100,
                           ymn = 50, ymx = 50 + 3 * 100,
                           res = c(100, 100), val = 1)
   pixelGroupMap[1] <- -1
-  pixelGroupMap[2:6] <- 2
+  pixelGroupMap[2:6] <- 0
 
   cohortData <- data.table(pixelGroup = 1, ecoregionGroup = 1L,
                            speciesCode = 7, age = 31, B = 30000L,
@@ -23,7 +30,7 @@ test_that("test add new cohort function",{
     , ':='(speciesCode = c(7, 1, 7), pixelIndex = c(7:8, 8))]
   newcohortdata2 <- rbindlist(list(cohortData, cohortData,
                                    cohortData, cohortData))[
-                                     ,':='(pixelGroup = 2,
+                                     ,':='(pixelGroup = 0,
                                            speciesCode = c(7, 7, 7, 1),
                                            pixelIndex = c(4:6, 6))]
   newcohortdata <- rbindlist(list(newcohortdata1,newcohortdata2))[
@@ -42,18 +49,25 @@ test_that("test add new cohort function",{
                    modules = module,
                    objects = objects,
                    paths = path)
-  output <- addNewCohorts(newcohortdata, cohortData, pixelGroupMap,
-                          currentTime = time(mySim), speciesEcoregion = speciesEcoregion)
+  output <- updateCohortData(newcohortdata, cohortData, pixelGroupMap,
+                          currentTime = time(mySim), speciesEcoregion = speciesEcoregion,
+                          successionTimestep = successionTimestep)
   mapOutput <- getValues(output$pixelGroupMap)
   mapOutput_compared <- c(-1, 2, 2, 4, 4, 6, 3, 5, 1)
+  mapOutput_compared <- c(-1, 0, 0, 1, 1, 2, 5, 4, 3)
   expect_equal(mapOutput,mapOutput_compared)
 
   cohortdataOutput <- output$cohortData[,.(pixelGroup, speciesCode, age, B)]
 
-  cohortdataOutput_compared <- data.table(pixelGroup = c(1, 3, 3, 4, 5, 5, 5, 6, 6),
-                                          speciesCode = c(rep(7, 4), 1, 7, 7, 1, 7),
-                                          age = c(31, 1, 31, 1, 1, 1, 31, 1, 1),
-                                          B = c(30000L, 225L, 30000L, 969L, 205L, 225L, 30000L, 886L, 969L))
+  cohortdataOutput_compared <-
+    setDT(list(pixelGroup = c(5L, 4L, 3L, 5L, 4L, 4L, 1L, 2L, 2L),
+               speciesCode = c(7, 7, 7, 7, 1, 7, 7, 7, 1),
+               age = c(31, 31, 31, 1, 1, 1, 1, 1, 1),
+               B = c(30000L, 30000L, 30000L, 225L, 206L, 225L, 969L, 969L, 886L)))
+  # cohortdataOutput_compared <- data.table(pixelGroup = c(1, 3, 3, 4, 5, 5, 5, 6, 6),
+  #                                         speciesCode = c(rep(7, 4), 1, 7, 7, 1, 7),
+  #                                         age = c(31, 1, 31, 1, 1, 1, 31, 1, 1),
+  #                                         B = c(30000L, 225L, 30000L, 969L, 205L, 225L, 30000L, 886L, 969L))
 
   expect_equal(cohortdataOutput,cohortdataOutput_compared)
 })
diff --git a/tests/testthat/test-agereclassification.R b/tests/testthat/test-agereclassification.R
index 76a9a8d..de353ce 100644
--- a/tests/testthat/test-agereclassification.R
+++ b/tests/testthat/test-agereclassification.R
@@ -1,11 +1,21 @@
 test_that("test process of age reclassification",{
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
+
   successionTimestep <- 10
   objects <- list()
   mySim <- simInit(times=list(start=0, end=1),
@@ -39,7 +49,7 @@ test_that("test process of age reclassification",{
 
   cohortData_output <- setkey(output,age)
   cohortData_output_compared <- setkey(data.table(pixelGroup = 1, ecoregionGroup = 1,
-                                speciesCode = 16, age = c(9,49), B = c(55,11),
+                                speciesCode = 16, age = c(11,49), B = c(55,11),
                                 mortality = c(1500,150), aNPPAct = c(9990,999)),age)
 
   expect_equal(cohortData_output,cohortData_output_compared)
diff --git a/tests/testthat/test-assignLightProb.R b/tests/testthat/test-assignLightProb.R
index 98e81b8..fd90e1b 100644
--- a/tests/testthat/test-assignLightProb.R
+++ b/tests/testthat/test-assignLightProb.R
@@ -1,12 +1,21 @@
 test_that("test assign light probability for a given species tolerance and site shade",{
-  library(SpaDES)
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
+
   sufficientLight <- data.frame(speciesshadetolerance=1:5,
                                 X0=seq(1,0.8,length=5),
                                 X1=seq(0.8,0.6,length=5),
diff --git a/tests/testthat/test-cacheSpinUpFunction.R b/tests/testthat/test-cacheSpinUpFunction.R
index 6c57b2c..d3364b9 100644
--- a/tests/testthat/test-cacheSpinUpFunction.R
+++ b/tests/testthat/test-cacheSpinUpFunction.R
@@ -1,10 +1,21 @@
 test_that("test cache function for spinUp",{
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
+
   useCache <- TRUE
   objects <- list("useCache"=useCache)
   mySim <- simInit(times=list(start=0, end=2),
diff --git a/tests/testthat/test-calcSiteShade.R b/tests/testthat/test-calcSiteShade.R
index 235f181..2db9e38 100644
--- a/tests/testthat/test-calcSiteShade.R
+++ b/tests/testthat/test-calcSiteShade.R
@@ -1,12 +1,21 @@
 test_that("test site shade calculation",{
-  library(SpaDES)
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
+
   minRelativeB <- data.table(ecoregion = c("eco1", "eco2"), X1 = c(0.145, 0.15),
                              X2 = c(0.217, 0.25), X3 = c(0.288, 0.5), X4 = c(0.359, 0.8),
                              X5 = c(0.430, 0.95), ecoregionGroup = 1:2)
diff --git a/tests/testthat/test-calculateANPP.R b/tests/testthat/test-calculateANPP.R
index 5525749..9e82de5 100644
--- a/tests/testthat/test-calculateANPP.R
+++ b/tests/testthat/test-calculateANPP.R
@@ -1,12 +1,21 @@
 test_that("test ANPP calculation",{
-  library(SpaDES)
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
+
   objects <- list()
   mySim <- simInit(times=list(start=0, end=1),
                    params=parameters,
diff --git a/tests/testthat/test-calculateAgeMortality.R b/tests/testthat/test-calculateAgeMortality.R
index 63212a1..bc0727f 100644
--- a/tests/testthat/test-calculateAgeMortality.R
+++ b/tests/testthat/test-calculateAgeMortality.R
@@ -1,12 +1,21 @@
 test_that("test process of age reclassification",{
-  library(SpaDES)
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
+
   spinupMortalityfraction <- 0.001
   objects <- list()
   mySim <- simInit(times=list(start=0, end=1),
diff --git a/tests/testthat/test-calculateCompetition.R b/tests/testthat/test-calculateCompetition.R
index a128acf..d7862ab 100644
--- a/tests/testthat/test-calculateCompetition.R
+++ b/tests/testthat/test-calculateCompetition.R
@@ -1,12 +1,21 @@
 test_that("test competition calculation at both spinup stage and main simulation stage",{
-  library(SpaDES)
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
+
   objects <- list()
   mySim <- simInit(times=list(start=0, end=1),
                    params=parameters,
diff --git a/tests/testthat/test-calculateGrowthMortality.R b/tests/testthat/test-calculateGrowthMortality.R
index 026e324..7dff4f9 100644
--- a/tests/testthat/test-calculateGrowthMortality.R
+++ b/tests/testthat/test-calculateGrowthMortality.R
@@ -1,12 +1,20 @@
 test_that("test Growth-related Mortality calculation",{
-  library(SpaDES)
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
   objects <- list()
   mySim <- simInit(times=list(start=0, end=1),
                    params=parameters,
diff --git a/tests/testthat/test-calculateSumB.R b/tests/testthat/test-calculateSumB.R
index 7993579..b571b2f 100644
--- a/tests/testthat/test-calculateSumB.R
+++ b/tests/testthat/test-calculateSumB.R
@@ -1,12 +1,20 @@
 test_that("test total site biomass (sumB) ",{
-  library(SpaDES)
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
   successionTimestep <- 10
   objects <- list()
   mySim <- simInit(times=list(start=0, end=1),
diff --git a/tests/testthat/test-spinUp.R b/tests/testthat/test-spinUp.R
index ed0cbe4..a1c0a73 100644
--- a/tests/testthat/test-spinUp.R
+++ b/tests/testthat/test-spinUp.R
@@ -1,11 +1,20 @@
 test_that("test growth and mortality at spinup stage",{
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
   spinupMortalityfraction <- 0.001
   calibrate <- TRUE
   successionTimestep <- 1
diff --git a/tests/testthat/test-timeScheduleofAgeReclassificationandSumB.R b/tests/testthat/test-timeScheduleofAgeReclassificationandSumB.R
index ba95045..00b09af 100644
--- a/tests/testthat/test-timeScheduleofAgeReclassificationandSumB.R
+++ b/tests/testthat/test-timeScheduleofAgeReclassificationandSumB.R
@@ -1,12 +1,20 @@
 test_that("test scheduled time for age classification and sumB at spinup stage",{
-  library(SpaDES)
-  # define the module and path
-  module <- list("Biomass_core")
-  path <- list(modulePath="..",
-               outputPath="~/output")
-  parameters <- list(.progress=list(type="graphical", interval=1),
-                     .globals=list(verbose=FALSE),
-                     Biomass_core=list( .saveInitialTime=NA))
+  opts <- options(reproducible.useGDAL = FALSE,
+                  spades.moduleCodeChecks = FALSE,
+                  reproducible.useMemoise = TRUE,
+                  spades.useRequire = FALSE,
+                  LandR.assertions = FALSE,
+                  spades.recoveryMode = FALSE)
+  on.exit(options(opts))
+  require("raster")
+  require("data.table")
+  module <- "Biomass_core"
+  modulePath <- getwd()
+  while( grepl(module, modulePath)) modulePath <- dirname(modulePath)
+  outputPath <- checkPath(file.path(tempdir(), rndstr(1)), create = TRUE)
+  path <- list(modulePath = modulePath, # TODO: use general path
+               outputPath = outputPath) # TODO: use general path
+  parameters <- list(Biomass_core = list(.saveInitialTime = NA))
   spinupMortalityfraction <- 0.001
   calibrate <- TRUE
   successionTimestep <- 1
diff --git a/tests/unitTests.R b/tests/unitTests.R
index 26a0cc1..cd2e4bf 100644
--- a/tests/unitTests.R
+++ b/tests/unitTests.R
@@ -3,7 +3,5 @@
 # please specify the package you need to run the sim function in the test files.
 
 # to test all the test files in the tests folder:
-test_dir("../Biomass_core/tests/testthat")
+try(test_dir("../Biomass_core/tests/testthat"))
 
-# Alternative, you can use test_file to test individual test file, e.g.:
-test_file("../Biomass_core/tests/testthat/test-template.R")