diff --git a/cime_config/SystemTests/rxcropmaturity.py b/cime_config/SystemTests/rxcropmaturity.py
index a0eced83c5..589fe38fab 100644
--- a/cime_config/SystemTests/rxcropmaturity.py
+++ b/cime_config/SystemTests/rxcropmaturity.py
@@ -448,6 +448,7 @@ def _run_generate_gdds(self, case_gddgen):
                 f"--hdates-file {hdates_file}",
                 f"--output-dir generate_gdds_out",
                 f"--skip-crops miscanthus,irrigated_miscanthus,switchgrass,irrigated_switchgrass",
+                "--max-season-length-from-hdates-file",
             ]
         )
         stu.run_python_script(
diff --git a/python/ctsm/crop_calendars/cropcal_module.py b/python/ctsm/crop_calendars/cropcal_module.py
index 393e3a9cd5..976d022e85 100644
--- a/python/ctsm/crop_calendars/cropcal_module.py
+++ b/python/ctsm/crop_calendars/cropcal_module.py
@@ -2,8 +2,8 @@
 Helper functions for various crop calendar stuff
 """
 
-import os
-import glob
+import warnings
+
 import numpy as np
 import xarray as xr
 
@@ -203,16 +203,12 @@ def get_gs_len_da(this_da):
     return this_da
 
 
-def import_max_gs_length(paramfile_dir, my_clm_ver, my_clm_subver):
+def import_max_gs_length(paramfile):
     """
     Import maximum growing season length
     """
     # Get parameter file
-    pattern = os.path.join(paramfile_dir, f"*{my_clm_ver}_params.{my_clm_subver}.nc")
-    paramfile = glob.glob(pattern)
-    if len(paramfile) != 1:
-        raise RuntimeError(f"Expected to find 1 match of {pattern}; found {len(paramfile)}")
-    paramfile_ds = xr.open_dataset(paramfile[0])
+    paramfile_ds = xr.open_dataset(paramfile)
 
     # Import max growing season length (stored in netCDF as nanoseconds!)
     paramfile_mxmats = paramfile_ds["mxmat"].values / np.timedelta64(1, "D")
@@ -234,6 +230,56 @@ def import_max_gs_length(paramfile_dir, my_clm_ver, my_clm_subver):
     return mxmat_dict
 
 
+def cushion_gs_length(mxmat_dict_in, cushion, *, min_mxmat=1, max_mxmat=365):
+    """
+    Given a dictionary of maximum growing season lengths, apply a "cushion". This is useful for
+    generating crop maturity requirements: If observed growing season is longer than maximum,
+    and you're limiting growing seasons based on the maximum, then you'd expect about 50% of
+    growing seasons to fail to reach maturity (assuming a normal distribution of seasonal
+    growing degree-days).
+    """
+
+    mxmat_dict_out = mxmat_dict_in.copy()
+
+    for pftname, mxmat in mxmat_dict_in.items():
+        # Skip PFTs without max growing season length
+        if np.isinf(mxmat):
+            continue
+
+        assert mxmat >= min_mxmat, f"{pftname} input mxmat ({mxmat}) is < min_mxmat ({min_mxmat})"
+        assert mxmat <= max_mxmat, f"{pftname} input mxmat ({mxmat}) is > max_mxmat ({max_mxmat})"
+
+        new_mxmat = mxmat - cushion
+
+        # Apply limits
+        msg = None
+        if new_mxmat < min_mxmat:
+            msg = (
+                f"Applying cushion of {cushion} to {pftname}'s mxmat ({mxmat}) resulted in new"
+                f" mxmat of {new_mxmat}; increasing that to min_mxmat {min_mxmat}"
+            )
+            new_mxmat = min_mxmat
+        elif new_mxmat > max_mxmat:
+            msg = (
+                f"Applying cushion of {cushion} to {pftname}'s mxmat ({mxmat}) resulted in new"
+                f" mxmat of {new_mxmat}; decreasing that to max_mxmat {max_mxmat}"
+            )
+            new_mxmat = max_mxmat
+        if msg:
+            warnings.warn(msg, RuntimeWarning)
+
+        assert (
+            new_mxmat >= min_mxmat
+        ), f"{pftname} new_mxmat ({new_mxmat}) < min_mxmat ({min_mxmat})"
+        assert (
+            new_mxmat <= max_mxmat
+        ), f"{pftname} new_mxmat ({new_mxmat}) > max_mxmat ({max_mxmat})"
+
+        mxmat_dict_out[pftname] = new_mxmat
+
+    return mxmat_dict_out
+
+
 def unexpected_negative_rx_gdd(data_array):
     """
     Return True if there's a negative value not matching the designated missing value
diff --git a/python/ctsm/crop_calendars/generate_gdds.py b/python/ctsm/crop_calendars/generate_gdds.py
index 0408792beb..7af82f9fa1 100644
--- a/python/ctsm/crop_calendars/generate_gdds.py
+++ b/python/ctsm/crop_calendars/generate_gdds.py
@@ -26,10 +26,20 @@
 # fixed. For now, we'll just disable the warning.
 # pylint: disable=too-many-positional-arguments
 
-# Global constants
-PARAMFILE_DIR = "/glade/campaign/cesm/cesmdata/cseg/inputdata/lnd/clm2/paramdata"
-MY_CLM_VER = 51
-MY_CLM_SUBVER = "c211112"
+
+def _get_max_growing_season_lengths(max_season_length_from_hdates_file, paramfile, cushion):
+    """
+    Import maximum growing season lengths from paramfile, if doing so.
+    """
+    if max_season_length_from_hdates_file:
+        return None
+
+    mxmats = cc.import_max_gs_length(paramfile)
+
+    if cushion:
+        mxmats = cc.cushion_gs_length(mxmats, cushion)
+
+    return mxmats
 
 
 def main(
@@ -47,10 +57,12 @@ def main(
     land_use_file=None,
     first_land_use_year=None,
     last_land_use_year=None,
-    unlimited_season_length=False,
+    max_season_length_from_hdates_file=False,
     skip_crops=None,
     logger=None,
     no_pickle=None,
+    paramfile=None,
+    max_season_length_cushion=None,
 ):  # pylint: disable=missing-function-docstring,too-many-statements
     # Directories to save output files and figures
     if not output_dir:
@@ -58,7 +70,7 @@ def main(
             output_dir = input_dir
         else:
             output_dir = os.path.join(input_dir, "generate_gdds")
-            if not unlimited_season_length:
+            if not max_season_length_from_hdates_file:
                 output_dir += ".mxmat"
             output_dir += "." + dt.datetime.now().strftime("%Y-%m-%d-%H%M%S")
     if not os.path.exists(output_dir):
@@ -146,10 +158,9 @@ def main(
         sdates_rx = sdates_file
         hdates_rx = hdates_file
 
-        if not unlimited_season_length:
-            mxmats = cc.import_max_gs_length(PARAMFILE_DIR, MY_CLM_VER, MY_CLM_SUBVER)
-        else:
-            mxmats = None
+        mxmats = _get_max_growing_season_lengths(
+            max_season_length_from_hdates_file, paramfile, max_season_length_cushion
+        )
 
         h1_instantaneous = None
         for yr_index, this_yr in enumerate(np.arange(first_season + 1, last_season + 3)):
@@ -390,11 +401,34 @@ def add_attrs_to_map_ds(
         )
 
 
-if __name__ == "__main__":
-    ###############################
-    ### Process input arguments ###
-    ###############################
-    parser = argparse.ArgumentParser(description="ADD DESCRIPTION HERE")
+def _parse_args(argv):
+    parser = argparse.ArgumentParser(
+        description=(
+            "A script to generate maturity requirements for CLM crops in units of growing degree-"
+            "days (GDDs)."
+        )
+    )
+
+    # Required but mutually exclusive
+    max_growing_season_length_group = parser.add_mutually_exclusive_group(required=True)
+    max_growing_season_length_group.add_argument(
+        "--paramfile",
+        help=(
+            "Path to parameter file with maximum growing season lengths (mxmat)."
+            " Mutually exclusive with --max-season-length-from-hdates-file."
+        ),
+    )
+    max_growing_season_length_group.add_argument(
+        "--max-season-length-from-hdates-file",
+        help=(
+            "Rather than limiting growing season length based on mxmat values from a CLM parameter"
+            " file, use the season lengths from --hdates-file. Not recommended unless you use the"
+            "results of this script in a run with sufficiently long mxmat values!"
+            " Mutually exclusive with --paramfile."
+        ),
+        action="store_true",
+        default=False,
+    )
 
     # Required
     parser.add_argument(
@@ -467,12 +501,6 @@ def add_attrs_to_map_ds(
         default=None,
         type=int,
     )
-    parser.add_argument(
-        "--unlimited-season-length",
-        help="Limit mean growing season length based on CLM CFT parameter mxmat.",
-        action="store_true",
-        default=False,
-    )
     parser.add_argument(
         "--skip-crops",
         help="Skip processing of these crops. Comma- or space-separated list.",
@@ -485,12 +513,40 @@ def add_attrs_to_map_ds(
         action="store_true",
         default=False,
     )
+    parser.add_argument(
+        "--max-season-length-cushion",
+        help=(
+            "How much to reduce the maximum growing season length (mxmat) for each crop in the"
+            " parameter file. This might be useful for helping avoid high rates of immature"
+            " harvests for gridcells where the observed harvest date is longer than mxmat."
+            " Incompatible with --max-season-length-from-hdates-file."
+        ),
+        default=0,
+        type=int,
+    )
 
     # Get arguments
-    args = parser.parse_args(sys.argv[1:])
-    for k, v in sorted(vars(args).items()):
+    args_parsed = parser.parse_args(argv)
+    for k, v in sorted(vars(args_parsed).items()):
         print(f"{k}: {v}")
 
+    # Check arguments
+    if args_parsed.max_season_length_from_hdates_file and args_parsed.max_season_length_cushion:
+        raise argparse.ArgumentError(
+            None,
+            "--max-season-length-from-hdates-file is incompatible with --max-season-length-cushion"
+            " ≠ 0.",
+        )
+
+    return args_parsed
+
+
+if __name__ == "__main__":
+    ###############################
+    ### Process input arguments ###
+    ###############################
+    args = _parse_args(sys.argv[1:])
+
     # Call main()
     main(
         input_dir=args.input_dir,
@@ -506,7 +562,9 @@ def add_attrs_to_map_ds(
         land_use_file=args.land_use_file,
         first_land_use_year=args.first_land_use_year,
         last_land_use_year=args.last_land_use_year,
-        unlimited_season_length=args.unlimited_season_length,
+        max_season_length_from_hdates_file=args.max_season_length_from_hdates_file,
         skip_crops=args.skip_crops,
         no_pickle=args.no_pickle,
+        paramfile=args.paramfile,
+        max_season_length_cushion=args.max_season_length_cushion,
     )
diff --git a/python/ctsm/crop_calendars/generate_gdds_functions.py b/python/ctsm/crop_calendars/generate_gdds_functions.py
index ec60a85beb..4f6dd6b966 100644
--- a/python/ctsm/crop_calendars/generate_gdds_functions.py
+++ b/python/ctsm/crop_calendars/generate_gdds_functions.py
@@ -22,6 +22,9 @@
 # fixed. For now, we'll just disable the warning.
 # pylint: disable=too-many-positional-arguments
 
+# Tolerance (degrees) for checking lat/lon grid matches
+GRID_TOL_DEG = 1e-6
+
 CAN_PLOT = True
 try:
     # pylint: disable=wildcard-import,unused-wildcard-import
@@ -54,6 +57,31 @@
     CAN_PLOT = False
 
 
+def check_grid_match(grid0, grid1, tol=GRID_TOL_DEG):
+    """
+    Check whether latitude or longitude values match
+    """
+    if grid0.shape != grid1.shape:
+        return False, None
+
+    if hasattr(grid0, "values"):
+        grid0 = grid0.values
+    if hasattr(grid1, "values"):
+        grid1 = grid1.values
+
+    abs_diff = np.abs(grid1 - grid0)
+    if np.any(np.isnan(abs_diff)):
+        if np.any(np.isnan(grid0) != np.isnan(grid1)):
+            warnings.warn("NaN(s) in grid don't match", RuntimeWarning)
+            return False, None
+        warnings.warn("NaN(s) in grid", RuntimeWarning)
+
+    max_abs_diff = np.nanmax(abs_diff)
+    match = max_abs_diff < tol
+
+    return match, max_abs_diff
+
+
 def check_sdates(dates_ds, sdates_rx, outdir_figs, logger, verbose=False):
     """
     Checking that input and output sdates match
@@ -62,6 +90,29 @@ def check_sdates(dates_ds, sdates_rx, outdir_figs, logger, verbose=False):
 
     sdates_grid = grid_one_variable(dates_ds, "SDATES")
 
+    # In this script, we assume that you used prescribed sowing dates on the same grid as the
+    # CLM run
+    # Check that latitudes match, coercing a match if difference is acceptable
+    match, max_abs_diff = check_grid_match(sdates_rx["lat"], sdates_grid["lat"])
+    assert bool(
+        match
+    ), f"CLM lat grid doesn't match rx sdates's within {GRID_TOL_DEG}; max abs diff {max_abs_diff}"
+    if max_abs_diff > 0:
+        log(logger, f"Max lat abs diff: {max_abs_diff}. Coercing match.")
+        # sdates_grid comes from the CLM outputs, so we coerce the prescribed sowing date file
+        # coordinate to match, because we want the outputs of this script to have CLM's coordinates.
+        sdates_rx["lat"] = sdates_grid["lat"]
+    # Check that longitudes match, coercing a match if difference is acceptable
+    match, max_abs_diff = check_grid_match(sdates_rx["lon"], sdates_grid["lon"])
+    assert bool(
+        check_grid_match(sdates_rx["lon"], sdates_grid["lon"])
+    ), f"CLM lon grid doesn't match rx sdates's within {GRID_TOL_DEG}; max abs diff {max_abs_diff}"
+    if max_abs_diff > 0:
+        log(logger, f"Max lon abs diff: {max_abs_diff}. Coercing match.")
+        # sdates_grid comes from the CLM outputs, so we coerce the prescribed sowing date file
+        # coordinate to match, because we want the outputs of this script to have CLM's coordinates.
+        sdates_rx["lon"] = sdates_grid["lon"]
+
     all_ok = True
     any_found = False
     vegtypes_skipped = []
@@ -83,8 +134,16 @@ def check_sdates(dates_ds, sdates_rx, outdir_figs, logger, verbose=False):
         # Output
         out_map = sdates_grid.sel(ivt_str=vegtype_str).squeeze(drop=True)
 
-        # Check for differences
+        # Calculate differences
         diff_map = out_map - in_map
+        assert (
+            diff_map.shape == in_map.shape
+        ), f"Diff map shape {diff_map.shape} doesn't match in_map shape {in_map.shape}"
+        assert (
+            diff_map.shape == out_map.shape
+        ), f"Diff map shape {diff_map.shape} doesn't match out_map shape {out_map.shape}"
+
+        # Check for differences
         diff_map_notnan = diff_map.values[np.invert(np.isnan(diff_map.values))]
         if np.any(diff_map_notnan):
             log(logger, f"Difference(s) found in {vegtype_str}")
@@ -496,6 +555,29 @@ def import_and_process_1yr(
         "h", hdates_rx, incl_patches1d_itype_veg, mxsowings, logger
     )  # Yes, mxsowings even when importing harvests
 
+    # In this script, we assume that you have prescribed harvest dates on the same grid as the
+    # CLM run
+    # Check that latitudes match, coercing a match if difference is acceptable
+    match, max_abs_diff = check_grid_match(hdates_rx_orig["lat"], dates_incl_ds["lat"])
+    assert bool(
+        match
+    ), f"CLM lat grid doesn't match rx hdates's within {GRID_TOL_DEG}; max abs diff {max_abs_diff}"
+    if max_abs_diff > 0:
+        log(logger, f"Max lat abs diff: {max_abs_diff}. Coercing match.")
+        # dates_incl_ds comes from the CLM outputs, so we coerce the prescribed harvest date file
+        # coordinate to match, because we want the outputs of this script to have CLM's coordinates.
+        hdates_rx_orig["lat"] = dates_incl_ds["lat"]
+    # Check that longitudes match, coercing a match if difference is acceptable
+    match, max_abs_diff = check_grid_match(hdates_rx_orig["lon"], dates_incl_ds["lon"])
+    assert bool(
+        match
+    ), f"CLM lon grid doesn't match rx hdates's within {GRID_TOL_DEG}; max abs diff {max_abs_diff}"
+    if max_abs_diff > 0:
+        log(logger, f"Max lon abs diff: {max_abs_diff}. Coercing match.")
+        # dates_incl_ds comes from the CLM outputs, so we coerce the prescribed harvest date file
+        # coordinate to match, because we want the outputs of this script to have CLM's coordinates.
+        hdates_rx_orig["lon"] = dates_incl_ds["lon"]
+
     # Limit growing season to CLM max growing season length, if needed
     if mxmats and (imported_sdates or imported_hdates):
         print("   Limiting growing season length...")
diff --git a/python/ctsm/test/test_sys_cropcal_module.py b/python/ctsm/test/test_sys_cropcal_module.py
new file mode 100755
index 0000000000..5614fbc211
--- /dev/null
+++ b/python/ctsm/test/test_sys_cropcal_module.py
@@ -0,0 +1,45 @@
+#!/usr/bin/env python3
+
+"""
+System tests for cropcal_module.py
+"""
+
+import unittest
+import os
+
+from ctsm import unit_testing
+from ctsm.crop_calendars.cropcal_module import import_max_gs_length
+
+# Allow test names that pylint doesn't like; otherwise hard to make them
+# readable
+# pylint: disable=invalid-name
+
+# pylint: disable=protected-access
+
+## Too many instant variables as part of the class (too many self.<varible> in the SetUp)
+# pylint: disable=too-many-instance-attributes
+
+
+class TestImportMaxGsLength(unittest.TestCase):
+    """Tests of import_max_gs_length()"""
+
+    def setUp(self):
+        self._paramfile_51 = os.path.join(
+            os.path.dirname(__file__), "testinputs", "ctsm51_params.c211112.nc"
+        )
+        self._paramfile_60 = os.path.join(
+            os.path.dirname(__file__), "testinputs", "ctsm60_params_cal115_c250813.nc"
+        )
+
+    def test_import_max_gs_length_ctsm51(self):
+        """Basic test with ctsm51 paramfile"""
+        import_max_gs_length(self._paramfile_51)
+
+    def test_import_max_gs_length_ctsm60(self):
+        """Basic test with ctsm60 paramfile"""
+        import_max_gs_length(self._paramfile_60)
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_unit_cropcal_module.py b/python/ctsm/test/test_unit_cropcal_module.py
new file mode 100755
index 0000000000..cff339fdcb
--- /dev/null
+++ b/python/ctsm/test/test_unit_cropcal_module.py
@@ -0,0 +1,125 @@
+#!/usr/bin/env python3
+
+"""Unit tests for cropcal_module.py"""
+
+import unittest
+
+import numpy as np
+
+from ctsm import unit_testing
+from ctsm.crop_calendars import cropcal_module as ccu
+
+# Allow names that pylint doesn't like, because otherwise it's hard
+# to make readable unit test names
+# pylint: disable=invalid-name
+
+
+class TestCushionGsLength(unittest.TestCase):
+    """Tests of cushion_gs_length()"""
+
+    def setUp(self):
+        self._mxmat_dict = {
+            "cropA": 55,
+            "cropB": 57,
+            "cropC": 87,
+            "cropD": 86,
+        }
+
+    def test_cushion_gs_length_1(self):
+        """
+        Test cushion_gs_length() for a cushion of 1.
+        """
+        cushion = 1
+        result = ccu.cushion_gs_length(self._mxmat_dict, cushion)
+        for crop, mxmat_orig in self._mxmat_dict.items():
+            self.assertEqual(result[crop], mxmat_orig - cushion)
+
+    def test_cushion_gs_length_neg1(self):
+        """
+        Test cushion_gs_length() for a cushion of -1.
+        """
+        cushion = -1
+        result = ccu.cushion_gs_length(self._mxmat_dict, cushion)
+        for crop, mxmat_orig in self._mxmat_dict.items():
+            self.assertEqual(result[crop], mxmat_orig - cushion)
+
+    def test_cushion_gs_length_inf(self):
+        """
+        As test_cushion_gs_length_1 but with an infinite value in original.
+        """
+        self._mxmat_dict["cropB"] = np.inf
+        cushion = 1
+        result = ccu.cushion_gs_length(self._mxmat_dict, cushion)
+        for crop, mxmat_orig in self._mxmat_dict.items():
+            if np.isinf(mxmat_orig):
+                self.assertTrue(np.isinf(result[crop]))
+            else:
+                self.assertEqual(result[crop], mxmat_orig - cushion)
+
+    def test_cushion_gs_length_limit_range_default_min1(self):
+        """
+        Test that cushion_gs_length() limits output range to a minimum of 1 by default
+        """
+        cushion = np.inf
+
+        with self.assertWarnsRegex(RuntimeWarning, "increasing that to"):
+            result = ccu.cushion_gs_length(self._mxmat_dict, cushion)
+
+        self.assertTrue(all(x == 1 for x in result.values()))
+
+    def test_cushion_gs_length_limit_range_default_max365(self):
+        """
+        Test that cushion_gs_length() limits output range to a maximum of 365 by default
+        """
+        cushion = -np.inf
+
+        with self.assertWarnsRegex(RuntimeWarning, "decreasing that to"):
+            result = ccu.cushion_gs_length(self._mxmat_dict, cushion)
+
+        print(result)
+        self.assertTrue(all(x == 365 for x in result.values()))
+
+    def test_cushion_gs_length_limit_range_custom_min(self):
+        """
+        Test cushion_gs_length() with custom min
+        """
+        cushion = np.inf
+        min_mxmat = -1
+
+        with self.assertWarnsRegex(RuntimeWarning, "increasing that to"):
+            result = ccu.cushion_gs_length(self._mxmat_dict, cushion, min_mxmat=min_mxmat)
+
+        self.assertTrue(all(x == min_mxmat for x in result.values()))
+
+    def test_cushion_gs_length_limit_range_custom_max(self):
+        """
+        Test cushion_gs_length() with custom max
+        """
+        cushion = -np.inf
+        max_mxmat = 400
+
+        with self.assertWarnsRegex(RuntimeWarning, "decreasing that to"):
+            result = ccu.cushion_gs_length(self._mxmat_dict, cushion, max_mxmat=max_mxmat)
+
+        self.assertTrue(all(x == max_mxmat for x in result.values()))
+
+    def test_cushion_gs_length_assert_orig_too_low(self):
+        """
+        Test that cushion_gs_length() errors if original value is too low
+        """
+        self._mxmat_dict["cropC"] = -999
+        with self.assertRaisesRegex(AssertionError, "is < min_mxmat"):
+            ccu.cushion_gs_length(self._mxmat_dict, 1)
+
+    def test_cushion_gs_length_assert_orig_too_high(self):
+        """
+        Test that cushion_gs_length() errors if original value is too high
+        """
+        self._mxmat_dict["cropC"] = 999
+        with self.assertRaisesRegex(AssertionError, "is > max_mxmat"):
+            ccu.cushion_gs_length(self._mxmat_dict, 1)
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_unit_generate_gdds.py b/python/ctsm/test/test_unit_generate_gdds.py
new file mode 100755
index 0000000000..4976097b7d
--- /dev/null
+++ b/python/ctsm/test/test_unit_generate_gdds.py
@@ -0,0 +1,330 @@
+#!/usr/bin/env python3
+
+"""
+Unit tests for generate_gdds.py and generate_gdds_functions.py
+"""
+
+import unittest
+import os
+import argparse
+
+import numpy as np
+import xarray as xr
+
+from ctsm import unit_testing
+from ctsm.crop_calendars import generate_gdds as gg
+from ctsm.crop_calendars import generate_gdds_functions as gf
+
+# Allow test names that pylint doesn't like; otherwise hard to make them
+# readable
+# pylint: disable=invalid-name
+
+# pylint: disable=protected-access
+
+## Too many instant variables as part of the class (too many self.<varible> in the SetUp)
+# pylint: disable=too-many-instance-attributes
+
+
+class TestGenerateGddsArgs(unittest.TestCase):
+    """Tests the generate_gdds.py argument parsing"""
+
+    def setUp(self):
+        self._input_dir = os.path.join("dummy", "path", "to", "inputdir")
+        self._sdates_file = os.path.join("dummy", "path", "to", "sdates")
+        self._hdates_file = os.path.join("dummy", "path", "to", "hdates")
+        self._paramfile = os.path.join("dummy", "path", "to", "paramfile")
+
+    def test_generate_gdds_args_reqd_shortnames(self):
+        """Basic test with all required inputs, short arg names"""
+        args = [
+            "-i",
+            self._input_dir,
+            "-1",
+            "1986",
+            "-n",
+            "1987",
+            "-sd",
+            self._sdates_file,
+            "-hd",
+            self._hdates_file,
+            "--paramfile",
+            self._paramfile,
+        ]
+        gg._parse_args(args)
+
+        # Again, with capital -N option
+        args = [
+            "-i",
+            self._input_dir,
+            "-1",
+            "1986",
+            "-N",
+            "1987",
+            "-sd",
+            self._sdates_file,
+            "-hd",
+            self._hdates_file,
+            "--paramfile",
+            self._paramfile,
+        ]
+        gg._parse_args(args)
+
+    def test_generate_gdds_args_reqd_longnames(self):
+        """Basic test with all required inputs, long arg names"""
+        args = [
+            "--input-dir",
+            self._input_dir,
+            "--first-season",
+            "1986",
+            "--last-season",
+            "1987",
+            "--sdates-file",
+            self._sdates_file,
+            "--hdates-file",
+            self._hdates_file,
+            "--paramfile",
+            self._paramfile,
+        ]
+        gg._parse_args(args)
+
+    def test_generate_gdds_args_mxmat_from_hdatefile(self):
+        """Test with option to get max season length from hdates file"""
+        args = [
+            "--input-dir",
+            self._input_dir,
+            "--first-season",
+            "1986",
+            "--last-season",
+            "1987",
+            "--sdates-file",
+            self._sdates_file,
+            "--hdates-file",
+            self._hdates_file,
+            "--max-season-length-from-hdates-file",
+        ]
+        gg._parse_args(args)
+
+    def test_generate_gdds_args_error_with_paramfile_and_nomxmat(self):
+        """Should error if both --paramfile and --max-season-length-from-hdates-file are given"""
+        args = [
+            "--input-dir",
+            self._input_dir,
+            "--first-season",
+            "1986",
+            "--last-season",
+            "1987",
+            "--sdates-file",
+            self._sdates_file,
+            "--hdates-file",
+            self._hdates_file,
+            "--max-season-length-from-hdates-file",
+            "--paramfile",
+            self._paramfile,
+        ]
+        with self.assertRaises(SystemExit):
+            gg._parse_args(args)
+
+    def test_generate_gdds_args_error_with_nomxmat_and_cushion(self):
+        """Should error if both --max-season-length-cushion and --max-season-length-from-hdates-file
+        are given"""
+        args = [
+            "--input-dir",
+            self._input_dir,
+            "--first-season",
+            "1986",
+            "--last-season",
+            "1987",
+            "--sdates-file",
+            self._sdates_file,
+            "--hdates-file",
+            self._hdates_file,
+            "--max-season-length-from-hdates-file",
+            "--max-season-length-cushion",
+            "14",
+        ]
+        with self.assertRaises(argparse.ArgumentError):
+            gg._parse_args(args)
+
+    def test_generate_gdds_args_ok_with_nomxmat_and_cushion0(self):
+        """As test_generate_gdds_args_error_with_nomxmat_and_cushion, but cushion 0 is ok"""
+        args = [
+            "--input-dir",
+            self._input_dir,
+            "--first-season",
+            "1986",
+            "--last-season",
+            "1987",
+            "--sdates-file",
+            self._sdates_file,
+            "--hdates-file",
+            self._hdates_file,
+            "--max-season-length-from-hdates-file",
+            "--max-season-length-cushion",
+            "0",
+        ]
+        gg._parse_args(args)
+
+
+class TestGetMaxGsLengths(unittest.TestCase):
+    """Tests get_max_growing_season_lengths()"""
+
+    def setUp(self):
+        self._paramfile_51 = os.path.join(
+            os.path.dirname(__file__), "testinputs", "ctsm51_params.c211112.nc"
+        )
+        self._paramfile_60 = os.path.join(
+            os.path.dirname(__file__), "testinputs", "ctsm60_params_cal115_c250813.nc"
+        )
+
+        # Default arguments
+        self.no_mxmats = False
+        self.paramfile = self._paramfile_60
+        self.cushion = 0
+
+    def test_generate_gdds_get_mxmats_ctsm51(self):
+        """Test importing from a ctsm51 paramfile (no fail)"""
+        paramfile = self._paramfile_51
+        gg._get_max_growing_season_lengths(self.no_mxmats, paramfile, self.cushion)
+
+    def test_generate_gdds_get_mxmats_ctsm60(self):
+        """Test importing from a ctsm60 paramfile (no fail)"""
+        paramfile = self._paramfile_60
+        gg._get_max_growing_season_lengths(self.no_mxmats, paramfile, self.cushion)
+
+    def test_generate_gdds_get_mxmats_none(self):
+        """Test not importing from a paramfile (should return None)"""
+        max_season_length_from_hdates_file = True
+        paramfile = None
+        result = gg._get_max_growing_season_lengths(
+            max_season_length_from_hdates_file, paramfile, self.cushion
+        )
+        self.assertIsNone(result)
+
+    def test_generate_gdds_get_mxmats_values(self):
+        """Check values with no cushion"""
+        mxmats = gg._get_max_growing_season_lengths(self.no_mxmats, self.paramfile, self.cushion)
+
+        # Check values
+        self.assertTrue(np.isinf(mxmats["needleleaf_evergreen_temperate_tree"]))
+        self.assertEqual(mxmats["temperate_corn"], 165)
+        self.assertEqual(mxmats["miscanthus"], 210)
+
+    def test_generate_gdds_get_mxmats_cushion14(self):
+        """As test_generate_gdds_get_mxmats_values, but cushion 14"""
+        cushion = 14
+        mxmats = gg._get_max_growing_season_lengths(self.no_mxmats, self.paramfile, cushion)
+
+        # Check values
+        self.assertTrue(np.isinf(mxmats["needleleaf_evergreen_temperate_tree"]))
+        self.assertEqual(mxmats["temperate_corn"], 165 - cushion)
+        self.assertEqual(mxmats["miscanthus"], 210 - cushion)
+
+    def test_generate_gdds_get_mxmats_cushionneg14(self):
+        """As test_generate_gdds_get_mxmats_values, but cushion -14"""
+        cushion = -14
+        mxmats = gg._get_max_growing_season_lengths(self.no_mxmats, self.paramfile, cushion)
+
+        # Check values
+        self.assertTrue(np.isinf(mxmats["needleleaf_evergreen_temperate_tree"]))
+        self.assertEqual(mxmats["temperate_corn"], 165 - cushion)
+        self.assertEqual(mxmats["miscanthus"], 210 - cushion)
+
+
+class TestCheckGridMatch(unittest.TestCase):
+    """Tests check_grid_match()"""
+
+    def test_check_grid_match_true_npnp(self):
+        """Test check_grid_match() with two matching numpy arrays"""
+        np0 = np.array([0, 1, 2, np.pi])
+        match, max_abs_diff = gf.check_grid_match(np0, np0)
+        self.assertTrue(match)
+        self.assertEqual(max_abs_diff, 0.0)
+
+    def test_check_grid_match_true_dada(self):
+        """Test check_grid_match() with two matching DataArrays"""
+        np0 = np.array([0, 1, 2, np.pi])
+        da0 = xr.DataArray(data=np0)
+        match, max_abs_diff = gf.check_grid_match(da0, da0)
+        self.assertTrue(match)
+        self.assertEqual(max_abs_diff, 0.0)
+
+    def test_check_grid_match_false_npnp(self):
+        """Test check_grid_match() with two non-matching numpy arrays"""
+        np0 = np.array([0, 1, 2, np.pi])
+        np1 = np0.copy()
+        diff = 2 * gf.GRID_TOL_DEG
+        np1[0] = np0[0] + diff
+        match, max_abs_diff = gf.check_grid_match(np0, np1)
+        self.assertFalse(match)
+        self.assertEqual(max_abs_diff, diff)
+
+    def test_check_grid_match_false_dada(self):
+        """Test check_grid_match() with two non-matching DataArrays"""
+        np0 = np.array([0, 1, 2, np.pi])
+        np1 = np0.copy()
+        diff = 2 * gf.GRID_TOL_DEG
+        np1[0] = np0[0] + diff
+        da0 = xr.DataArray(data=np0)
+        da1 = xr.DataArray(data=np1)
+        match, max_abs_diff = gf.check_grid_match(da0, da1)
+        self.assertFalse(match)
+        self.assertEqual(max_abs_diff, diff)
+
+    def test_check_grid_match_falseneg_npnp(self):
+        """As test_check_grid_match_false_npnp, but with diff in negative direction"""
+        np0 = np.array([0, 1, 2, np.pi])
+        np1 = np0.copy()
+        diff = -2 * gf.GRID_TOL_DEG
+        np1[0] = np0[0] + diff
+        match, max_abs_diff = gf.check_grid_match(np0, np1)
+        self.assertFalse(match)
+        self.assertEqual(max_abs_diff, abs(diff))
+
+    def test_check_grid_match_matchnans_true_npnp(self):
+        """Test check_grid_match() with two numpy arrays that have nans and match"""
+        np0 = np.array([np.nan, 1, 2, np.pi])
+        with self.assertWarnsRegex(RuntimeWarning, r"NaN\(s\) in grid"):
+            match, max_abs_diff = gf.check_grid_match(np0, np0)
+        self.assertTrue(match)
+        self.assertEqual(max_abs_diff, 0.0)
+
+    def test_check_grid_match_matchnans_true_dada(self):
+        """Test check_grid_match() with two DataArrays that have nans and match"""
+        np0 = np.array([np.nan, 1, 2, np.pi])
+        da0 = xr.DataArray(data=np0)
+        with self.assertWarnsRegex(RuntimeWarning, r"NaN\(s\) in grid"):
+            match, max_abs_diff = gf.check_grid_match(da0, da0)
+        self.assertTrue(match)
+        self.assertEqual(max_abs_diff, 0.0)
+
+    def test_check_grid_match_matchnans_false_npnp(self):
+        """Test check_grid_match() with two numpy arrays with nans that DON'T match"""
+        np0 = np.array([np.nan, 1, 2, np.pi])
+        np1 = np.array([np.nan, 1, np.nan, np.pi])
+        with self.assertWarnsRegex(RuntimeWarning, r"NaN\(s\) in grid don't match"):
+            match, max_abs_diff = gf.check_grid_match(np0, np1)
+        self.assertFalse(match)
+        self.assertIsNone(max_abs_diff)
+
+    def test_check_grid_match_matchnans_falseshape_npnp(self):
+        """Test check_grid_match() with two numpy arrays that have different shapes"""
+        np0 = np.array([0, 1, 2, np.pi])
+        np1 = np.array([0, 1, 2, np.pi, 4])
+        match, max_abs_diff = gf.check_grid_match(np0, np1)
+        self.assertFalse(match)
+        self.assertIsNone(max_abs_diff)
+
+    def test_check_grid_match_matchnans_falseshape_dada(self):
+        """Test check_grid_match() with two DataArrays that have different shapes"""
+        np0 = np.array([0, 1, 2, np.pi])
+        np1 = np.array([0, 1, 2, np.pi, 4])
+        da0 = xr.DataArray(data=np0)
+        da1 = xr.DataArray(data=np1)
+        match, max_abs_diff = gf.check_grid_match(da0, da1)
+        self.assertFalse(match)
+        self.assertIsNone(max_abs_diff)
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/testinputs/ctsm51_params.c211112.nc b/python/ctsm/test/testinputs/ctsm51_params.c211112.nc
new file mode 100644
index 0000000000..ce6c2676e2
--- /dev/null
+++ b/python/ctsm/test/testinputs/ctsm51_params.c211112.nc
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e1cf2f3455e4cee4ea34459ef28e57f546bea595f865300b578d889afdb208cf
+size 228464
diff --git a/python/ctsm/test/testinputs/ctsm60_params_cal115_c250813.nc b/python/ctsm/test/testinputs/ctsm60_params_cal115_c250813.nc
new file mode 100644
index 0000000000..6fceb6913d
--- /dev/null
+++ b/python/ctsm/test/testinputs/ctsm60_params_cal115_c250813.nc
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7ca369ca10d893f49c4d5223b0cd705dfab699f8e1f6e29d55d878e6dcb8c22d
+size 224108