diff --git a/NAMESPACE b/NAMESPACE
index fa41828..972ff30 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -7,6 +7,7 @@ export(ac_to_ha)
 export(calc_crown_overlay)
 export(calc_crwidth)
 export(calc_ht_metrics)
+export(calc_landfire_stand_ht)
 export(calc_tcc_metrics)
 export(cm_to_in)
 export(create_fia_owin)
diff --git a/R/calc_ht_metrics.R b/R/calc_ht_metrics.R
index dc61dd5..55d2391 100644
--- a/R/calc_ht_metrics.R
+++ b/R/calc_ht_metrics.R
@@ -1,12 +1,24 @@
 #' Calculate stand height metrics from tree list data
 #'
-#' `calc_ht_metrics()` computes several stand height metrics for a given tree
-#' list.
+#' These functions compute several stand height metrics from tree list data.
 #'
+#' @name calc_ht_metrics
 #' @details
+#'
+#' `calc_ht_metrics()` computes several stand height metrics for a given tree
+#' list. The return value is a named list as described below.
+#'
+#' `calc_landfire_stand_ht()` computes LANDFIRE stand height. This metric is
+#' computed separately since it depends on canopy cover estimates for the
+#' sapling and overstory layers derived by overlaying modeled crowns. The input
+#' data are given as vectors of values for one or more plots (all input vectors
+#' must have the same length). The return value is a numeric vector of stand
+#' heights, with length equal to the number of elements in each of the input
+#' vectors.
+#'
 #' Stand height metrics are based on live trees (`STATUSCD == 1`), and are
-#' are assigned `0` by definition if no live trees are present. Height metrics
-#' are returned in a named list with the following elements:
+#' are assigned `0` by definition if no live trees are present.
+#' `calc_ht_metrics()`returns a named list with the following elements:
 #' * `$numTrees`: number of live trees `>= 5.0` in. (`12.7` cm) diameter
 #' * `$meanTreeHt`: mean height of trees `>= 5.0` in. (`12.7` cm) diameter
 #' * `$meanTreeHtBAW`: basal-area weighted mean height of trees `>= 5.0` in.
@@ -30,15 +42,37 @@
 #' (co-dominant), but exclude trees with `CCLCD` of `4` (intermediate) or `5`
 #' (over-topped).
 #'
+#' LANDFIRE stand height is a metric computed for FIA plots used as reference
+#' data supporting development of the Existing Vegetation Height (EVH) raster
+#' product (\url{https://www.landfire.gov/vegetation/evh}). It is generally the
+#' basal-area weighted mean height of canopy dominant/co-dominant trees (i.e.,
+#' `meanTreeHtDomBAW`). However, this metric attempts to identify plots that may
+#' be best characterized as sapling stage, in which case stand height is the
+#' mean height of saplings (`meanSapHt`). Sapling-stage plots are defined using
+#' arbitrary thresholds of canopy cover, estimated separately for the sapling
+#' layer based on microplot data, and the overstory tree layer based on subplot
+#' measurements. See `calc_tcc_metrics()` for variable definitions. A plot is
+#' considered sapling stage if `subp_overlay_mean <= 10` and
+#' `micr_overlay_mean >= 3 * subp_overlay_mean`.
+#'
 #' @param tree_list A data frame with tree records for one FIA plot.  Must have
 #' columns `DIA` (tree diameter), `HT` (tree height), `ACTUALHT` (tree actual
 #' height, `ACTUALHT < HT` indicating a broken top), `CCLCD` (FIA crown class
 #' code), `TPA_UNADJ` (trees per acre).
 #' @param digits Optional integer indicating the number of digits to keep in the
 #' return values (defaults to `1`).
-#' @return
-#' A named list of computed height metrics for the input tree list, as described
-#' in Details.
+#' @param subp_overlay_mean A numeric vector, value(s) of `subp_overlay_mean`
+#' from the output of `calc_tcc_metrics()`.
+#' @param micr_overlay_mean A numeric vector, value(s) of `micr_overlay_mean`
+#' from the output of `calc_tcc_metrics()`.
+#' @param numTrees A numeric vector, value(s) of `numTrees` from the output of
+#' `calc_ht_metrics()`.
+#' @param meanTreeHtDomBAW A numeric vector, value(s) of `meanTreeHtDomBAW` from
+#' the output of `calc_ht_metrics()`.
+#' @param meanTreeHtBAW A numeric vector, value(s) of `meanTreeHtBAW` from the
+#' output of `calc_ht_metrics()`.
+#' @param meanSapHt A numeric vector, value(s) of `meanTreeHmeanSapHttBAW` from
+#' the output of `calc_ht_metrics()`.
 #'
 #' @examples
 #' calc_ht_metrics(plantation)
@@ -141,3 +175,43 @@ calc_ht_metrics <- function(tree_list, digits = 1) {
 
     return(ht_metrics)
 }
+
+#' @name calc_ht_metrics
+#' @export
+calc_landfire_stand_ht <- function(subp_overlay_mean, micr_overlay_mean,
+                                   numTrees, meanTreeHtDomBAW, meanTreeHtBAW,
+                                   meanSapHt) {
+
+    if (!(is.numeric(subp_overlay_mean) && is.numeric(micr_overlay_mean) &&
+         is.numeric(numTrees) && is.numeric(meanTreeHtDomBAW) &&
+         is.numeric(meanTreeHtBAW) && is.numeric(meanSapHt))) {
+
+            stop("all arguments are required and must be numeric vectors",
+                 call. = FALSE)
+    }
+
+    input_vectors <- list(subp_overlay_mean, micr_overlay_mean, numTrees,
+                          meanTreeHtDomBAW, meanTreeHtBAW, meanSapHt)
+    if (length(unique(sapply(input_vectors, length))) != 1) {
+        stop("all inputs must have the same length", call. = FALSE)
+    }
+
+    standHt <- rep(NA_real_, length(subp_overlay_mean))
+
+    for (i in seq_along(standHt)) {
+        if (numTrees[i] > 0) {
+            if (subp_overlay_mean[i] <= 10 &&
+                micr_overlay_mean[i] >= 3 * subp_overlay_mean[i]) {
+                standHt[i] = meanSapHt[i]
+            } else if (meanTreeHtDomBAW[i] > meanTreeHtBAW[i]) {
+                standHt[i] = meanTreeHtDomBAW[i]
+            } else {
+                standHt[i] = meanTreeHtBAW[i]
+            }
+        } else {
+            standHt[i] = meanSapHt[i]
+        }
+    }
+
+    return(standHt)
+}
diff --git a/inst/WORDLIST b/inst/WORDLIST
index 8208186..cc552b0 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -21,6 +21,7 @@ Daryn
 Derwin
 DiTommaso
 Dockter
+EVH
 FIA
 FIADB
 FVS
diff --git a/inst/extdata/test-tree_data.csv b/inst/extdata/test-tree_data.csv
index 00f911b..4e2cfde 100644
--- a/inst/extdata/test-tree_data.csv
+++ b/inst/extdata/test-tree_data.csv
@@ -143,3 +143,24 @@ PLT_CN,SUBP,TREE,AZIMUTH,DIST,STATUSCD,SPCD,SPGRPCD,DIA,HT,ACTUALHT,CCLCD,TPA_UN
 3,4,14,291,6.2,1,6748,55,2.5,20,20,5,74.965282
 3,4,15,291,6.2,1,6748,55,1.8,18,18,5,74.965282
 3,4,16,353,6.3,1,313,41,3.1,25,25,5,74.965282
+4,3,19,33,5.4,1,935,43,1.7,20,20,5,74.965282
+4,3,18,30,5.5,1,935,43,1.5,21,21,4,74.965282
+4,3,21,155,2.5,1,375,41,1.8,23,23,4,74.965282
+4,3,22,160,5.3,1,371,30,2.7,25,25,3,74.965282
+4,3,23,178,5.9,1,761,43,2,20,20,4,74.965282
+4,3,24,220,5.8,1,375,41,1.6,25,25,4,74.965282
+4,3,25,246,5.7,1,761,43,1.9,24,24,3,74.965282
+4,3,26,260,5.8,1,761,43,1.2,20,20,5,74.965282
+4,3,27,13,6.5,1,375,41,1.5,23,23,3,74.965282
+4,3,28,50,4.7,1,761,43,1,17,17,5,74.965282
+4,3,29,61,4.2,1,371,30,1.4,21,21,4,74.965282
+4,3,30,182,6.4,1,371,30,1.7,21,21,5,74.965282
+4,3,31,200,2.9,1,375,41,1.1,16,16,5,74.965282
+4,3,32,227,5.5,1,375,41,1.1,20,20,5,74.965282
+4,3,33,265,1.1,1,761,43,1.5,17,17,3,74.965282
+4,3,34,280,5.8,1,371,30,1.2,18,18,5,74.965282
+4,3,35,323,2.3,1,12,6,1.3,10,10,5,74.965282
+4,4,17,28,22.8,1,12,6,8,30,30,2,6.018046
+4,3,16,14,4.5,1,371,30,1.8,20,20,4,74.965282
+4,3,17,27,4.1,1,761,43,2.4,22,22,3,74.965282
+4,3,20,90,5.7,1,371,30,1.6,24,24,4,74.965282
diff --git a/man/calc_ht_metrics.Rd b/man/calc_ht_metrics.Rd
index f7cf6eb..49b57b9 100644
--- a/man/calc_ht_metrics.Rd
+++ b/man/calc_ht_metrics.Rd
@@ -2,9 +2,19 @@
 % Please edit documentation in R/calc_ht_metrics.R
 \name{calc_ht_metrics}
 \alias{calc_ht_metrics}
+\alias{calc_landfire_stand_ht}
 \title{Calculate stand height metrics from tree list data}
 \usage{
 calc_ht_metrics(tree_list, digits = 1)
+
+calc_landfire_stand_ht(
+  subp_overlay_mean,
+  micr_overlay_mean,
+  numTrees,
+  meanTreeHtDomBAW,
+  meanTreeHtBAW,
+  meanSapHt
+)
 }
 \arguments{
 \item{tree_list}{A data frame with tree records for one FIA plot.  Must have
@@ -14,19 +24,43 @@ code), \code{TPA_UNADJ} (trees per acre).}
 
 \item{digits}{Optional integer indicating the number of digits to keep in the
 return values (defaults to \code{1}).}
-}
-\value{
-A named list of computed height metrics for the input tree list, as described
-in Details.
+
+\item{subp_overlay_mean}{A numeric vector, value(s) of \code{subp_overlay_mean}
+from the output of \code{calc_tcc_metrics()}.}
+
+\item{micr_overlay_mean}{A numeric vector, value(s) of \code{micr_overlay_mean}
+from the output of \code{calc_tcc_metrics()}.}
+
+\item{numTrees}{A numeric vector, value(s) of \code{numTrees} from the output of
+\code{calc_ht_metrics()}.}
+
+\item{meanTreeHtDomBAW}{A numeric vector, value(s) of \code{meanTreeHtDomBAW} from
+the output of \code{calc_ht_metrics()}.}
+
+\item{meanTreeHtBAW}{A numeric vector, value(s) of \code{meanTreeHtBAW} from the
+output of \code{calc_ht_metrics()}.}
+
+\item{meanSapHt}{A numeric vector, value(s) of \code{meanTreeHmeanSapHttBAW} from
+the output of \code{calc_ht_metrics()}.}
 }
 \description{
-\code{calc_ht_metrics()} computes several stand height metrics for a given tree
-list.
+These functions compute several stand height metrics from tree list data.
 }
 \details{
+\code{calc_ht_metrics()} computes several stand height metrics for a given tree
+list. The return value is a named list as described below.
+
+\code{calc_landfire_stand_ht()} computes LANDFIRE stand height. This metric is
+computed separately since it depends on canopy cover estimates for the
+sapling and overstory layers derived by overlaying modeled crowns. The input
+data are given as vectors of values for one or more plots (all input vectors
+must have the same length). The return value is a numeric vector of stand
+heights, with length equal to the number of elements in each of the input
+vectors.
+
 Stand height metrics are based on live trees (\code{STATUSCD == 1}), and are
-are assigned \code{0} by definition if no live trees are present. Height metrics
-are returned in a named list with the following elements:
+are assigned \code{0} by definition if no live trees are present.
+\code{calc_ht_metrics()}returns a named list with the following elements:
 \itemize{
 \item \verb{$numTrees}: number of live trees \verb{>= 5.0} in. (\code{12.7} cm) diameter
 \item \verb{$meanTreeHt}: mean height of trees \verb{>= 5.0} in. (\code{12.7} cm) diameter
@@ -51,6 +85,19 @@ For the purpose of height calculations, metrics based on
 with FIA crown class codes \code{CCLCD} of \code{1} (open grown), \code{2} (dominant) or \code{3}
 (co-dominant), but exclude trees with \code{CCLCD} of \code{4} (intermediate) or \code{5}
 (over-topped).
+
+LANDFIRE stand height is a metric computed for FIA plots used as reference
+data supporting development of the Existing Vegetation Height (EVH) raster
+product (\url{https://www.landfire.gov/vegetation/evh}). It is generally the
+basal-area weighted mean height of canopy dominant/co-dominant trees (i.e.,
+\code{meanTreeHtDomBAW}). However, this metric attempts to identify plots that may
+be best characterized as sapling stage, in which case stand height is the
+mean height of saplings (\code{meanSapHt}). Sapling-stage plots are defined using
+arbitrary thresholds of canopy cover, estimated separately for the sapling
+layer based on microplot data, and the overstory tree layer based on subplot
+measurements. See \code{calc_tcc_metrics()} for variable definitions. A plot is
+considered sapling stage if \code{subp_overlay_mean <= 10} and
+\code{micr_overlay_mean >= 3 * subp_overlay_mean}.
 }
 \examples{
 calc_ht_metrics(plantation)
diff --git a/tests/testthat/test-calc_ht_metrics.R b/tests/testthat/test-calc_ht_metrics.R
index ccf46a0..71a8ffa 100644
--- a/tests/testthat/test-calc_ht_metrics.R
+++ b/tests/testthat/test-calc_ht_metrics.R
@@ -1,4 +1,4 @@
-test_that("calc_ht_metrics works", {
+test_that("calc_ht_metrics and calc_landfire_stand_ht", {
     expected <- vector(mode = "list", length = 10)
     names(expected) <- c("numTrees", "meanTreeHt", "meanTreeHtBAW",
                          "meanTreeHtDom", "meanTreeHtDomBAW", "maxTreeHt",
@@ -16,4 +16,30 @@ test_that("calc_ht_metrics works", {
     expected$meanSapHt <- 34.5
     expected$maxSapHt <- 43
     expect_equal(calc_ht_metrics(plantation), expected, tolerance = 1e-3)
+
+    # LANDFIRE stand height
+    tcc_pred <- calc_tcc_metrics(plantation, digits = 3)
+    standHt <- calc_landfire_stand_ht(tcc_pred$subp_overlay_mean,
+                                      tcc_pred$micr_overlay_mean,
+                                      tcc_pred$numTrees,
+                                      tcc_pred$meanTreeHtDomBAW,
+                                      tcc_pred$meanTreeHtBAW,
+                                      tcc_pred$meanSapHt)
+    expect_equal(standHt, tcc_pred$meanTreeHtDomBAW, tolerance = 1e-3)
+
+    # sapling plot
+    f <- system.file("extdata/test-tree_data.csv", package="FIAstemmap")
+    tree_tbl <- load_tree_data(f)
+
+    tree_list_4 <- tree_tbl[tree_tbl$PLT_CN == "4", ]
+    expect_equal(nrow(tree_list_4), 21)
+    expect_no_error(
+        tcc_pred <- calc_tcc_metrics(tree_list_4, digits = 3))
+    standHt <- calc_landfire_stand_ht(tcc_pred$subp_overlay_mean,
+                                      tcc_pred$micr_overlay_mean,
+                                      tcc_pred$numTrees,
+                                      tcc_pred$meanTreeHtDomBAW,
+                                      tcc_pred$meanTreeHtBAW,
+                                      tcc_pred$meanSapHt)
+    expect_equal(standHt, tcc_pred$meanSapHt, tolerance = 1e-3)
 })