diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs
index 6cffe9dd35..518de3672c 100644
--- a/.git-blame-ignore-revs
+++ b/.git-blame-ignore-revs
@@ -68,3 +68,4 @@ cdf40d265cc82775607a1bf25f5f527bacc97405
 3dd489af7ebe06566e2c6a1c7ade18550f1eb4ba
 742cfa606039ab89602fde5fef46458516f56fd4
 4ad46f46de7dde753b4653c15f05326f55116b73
+75db098206b064b8b7b2a0604d3f0bf8fdb950cc
diff --git a/python/ctsm/pft_utils.py b/python/ctsm/pft_utils.py
new file mode 100644
index 0000000000..40ab8b9f23
--- /dev/null
+++ b/python/ctsm/pft_utils.py
@@ -0,0 +1,21 @@
+"""
+Constants and functions relating to PFTs
+"""
+
+MIN_PFT = 0  # bare ground
+MIN_NAT_PFT = 1  # minimum natural pft (not including bare ground)
+MAX_NAT_PFT = 14  # maximum natural pft
+MAX_PFT_GENERICCROPS = 16  # for runs with generic crops
+MAX_PFT_MANAGEDCROPS = 78  # for runs with explicit crops
+
+
+def is_valid_pft(pft_num, managed_crops):
+    """
+    Given a number, check whether it represents a valid PFT (bare ground OK)
+    """
+    if managed_crops:
+        max_allowed_pft = MAX_PFT_MANAGEDCROPS
+    else:
+        max_allowed_pft = MAX_PFT_GENERICCROPS
+
+    return MIN_PFT <= pft_num <= max_allowed_pft
diff --git a/python/ctsm/site_and_regional/plumber2_shared.py b/python/ctsm/site_and_regional/plumber2_shared.py
new file mode 100644
index 0000000000..d4ab9d00b3
--- /dev/null
+++ b/python/ctsm/site_and_regional/plumber2_shared.py
@@ -0,0 +1,21 @@
+"""
+Things shared between plumber2 scripts
+"""
+
+import os
+import pandas as pd
+from ctsm.path_utils import path_to_ctsm_root
+
+PLUMBER2_SITES_CSV = os.path.join(
+    path_to_ctsm_root(),
+    "tools",
+    "site_and_regional",
+    "PLUMBER2_sites.csv",
+)
+
+
+def read_plumber2_sites_csv(file=PLUMBER2_SITES_CSV):
+    """
+    Read PLUMBER2_sites.csv using pandas
+    """
+    return pd.read_csv(file, skiprows=5)
diff --git a/python/ctsm/site_and_regional/plumber2_surf_wrapper.py b/python/ctsm/site_and_regional/plumber2_surf_wrapper.py
index 022914d17e..cedc6b25e0 100755
--- a/python/ctsm/site_and_regional/plumber2_surf_wrapper.py
+++ b/python/ctsm/site_and_regional/plumber2_surf_wrapper.py
@@ -22,16 +22,18 @@
 
 import argparse
 import logging
-import os
-import subprocess
+import sys
 import tqdm
 
-import pandas as pd
+# pylint:disable=wrong-import-position
+from ctsm.site_and_regional.plumber2_shared import PLUMBER2_SITES_CSV, read_plumber2_sites_csv
+from ctsm import subset_data
+from ctsm.pft_utils import MAX_PFT_MANAGEDCROPS, is_valid_pft
 
 
-def get_parser():
+def get_args():
     """
-    Get parser object for this script.
+    Get arguments for this script.
     """
     parser = argparse.ArgumentParser(
         description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
@@ -45,39 +47,44 @@ def get_parser():
         help="Verbose mode will print more information. ",
         action="store_true",
         dest="verbose",
-        default=False,
     )
 
     parser.add_argument(
-        "--16pft",
-        help="Create and/or modify 16-PFT surface datasets (e.g. for a FATES run) ",
+        "--crop",
+        help=f"Create and/or modify {MAX_PFT_MANAGEDCROPS}-PFT "
+        "surface datasets (e.g. for a non-FATES run)",
         action="store_true",
-        dest="pft_16",
-        default=True,
+        dest="use_managed_crops",
     )
 
-    return parser
+    parser.add_argument(
+        "--overwrite",
+        help="Overwrite any existing files",
+        action="store_true",
+    )
+
+    parser.add_argument(
+        "--plumber2-sites-csv",
+        help=f"Comma-separated value (CSV) file with Plumber2 sites. Default: {PLUMBER2_SITES_CSV}",
+        default=PLUMBER2_SITES_CSV,
+    )
+
+    return parser.parse_args()
 
 
 def execute(command):
     """
-    Function for running a command on shell.
+    Runs subset_data with given arguments.
     Args:
-        command (str):
-            command that we want to run.
+        command (list):
+            list of args for command that we want to run.
     Raises:
-        Error with the return code from shell.
+        Whatever error subset_data gives, if any.
     """
     print("\n", " >>  ", *command, "\n")
 
-    try:
-        subprocess.check_call(command, stdout=open(os.devnull, "w"), stderr=subprocess.STDOUT)
-
-    except subprocess.CalledProcessError as err:
-        # raise RuntimeError("command '{}' return with error
-        # (code {}): {}".format(e.cmd, e.returncode, e.output))
-        # print (e.ouput)
-        print(err)
+    sys.argv = command
+    subset_data.main()
 
 
 def main():
@@ -85,97 +92,103 @@ def main():
     Read plumber2_sites from csv, iterate through sites, and add dominant PFT
     """
 
-    args = get_parser().parse_args()
+    args = get_args()
 
     if args.verbose:
         logging.basicConfig(level=logging.DEBUG)
 
-    plumber2_sites = pd.read_csv("PLUMBER2_sites.csv", skiprows=4)
+    plumber2_sites = read_plumber2_sites_csv(args.plumber2_sites_csv)
 
     for _, row in tqdm.tqdm(plumber2_sites.iterrows()):
         lat = row["Lat"]
         lon = row["Lon"]
         site = row["Site"]
+
+        clmsite = "1x1_PLUMBER2_" + site
+        print("Now processing site :", site)
+
+        # Set up part of subset_data command that is shared among all options
+        subset_command = [
+            "./subset_data",
+            "point",
+            "--lat",
+            str(lat),
+            "--lon",
+            str(lon),
+            "--site",
+            clmsite,
+            "--create-surface",
+            "--uniform-snowpack",
+            "--cap-saturation",
+            "--lon-type",
+            "180",
+        ]
+
+        # Read info for first PFT
         pft1 = row["pft1"]
+        if not is_valid_pft(pft1, args.use_managed_crops):
+            raise RuntimeError(f"pft1 must be a valid PFT; got {pft1}")
         pctpft1 = row["pft1-%"]
         cth1 = row["pft1-cth"]
         cbh1 = row["pft1-cbh"]
-        pft2 = row["pft2"]
-        pctpft2 = row["pft2-%"]
-        cth2 = row["pft2-cth"]
-        cbh2 = row["pft2-cbh"]
-        # overwrite missing values from .csv file
-        if pft1 == -999:
-            pft1 = 0
-            pctpft1 = 0
-            cth1 = 0
-            cbh1 = 0
-        if pft2 == -999:
-            pft2 = 0
-            pctpft2 = 0
-            cth2 = 0
-            cbh2 = 0
-        clmsite = "1x1_PLUMBER2_" + site
-        print("Now processing site :", site)
 
-        if args.pft_16:
-            # use surface dataset with 16 pfts, but overwrite to 100% 1 dominant PFT
-            # don't set crop flag
-            # set dominant pft
-            subset_command = [
-                "./subset_data",
-                "point",
-                "--lat",
-                str(lat),
-                "--lon",
-                str(lon),
-                "--site",
-                clmsite,
+        # Read info for second PFT, if a valid one is given in the .csv file
+        pft2 = row["pft2"]
+        if is_valid_pft(pft2, args.use_managed_crops):
+            pctpft2 = row["pft2-%"]
+            cth2 = row["pft2-cth"]
+            cbh2 = row["pft2-cbh"]
+
+        # Set dominant PFT(s)
+        if is_valid_pft(pft2, args.use_managed_crops):
+            subset_command += [
                 "--dompft",
                 str(pft1),
                 str(pft2),
                 "--pctpft",
                 str(pctpft1),
                 str(pctpft2),
-                "--cth",
-                str(cth1),
-                str(cth2),
-                "--cbh",
-                str(cbh1),
-                str(cbh2),
-                "--create-surface",
-                "--uniform-snowpack",
-                "--cap-saturation",
-                "--verbose",
-                "--overwrite",
             ]
         else:
-            # use surface dataset with 78 pfts, and overwrite to 100% 1 dominant PFT
-            # NOTE: FATES will currently not run with a 78-PFT surface dataset
-            # set crop flag
-            # set dominant pft
-            subset_command = [
-                "./subset_data",
-                "point",
-                "--lat",
-                str(lat),
-                "--lon",
-                str(lon),
-                "--site",
-                clmsite,
-                "--crop",
+            subset_command += [
                 "--dompft",
                 str(pft1),
-                str(pft2),
                 "--pctpft",
                 str(pctpft1),
-                str(pctpft2),
-                "--create-surface",
-                "--uniform-snowpack",
-                "--cap-saturation",
-                "--verbose",
-                "--overwrite",
             ]
+
+        if not args.use_managed_crops:
+            # use surface dataset with 78 pfts, but overwrite to 100% 1 dominant PFT
+            # don't set crop flag
+            # set canopy top and bottom heights
+            if is_valid_pft(pft2, args.use_managed_crops):
+                subset_command += [
+                    "--cth",
+                    str(cth1),
+                    str(cth2),
+                    "--cbh",
+                    str(cbh1),
+                    str(cbh2),
+                ]
+            else:
+                subset_command += [
+                    "--cth",
+                    str(cth1),
+                    "--cbh",
+                    str(cbh1),
+                ]
+        else:
+            # use surface dataset with 78 pfts, and overwrite to 100% 1 dominant PFT
+            # NOTE: FATES will currently not run with a 78-PFT surface dataset
+            # set crop flag
+            subset_command += ["--crop"]
+            # don't set canopy top and bottom heights
+
+        if args.verbose:
+            subset_command += ["--verbose"]
+        if args.overwrite:
+            subset_command += ["--overwrite"]
+
         execute(subset_command)
 
 
diff --git a/python/ctsm/site_and_regional/plumber2_usermods.py b/python/ctsm/site_and_regional/plumber2_usermods.py
index 7b7f294a24..7c8f37b1b5 100644
--- a/python/ctsm/site_and_regional/plumber2_usermods.py
+++ b/python/ctsm/site_and_regional/plumber2_usermods.py
@@ -13,7 +13,8 @@
 import os
 import tqdm
 
-import pandas as pd
+# pylint:disable=wrong-import-position
+from ctsm.site_and_regional.plumber2_shared import read_plumber2_sites_csv
 
 
 # Big ugly function to create usermod_dirs for each site
@@ -155,7 +156,7 @@ def main():
     """
 
     # For now we can just run the 'main' program as a loop
-    plumber2_sites = pd.read_csv("PLUMBER2_sites.csv", skiprows=4)
+    plumber2_sites = read_plumber2_sites_csv()
 
     for _, row in tqdm.tqdm(plumber2_sites.iterrows()):
         lat = row["Lat"]
diff --git a/python/ctsm/site_and_regional/single_point_case.py b/python/ctsm/site_and_regional/single_point_case.py
index e122b2c251..d71d014f36 100644
--- a/python/ctsm/site_and_regional/single_point_case.py
+++ b/python/ctsm/site_and_regional/single_point_case.py
@@ -16,13 +16,10 @@
 from ctsm.site_and_regional.base_case import BaseCase, USRDAT_DIR, DatmFiles
 from ctsm.utils import add_tag_to_filename, ensure_iterable
 from ctsm.longitude import detect_lon_type
+from ctsm.pft_utils import MAX_NAT_PFT, MAX_PFT_GENERICCROPS, MAX_PFT_MANAGEDCROPS
 
 logger = logging.getLogger(__name__)
 
-NAT_PFT = 15  # natural pfts
-NUM_PFT = 17  # for runs with generic crops
-MAX_PFT = 78  # for runs with explicit crops
-
 
 class SinglePointCase(BaseCase):
     """
@@ -190,20 +187,21 @@ def check_dom_pft(self):
           same range.
           e.g. If users specified multiple dom_pft, they should be
           either in :
-            - 0 - NAT_PFT-1 range
+            - 0 - MAX_NAT_PFT range
             or
-            - NAT_PFT - MAX_PFT range
+            - MAX_NAT_PFT+1 - MAX_PFT_MANAGEDCROPS range
             - give an error: mixed land units not possible
 
         -------------
         Raises:
             Error (ArgumentTypeError):
-                If any dom_pft is bigger than MAX_PFT.
+                If any dom_pft is bigger than MAX_PFT_MANAGEDCROPS.
             Error (ArgumentTypeError):
                 If any dom_pft is less than 1.
             Error (ArgumentTypeError):
                 If mixed land units are chosen.
-                dom_pft values are both in range of (0 - NAT_PFT-1) and (NAT_PFT - MAX_PFT).
+                dom_pft values are both in range of
+                (0 - MAX_NAT_PFT) and (MAX_NAT_PFT+1 - MAX_PFT_MANAGEDCROPS).
 
 
         """
@@ -217,27 +215,29 @@ def check_dom_pft(self):
             min_dom_pft = min(self.dom_pft)
             max_dom_pft = max(self.dom_pft)
 
-            # -- check dom_pft values should be between 0-MAX_PFT
-            if min_dom_pft < 0 or max_dom_pft > MAX_PFT:
-                err_msg = "values for --dompft should be between 1 and 78."
+            # -- check dom_pft values should be between 0-MAX_PFT_MANAGEDCROPS
+            if min_dom_pft < 0 or max_dom_pft > MAX_PFT_MANAGEDCROPS:
+                err_msg = f"values for --dompft should be between 1 and {MAX_PFT_MANAGEDCROPS}."
                 raise argparse.ArgumentTypeError(err_msg)
 
             # -- check dom_pft vs num_pft
             if max_dom_pft > self.num_pft:
-                err_msg = "Please use --crop flag when --dompft is above 16."
+                err_msg = f"Please use --crop flag when --dompft is above {MAX_PFT_GENERICCROPS}."
                 raise argparse.ArgumentTypeError(err_msg)
 
             # -- check dom_pft vs MAX_pft
-            if self.num_pft - 1 < max_dom_pft < NUM_PFT:
+            if self.num_pft - 1 < max_dom_pft <= MAX_PFT_GENERICCROPS:
                 logger.info(
-                    "WARNING, you trying to run with generic crops (16 PFT surface dataset)"
+                    "WARNING, you are trying to run with generic crops (%s PFT surface dataset)",
+                    MAX_PFT_GENERICCROPS,
                 )
 
             # -- check if all dom_pft are in the same range:
-            if min_dom_pft < NAT_PFT <= max_dom_pft:
+            if min_dom_pft <= MAX_NAT_PFT < max_dom_pft:
                 err_msg = (
                     "You are subsetting using mixed land units that have both "
-                    "natural pfts and crop cfts. Check your surface dataset. "
+                    "natural pfts and crop cfts. Check your surface dataset.\n"
+                    f"{min_dom_pft} <= {MAX_NAT_PFT} < {max_dom_pft}\n"
                 )
                 raise argparse.ArgumentTypeError(err_msg)
 
@@ -425,7 +425,7 @@ def modify_surfdata_atpoint(self, f_orig):
         if self.dom_pft is not None:
             max_dom_pft = max(self.dom_pft)
             # -- First initialize everything:
-            if max_dom_pft < NAT_PFT:
+            if max_dom_pft <= MAX_NAT_PFT:
                 f_mod["PCT_NAT_PFT"][:, :, :] = 0
             else:
                 f_mod["PCT_CFT"][:, :, :] = 0
@@ -444,10 +444,10 @@ def modify_surfdata_atpoint(self, f_orig):
                 if cth is not None:
                     f_mod["MONTHLY_HEIGHT_TOP"][:, :, :, dom_pft] = cth
                     f_mod["MONTHLY_HEIGHT_BOT"][:, :, :, dom_pft] = cbh
-                if dom_pft < NAT_PFT:
+                if dom_pft <= MAX_NAT_PFT:
                     f_mod["PCT_NAT_PFT"][:, :, dom_pft] = pct_pft
                 else:
-                    dom_pft = dom_pft - NAT_PFT
+                    dom_pft = dom_pft - (MAX_NAT_PFT + 1)
                     f_mod["PCT_CFT"][:, :, dom_pft] = pct_pft
 
         # -------------------------------
@@ -465,7 +465,7 @@ def modify_surfdata_atpoint(self, f_orig):
 
             if self.dom_pft is not None:
                 max_dom_pft = max(self.dom_pft)
-                if max_dom_pft < NAT_PFT:
+                if max_dom_pft <= MAX_NAT_PFT:
                     f_mod["PCT_NATVEG"][:, :] = 100
                     f_mod["PCT_CROP"][:, :] = 0
                 else:
diff --git a/python/ctsm/subset_data.py b/python/ctsm/subset_data.py
index 820391ff8b..de4e51db9b 100644
--- a/python/ctsm/subset_data.py
+++ b/python/ctsm/subset_data.py
@@ -70,6 +70,7 @@
 from ctsm.utils import abort
 from ctsm.config_utils import check_lon1_lt_lon2
 from ctsm.longitude import Longitude, detect_lon_type
+from ctsm.pft_utils import MAX_PFT_GENERICCROPS, MAX_PFT_MANAGEDCROPS
 
 # -- import ctsm logging flags
 from ctsm.ctsm_logging import (
@@ -597,9 +598,9 @@ def determine_num_pft(crop):
         num_pft (int) : number of pfts for surface dataset
     """
     if crop:
-        num_pft = "78"
+        num_pft = str(MAX_PFT_MANAGEDCROPS)
     else:
-        num_pft = "16"
+        num_pft = str(MAX_PFT_GENERICCROPS)
     logger.debug("crop_flag = %s => num_pft = %s", str(crop), num_pft)
     return num_pft
 
diff --git a/python/ctsm/test/test_sys_plumber2_surf_wrapper.py b/python/ctsm/test/test_sys_plumber2_surf_wrapper.py
new file mode 100755
index 0000000000..12ca561150
--- /dev/null
+++ b/python/ctsm/test/test_sys_plumber2_surf_wrapper.py
@@ -0,0 +1,147 @@
+#!/usr/bin/env python3
+
+"""System tests for plumber2_surf_wrapper"""
+
+import os
+import unittest
+import tempfile
+import shutil
+import sys
+
+from ctsm import unit_testing
+from ctsm.site_and_regional.plumber2_surf_wrapper import main
+from ctsm.site_and_regional.plumber2_shared import read_plumber2_sites_csv
+from ctsm.path_utils import path_to_ctsm_root
+
+# Allow test names that pylint doesn't like; otherwise hard to make them
+# readable
+# pylint: disable=invalid-name
+
+
+class TestSysPlumber2SurfWrapper(unittest.TestCase):
+    """
+    System tests for plumber2_surf_wrapper
+    """
+
+    def setUp(self):
+        """
+        Make tempdir for use by these tests.
+        """
+        self._previous_dir = os.getcwd()
+        self._tempdir = tempfile.mkdtemp()
+        os.chdir(self._tempdir)  # cd to tempdir
+
+        # Path to script
+        self.tool_path = os.path.join(
+            path_to_ctsm_root(),
+            "tools",
+            "site_and_regional",
+            "plumber2_surf_wrapper",
+        )
+
+        # Path to test inputs directory
+        self.test_inputs = os.path.join(
+            os.path.dirname(__file__), "testinputs", "plumber2_surf_wrapper"
+        )
+
+    def tearDown(self):
+        """
+        Remove temporary directory
+        """
+        os.chdir(self._previous_dir)
+        shutil.rmtree(self._tempdir, ignore_errors=True)
+
+    def test_plumber2_surf_wrapper(self):
+        """
+        Run the entire tool with default settings.
+        CAN ONLY RUN ON SYSTEMS WITH INPUTDATA
+        """
+
+        sys.argv = [self.tool_path]
+        main()
+
+        # How many files do we expect?
+        plumber2_csv = read_plumber2_sites_csv()
+        n_files_expected = len(plumber2_csv)
+
+        # How many files did we get?
+        file_list = os.listdir("subset_data_single_point")
+        n_files = len(file_list)
+
+        # Check
+        self.assertEqual(n_files_expected, n_files)
+
+    def test_plumber2_surf_wrapper_78pft(self):
+        """
+        Run the entire tool with --crop.
+        CAN ONLY RUN ON SYSTEMS WITH INPUTDATA
+        """
+
+        sys.argv = [self.tool_path, "--crop"]
+        main()
+
+        # How many files do we expect?
+        plumber2_csv = read_plumber2_sites_csv()
+        n_files_expected = len(plumber2_csv)
+
+        # How many files did we get?
+        file_list = os.listdir("subset_data_single_point")
+        n_files = len(file_list)
+
+        # Check
+        self.assertEqual(n_files_expected, n_files)
+
+    def test_plumber2_surf_wrapper_invalid_pft(self):
+        """
+        plumber2_surf_wrapper should error if invalid PFT is given
+        """
+
+        sys.argv = [
+            self.tool_path,
+            "--plumber2-sites-csv",
+            os.path.join(self.test_inputs, "PLUMBER2_sites_invalid_pft.csv"),
+        ]
+        with self.assertRaisesRegex(RuntimeError, "must be a valid PFT"):
+            main()
+
+    def test_plumber2_surf_wrapper_existing_no_overwrite_fails(self):
+        """
+        plumber2_surf_wrapper should fail if file exists but --overwrite isn't given
+        """
+
+        sys_argv_shared = [
+            self.tool_path,
+            "--plumber2-sites-csv",
+            os.path.join(self.test_inputs, "PLUMBER2_site_valid.csv"),
+        ]
+
+        # Run twice, expecting second to fail
+        sys.argv = sys_argv_shared
+        main()
+        sys.argv = sys_argv_shared
+        with self.assertRaisesRegex(SystemExit, "exists"):
+            main()
+
+    def test_plumber2_surf_wrapper_existing_overwrite_passes(self):
+        """
+        plumber2_surf_wrapper should pass if file exists and --overwrite is given
+        """
+
+        sys_argv_shared = [
+            self.tool_path,
+            "--plumber2-sites-csv",
+            os.path.join(self.test_inputs, "PLUMBER2_site_valid.csv"),
+        ]
+
+        # Run once to generate the files
+        sys.argv = sys_argv_shared
+        main()
+
+        # Run again with --overwrite, expecting pass
+        sys.argv = sys_argv_shared + ["--overwrite"]
+        main()
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_unit_plumber2_surf_wrapper.py b/python/ctsm/test/test_unit_plumber2_surf_wrapper.py
index 66f5578caa..4b84752edb 100755
--- a/python/ctsm/test/test_unit_plumber2_surf_wrapper.py
+++ b/python/ctsm/test/test_unit_plumber2_surf_wrapper.py
@@ -16,7 +16,7 @@
 
 # pylint: disable=wrong-import-position
 from ctsm import unit_testing
-from ctsm.site_and_regional.plumber2_surf_wrapper import get_parser
+from ctsm.site_and_regional.plumber2_surf_wrapper import get_args
 
 # pylint: disable=invalid-name
 
@@ -26,12 +26,60 @@ class TestPlumber2SurfWrapper(unittest.TestCase):
     Basic class for testing plumber2_surf_wrapper.py.
     """
 
-    def test_parser(self):
+    def setUp(self):
+        sys.argv = ["subset_data"]  # Could actually be anything
+
+    def test_parser_default_csv_exists(self):
+        """
+        Test that default PLUMBER2 sites CSV file exists
+        """
+
+        args = get_args()
+        self.assertTrue(os.path.exists(args.plumber2_sites_csv))
+
+    def test_parser_custom_csv(self):
+        """
+        Test that script accepts custom CSV file path
+        """
+
+        custom_path = "path/to/custom.csv"
+        sys.argv += ["--plumber2-sites-csv", custom_path]
+        args = get_args()
+        self.assertEqual(args.plumber2_sites_csv, custom_path)
+
+    def test_parser_verbose_false_default(self):
+        """
+        Test that script is not verbose by default
+        """
+
+        args = get_args()
+        self.assertFalse(args.verbose)
+
+    def test_parser_verbose_true(self):
+        """
+        Test that --verbose sets verbose to True
+        """
+
+        sys.argv += ["--verbose"]
+        args = get_args()
+        self.assertTrue(args.verbose)
+
+    def test_parser_78pft_false_default(self):
+        """
+        Test that script does not use 78pft mode by default
+        """
+
+        args = get_args()
+        self.assertFalse(args.use_managed_crops)
+
+    def test_parser_78pft_true(self):
         """
-        Test that parser has same defaults as expected
+        Test that --crop sets use_managed_crops to True
         """
 
-        self.assertEqual(get_parser().argument_default, None, "Parser not working as expected")
+        sys.argv += ["--crop"]
+        args = get_args()
+        self.assertTrue(args.use_managed_crops)
 
 
 if __name__ == "__main__":
diff --git a/python/ctsm/test/test_unit_singlept_data.py b/python/ctsm/test/test_unit_singlept_data.py
index 644af82588..bf29ced331 100755
--- a/python/ctsm/test/test_unit_singlept_data.py
+++ b/python/ctsm/test/test_unit_singlept_data.py
@@ -18,6 +18,7 @@
 # pylint: disable=wrong-import-position
 from ctsm import unit_testing
 from ctsm.site_and_regional.single_point_case import SinglePointCase
+from ctsm.pft_utils import MAX_PFT_GENERICCROPS, MAX_PFT_MANAGEDCROPS
 
 # pylint: disable=invalid-name
 
@@ -38,7 +39,7 @@ class TestSinglePointCase(unittest.TestCase):
     dom_pft = [8]
     evenly_split_cropland = False
     pct_pft = None
-    num_pft = 16
+    num_pft = MAX_PFT_GENERICCROPS
     cth = [0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9]
     cbh = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
     include_nonveg = False
@@ -131,7 +132,7 @@ def test_check_dom_pft_too_big(self):
             out_dir=self.out_dir,
             overwrite=self.overwrite,
         )
-        single_point.dom_pft = [16, 36, 79]
+        single_point.dom_pft = [MAX_PFT_GENERICCROPS, 36, 79]
         with self.assertRaisesRegex(argparse.ArgumentTypeError, "values for --dompft should*"):
             single_point.check_dom_pft()
 
@@ -161,7 +162,7 @@ def test_check_dom_pft_too_small(self):
             out_dir=self.out_dir,
             overwrite=self.overwrite,
         )
-        single_point.dom_pft = [16, 36, -1]
+        single_point.dom_pft = [MAX_PFT_GENERICCROPS, 36, -1]
         with self.assertRaisesRegex(argparse.ArgumentTypeError, "values for --dompft should*"):
             single_point.check_dom_pft()
 
@@ -192,7 +193,7 @@ def test_check_dom_pft_numpft(self):
             overwrite=self.overwrite,
         )
         single_point.dom_pft = [15, 53]
-        single_point.num_pft = 16
+        single_point.num_pft = MAX_PFT_GENERICCROPS
         with self.assertRaisesRegex(argparse.ArgumentTypeError, "Please use --crop*"):
             single_point.check_dom_pft()
 
@@ -223,7 +224,7 @@ def test_check_dom_pft_mixed_range(self):
             overwrite=self.overwrite,
         )
         single_point.dom_pft = [1, 5, 15]
-        single_point.num_pft = 78
+        single_point.num_pft = MAX_PFT_MANAGEDCROPS
         with self.assertRaisesRegex(
             argparse.ArgumentTypeError, "You are subsetting using mixed land*"
         ):
diff --git a/python/ctsm/test/test_unit_singlept_data_surfdata.py b/python/ctsm/test/test_unit_singlept_data_surfdata.py
index 2106799a4b..d163c29e4f 100755
--- a/python/ctsm/test/test_unit_singlept_data_surfdata.py
+++ b/python/ctsm/test/test_unit_singlept_data_surfdata.py
@@ -23,6 +23,7 @@
 # pylint: disable=wrong-import-position
 from ctsm import unit_testing
 from ctsm.site_and_regional.single_point_case import SinglePointCase
+from ctsm.pft_utils import MAX_PFT_GENERICCROPS, MAX_PFT_MANAGEDCROPS
 
 # pylint: disable=invalid-name
 # pylint: disable=too-many-lines
@@ -46,7 +47,7 @@ class TestSinglePointCaseSurfaceNoCrop(unittest.TestCase):
     dom_pft = [8]
     evenly_split_cropland = False
     pct_pft = None
-    num_pft = 16
+    num_pft = MAX_PFT_GENERICCROPS
     cth = 0.9
     cbh = 0.1
     include_nonveg = False
@@ -667,7 +668,7 @@ class TestSinglePointCaseSurfaceCrop(unittest.TestCase):
     dom_pft = [17]
     evenly_split_cropland = False
     pct_pft = None
-    num_pft = 78
+    num_pft = MAX_PFT_MANAGEDCROPS
     cth = 0.9
     cbh = 0.1
     include_nonveg = False
diff --git a/python/ctsm/test/testinputs/plumber2_surf_wrapper/PLUMBER2_site_valid.csv b/python/ctsm/test/testinputs/plumber2_surf_wrapper/PLUMBER2_site_valid.csv
new file mode 100644
index 0000000000..2c1580bc03
--- /dev/null
+++ b/python/ctsm/test/testinputs/plumber2_surf_wrapper/PLUMBER2_site_valid.csv
@@ -0,0 +1,7 @@
+#pftX-cth and pftX-cbh are the site=specific canopy top and bottom heights
+#start_year and end_year will be used to define DATM_YR_ALIGH, DATM_YR_START and DATM_YR_END, and STOP_N in units of nyears.
+#RUN_STARTDATE and START_TOD are specified because we are starting at GMT corresponding to local midnight.
+#ATM_NCPL is specified so that the time step of the model matches the time interval specified by the atm forcing data.
+#longitudes must be in the range [-180,180]
+,Site,Lat,Lon,pft1,pft1-%,pft1-cth,pft1-cbh,pft2,pft2-%,pft2-cth,pft2-cbh,start_year,end_year,RUN_STARTDATE,START_TOD,ATM_NCPL
+27,BE-Lon,50.551590, 4.746130,15,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2005,2014,2004-12-31,82800,48 
diff --git a/python/ctsm/test/testinputs/plumber2_surf_wrapper/PLUMBER2_sites_invalid_pft.csv b/python/ctsm/test/testinputs/plumber2_surf_wrapper/PLUMBER2_sites_invalid_pft.csv
new file mode 100644
index 0000000000..e8f0eb8fbb
--- /dev/null
+++ b/python/ctsm/test/testinputs/plumber2_surf_wrapper/PLUMBER2_sites_invalid_pft.csv
@@ -0,0 +1,8 @@
+#pftX-cth and pftX-cbh are the site=specific canopy top and bottom heights
+#start_year and end_year will be used to define DATM_YR_ALIGH, DATM_YR_START and DATM_YR_END, and STOP_N in units of nyears.
+#RUN_STARTDATE and START_TOD are specified because we are starting at GMT corresponding to local midnight.
+#ATM_NCPL is specified so that the time step of the model matches the time interval specified by the atm forcing data.
+#longitudes must be in the range [-180,180]
+,Site,Lat,Lon,pft1,pft1-%,pft1-cth,pft1-cbh,pft2,pft2-%,pft2-cth,pft2-cbh,start_year,end_year,RUN_STARTDATE,START_TOD,ATM_NCPL
+26,Invalid-Pft,51.309166, 4.520560,-1,19.22,21.00,10.50,7,80.78,21.00,12.08,2004,2014,2003-12-31,82800,48 
+27,BE-Lon,50.551590, 4.746130,15,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2005,2014,2004-12-31,82800,48 
diff --git a/python/ctsm/toolchain/gen_mksurfdata_namelist.py b/python/ctsm/toolchain/gen_mksurfdata_namelist.py
index 31fcbfe8ff..3a405bf5fa 100755
--- a/python/ctsm/toolchain/gen_mksurfdata_namelist.py
+++ b/python/ctsm/toolchain/gen_mksurfdata_namelist.py
@@ -15,6 +15,7 @@
 
 from ctsm.path_utils import path_to_ctsm_root, path_to_cime
 from ctsm.ctsm_logging import setup_logging_pre_config, add_logging_args, process_logging_args
+from ctsm.pft_utils import MAX_PFT_GENERICCROPS, MAX_PFT_MANAGEDCROPS
 
 logger = logging.getLogger(__name__)
 
@@ -306,9 +307,9 @@ def main():
 
     # Determine num_pft
     if nocrop_flag:
-        num_pft = "16"
+        num_pft = str(MAX_PFT_GENERICCROPS)
     else:
-        num_pft = "78"
+        num_pft = str(MAX_PFT_MANAGEDCROPS)
     logger.info("num_pft is %s", num_pft)
 
     # Write out if surface dataset will be created
diff --git a/tools/site_and_regional/PLUMBER2_sites.csv b/tools/site_and_regional/PLUMBER2_sites.csv
index f252fa1d61..1097568051 100644
--- a/tools/site_and_regional/PLUMBER2_sites.csv
+++ b/tools/site_and_regional/PLUMBER2_sites.csv
@@ -2,6 +2,7 @@
 #start_year and end_year will be used to define DATM_YR_ALIGH, DATM_YR_START and DATM_YR_END, and STOP_N in units of nyears.
 #RUN_STARTDATE and START_TOD are specified because we are starting at GMT corresponding to local midnight.
 #ATM_NCPL is specified so that the time step of the model matches the time interval specified by the atm forcing data.
+#longitudes must be in the range [-180,180]
 ,Site,Lat,Lon,pft1,pft1-%,pft1-cth,pft1-cbh,pft2,pft2-%,pft2-cth,pft2-cbh,start_year,end_year,RUN_STARTDATE,START_TOD,ATM_NCPL
 1,AR-SLu,-33.464802,-66.459808,5,50.00, 4.50, 0.13,7,50.00, 4.50, 2.59,2010,2010,2010-01-01,10800,48 
 2,AT-Neu,47.116669,11.317500,13,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2002,2012,2001-12-31,82800,48 
@@ -73,7 +74,7 @@
 68,DK-Sor,55.485870,11.644640,7,100.00,25.00,14.37,-999,-999.00,-999.00,-999.00,1997,2014,1996-12-31,82800,48 
 69,DK-ZaH,74.473282,-20.550293,12,100.00, 0.47, 0.01,-999,-999.00,-999.00,-999.00,2000,2013,2000-01-01,0,48 
 70,ES-ES1,39.345970,-0.318817,1,100.00, 7.50, 3.75,-999,-999.00,-999.00,-999.00,1999,2006,1998-12-31,82800,48 
-71,ES-ES2,39.275558,-0.315277,-999,-999.00,-999.00,-999.00,16,100.00, 0.50, 0.01,2005,2006,2004-12-31,82800,48 
+71,ES-ES2,39.275558,-0.315277,16,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2005,2006,2004-12-31,82800,48 
 72,ES-LgS,37.097935,-2.965820,10,30.00, 0.20, 0.04,13,70.00, 0.50, 0.01,2007,2007,2006-12-31,82800,48 
 73,ES-LMa,39.941502,-5.773346,7,30.00, 8.00, 4.60,14,70.00, 0.50, 0.01,2004,2006,2003-12-31,82800,48 
 74,ES-VDA,42.152180, 1.448500,7,30.00, 0.50, 0.29,13,70.00, 0.50, 0.01,2004,2004,2003-12-31,82800,48 
@@ -94,7 +95,7 @@
 89,IE-Ca1,52.858791,-6.918152,15,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2004,2006,2004-01-01,0,48 
 90,IE-Dri,51.986691,-8.751801,13,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2003,2005,2003-01-01,0,48 
 91,IT-Amp,41.904099,13.605160,13,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2003,2006,2002-12-31,82800,48 
-92,IT-BCi,40.523800,14.957440,-999,-999.00,-999.00,-999.00,16,100.00, 0.50, 0.01,2005,2010,2004-12-31,82800,48 
+92,IT-BCi,40.523800,14.957440,16,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2005,2010,2004-12-31,82800,48 
 93,IT-CA1,42.380409,12.026560,7,100.00, 5.50, 3.16,-999,-999.00,-999.00,-999.00,2012,2013,2011-12-31,82800,48 
 94,IT-CA2,42.377220,12.026040,15,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2012,2013,2011-12-31,82800,48 
 95,IT-CA3,42.380001,12.022200,7,100.00, 3.50, 2.01,-999,-999.00,-999.00,-999.00,2012,2013,2011-12-31,82800,48 
@@ -151,8 +152,8 @@
 146,US-MMS,39.323200,-86.413086,7,100.00,27.00,15.52,-999,-999.00,-999.00,-999.00,1999,2014,1999-01-01,18000,24 
 147,US-MOz,38.744110,-92.200012,7,100.00,24.00,13.80,-999,-999.00,-999.00,-999.00,2005,2006,2005-01-01,21600,48 
 148,US-Myb,38.049801,-121.765106,13,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2011,2014,2011-01-01,28800,48 
-149,US-Ne1,41.165100,-96.476593,-999,-999.00,-999.00,-999.00,16,100.00, 0.50, 0.01,2002,2012,2002-01-01,21600,24 
-150,US-Ne2,41.164902,-96.470093,-999,-999.00,-999.00,-999.00,16,100.00, 0.50, 0.01,2002,2012,2002-01-01,21600,24 
+149,US-Ne1,41.165100,-96.476593,16,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2002,2012,2002-01-01,21600,24 
+150,US-Ne2,41.164902,-96.470093,16,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2002,2012,2002-01-01,21600,24 
 151,US-Ne3,41.179699,-96.439697,15,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2002,2012,2002-01-01,21600,24 
 152,US-NR1,40.032902,-105.546402,1,100.00,12.00, 6.00,-999,-999.00,-999.00,-999.00,1999,2014,1999-01-01,25200,48 
 153,US-PFa,45.945900,-90.272308,1, 8.18,30.00,15.00,7,91.82,30.00,17.25,1995,2014,1995-01-01,21600,24 
@@ -165,7 +166,7 @@
 160,US-Syv,46.242001,-89.347717,1, 4.91,27.00,13.50,7,95.09,27.00,15.53,2002,2008,2002-01-01,21600,48 
 161,US-Ton,38.431599,-120.966003,7,70.00, 7.10, 4.08,14,30.00, 0.50, 0.01,2001,2014,2001-01-01,28800,48 
 162,US-Tw4,38.103001,-121.641403,13,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2014,2014,2014-01-01,28800,48 
-163,US-Twt,38.108700,-121.653107,-999,-999.00,-999.00,-999.00,16,100.00, 0.50, 0.01,2010,2014,2010-01-01,28800,48 
+163,US-Twt,38.108700,-121.653107,16,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2010,2014,2010-01-01,28800,48 
 164,US-UMB,45.559799,-84.713806,7,100.00,20.00,11.50,-999,-999.00,-999.00,-999.00,2000,2014,2000-01-01,18000,24 
 165,US-Var,38.413300,-120.950729,14,100.00, 0.50, 0.01,-999,-999.00,-999.00,-999.00,2001,2014,2001-01-01,28800,48 
 166,US-WCr,45.805901,-90.079895,7,100.00,24.00,13.80,-999,-999.00,-999.00,-999.00,1999,2006,1999-01-01,21600,48