diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs
index 9badbdf673..e56311635e 100644
--- a/.git-blame-ignore-revs
+++ b/.git-blame-ignore-revs
@@ -70,5 +70,8 @@ cdf40d265cc82775607a1bf25f5f527bacc97405
 4ad46f46de7dde753b4653c15f05326f55116b73
 75db098206b064b8b7b2a0604d3f0bf8fdb950cc
 84609494b54ea9732f64add43b2f1dd035632b4c
+7eb17f3ef0b9829fb55e0e3d7f02e157b0e41cfb
+62d7711506a0fb9a3ad138ceceffbac1b79a6caa
+49ad0f7ebe0b07459abc00a5c33c55a646f1e7e0
 ac03492012837799b7111607188acff9f739044a
 d858665d799690d73b56bcb961684382551193f4
diff --git a/cime_config/SystemTests/setparamfile.py b/cime_config/SystemTests/setparamfile.py
new file mode 100644
index 0000000000..f106a9505a
--- /dev/null
+++ b/cime_config/SystemTests/setparamfile.py
@@ -0,0 +1,95 @@
+"""
+CTSM-specific test that first runs the set_paramfile tool and then ensures that CTSM does not fail
+using the just-generated parameter file
+"""
+
+import os
+import sys
+import logging
+import re
+from CIME.SystemTests.system_tests_common import SystemTestsCommon
+
+# In case we need to import set_paramfile later
+_CTSM_PYTHON = os.path.join(
+    os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir, "python"
+)
+sys.path.insert(1, _CTSM_PYTHON)
+
+logger = logging.getLogger(__name__)
+
+
+class SETPARAMFILE(SystemTestsCommon):
+    def __init__(self, case):
+        """
+        initialize an object interface to the SMS system test
+        """
+        SystemTestsCommon.__init__(self, case)
+
+        # Create out-of-the-box lnd_in to obtain paramfile
+        case.create_namelists(component="lnd")
+
+        # Find the paramfile to modify
+        lnd_in_path = os.path.join(self._get_caseroot(), "CaseDocs", "lnd_in")
+        self._paramfile_in = None
+        with open(lnd_in_path, "r", encoding="utf-8") as lnd_in:
+            for line in lnd_in:
+                paramfile_in = re.match(r" *paramfile *= *'(.*)'", line)
+                if paramfile_in:
+                    self._paramfile_in = paramfile_in.group(1)
+                    break
+        if not self._paramfile_in:
+            raise RuntimeError(f"paramfile not found in {lnd_in_path}")
+
+        # Get the output file
+        self.paramfile_out = os.path.join(self._get_caseroot(), "paramfile.nc")
+
+        # Define set_paramfile command
+        self.set_paramfile_cmd = [
+            "set_paramfile",
+            "-i",
+            self._paramfile_in,
+            "-o",
+            self.paramfile_out,
+            # Change two parameters for one PFT
+            "-p",
+            "needleleaf_deciduous_boreal_tree",
+            "rswf_min=0.35",
+            "rswf_max=0.7",
+        ]
+
+    def build_phase(self, sharedlib_only=False, model_only=False):
+        """
+        Run set_paramfile and then build the model
+        """
+
+        # Run set_paramfile.
+        # build_phase gets called twice:
+        # - once with sharedlib_only = True and
+        # - once with model_only = True
+        # Because we only need set_paramfile run once, we only do it for the sharedlib_only call.
+        # We could also check for the existence of the set_paramfile outputs, but that might lead to
+        # a situation where the user expects set_paramfile to be called but it's not. Better to run
+        # unnecessarily (e.g., if you fixed some FORTRAN code and just need to rebuild).
+        if sharedlib_only:
+            self._run_set_paramfile()
+
+        # Do the build
+        self.build_indv(sharedlib_only=sharedlib_only, model_only=model_only)
+
+    def _run_set_paramfile(self):
+        """
+        Run set_paramfile
+        """
+        # Import set_paramfile. Do it here rather than at top because otherwise the import will
+        # be attempted even during RUN phase.
+        # pylint: disable=wrong-import-position,import-outside-toplevel
+        from ctsm.param_utils.set_paramfile import main as set_paramfile
+
+        # Run set_paramfile
+        sys.argv = self.set_paramfile_cmd
+        set_paramfile()
+
+        # Append
+        user_nl_clm_path = os.path.join(self._get_caseroot(), "user_nl_clm")
+        with open(user_nl_clm_path, "a", encoding="utf-8") as user_nl_clm:
+            user_nl_clm.write(f"paramfile = '{self.paramfile_out}'\n")
diff --git a/cime_config/config_tests.xml b/cime_config/config_tests.xml
index 12859b9131..6df0f357e4 100644
--- a/cime_config/config_tests.xml
+++ b/cime_config/config_tests.xml
@@ -155,6 +155,16 @@ This defines various CTSM-specific system tests
     <HIST_N>$STOP_N</HIST_N>
   </test>
 
+  <test NAME="SETPARAMFILE">
+    <DESC>Modify a copy of the paramfile and run with it.</DESC>
+    <INFO_DBUG>1</INFO_DBUG>
+    <DOUT_S>FALSE</DOUT_S>
+    <CONTINUE_RUN>FALSE</CONTINUE_RUN>
+    <REST_OPTION>never</REST_OPTION>
+    <HIST_OPTION>$STOP_OPTION</HIST_OPTION>
+    <HIST_N>$STOP_N</HIST_N>
+  </test>
+
 <!--
 SSP    smoke CLM spinup test (only valid for CLM compsets with CLM45)
        do an initial spin test (setting CLM_ACCELERATED_SPINUP to on)
diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml
index 478dc59bdd..9c60ed39fc 100644
--- a/cime_config/testdefs/testlist_clm.xml
+++ b/cime_config/testdefs/testlist_clm.xml
@@ -4495,5 +4495,15 @@
     </options>
   </test>
 
+  <test name="SETPARAMFILE_Ld5" grid="f10_f10_mg37" compset="I1850Clm60BgcCrujra" testmods="clm/default">
+    <machines>
+      <machine name="derecho" compiler="gnu" category="aux_clm"/>
+      <machine name="derecho" compiler="gnu" category="clm_pymods"/>
+    </machines>
+    <options>
+      <option name="wallclock">00:20:00</option>
+    </options>
+  </test>
+
 
 </testlist>
diff --git a/doc/source/users_guide/using-clm-tools/index.rst b/doc/source/users_guide/using-clm-tools/index.rst
index 9b5318e5a2..b721b3e7d6 100644
--- a/doc/source/users_guide/using-clm-tools/index.rst
+++ b/doc/source/users_guide/using-clm-tools/index.rst
@@ -22,3 +22,4 @@ Using CLM tools
    creating-domain-files.rst
    observational-sites-datasets.rst
    cprnc.rst
+   paramfile-tools.md
diff --git a/doc/source/users_guide/using-clm-tools/paramfile-tools.md b/doc/source/users_guide/using-clm-tools/paramfile-tools.md
new file mode 100644
index 0000000000..71d881e362
--- /dev/null
+++ b/doc/source/users_guide/using-clm-tools/paramfile-tools.md
@@ -0,0 +1,68 @@
+
+# Tools for working with parameter files
+
+This guide describes the features and usage of the `query_paramfile` and `set_paramfile` tools, located in `tools/param_utils/`. These utilities help users inspect and modify CLM parameter files.
+
+Note that you need to have the `ctsm_pylib` conda environment activated to use these tools. See Sect. :numref:`using-ctsm-pylib` for more information.
+
+## `query_paramfile`
+**Purpose:** Print the values of one or more parameters from a CTSM parameter file (NetCDF format).
+
+**Features:**
+- Print values for specified parameters or all.
+- Optionally filter output by Plant Functional Types (PFTs) for PFT-specific parameters.
+
+For more information, do `tools/param_utils/query_paramfile --help`.
+
+
+### Example usage
+
+Print all variables in a parameter file:
+```bash
+tools/param_utils/query_paramfile -i paramfile.nc
+```
+
+Print specific variables:
+```bash
+tools/param_utils/query_paramfile -i paramfile.nc jmaxha jmaxhd
+```
+
+Print values for specific PFTs:
+```bash
+tools/param_utils/query_paramfile -i paramfile.nc -p needleleaf_evergreen_temperate_tree,c4_grass medlynintercept medlynslope
+```
+
+## `set_paramfile`
+**Purpose:** Change values of one or more parameters in a CTSM parameter file (NetCDF format).
+
+**Features:**
+- Modify parameter values for all or selected PFTs.
+- Optionally drop PFTs not specified.
+- Set parameter values to fill (missing) values using `nan`.
+- Ensures safe file handling and checks for argument validity.
+
+Note that the output file must not already exist.
+
+For more information, do `tools/param_utils/set_paramfile --help`.
+
+### Example usage
+
+Change a scalar parameter:
+```bash
+tools/param_utils/set_paramfile -i paramfile.nc -o output.nc jmaxha=51000
+```
+
+Change a one-dimensional parameter (`mimics_fmet` has the `segment` dimension, length 4):
+```bash
+tools/param_utils/set_paramfile -i paramfile.nc -o output.nc mimics_fmet=0.1,0.2,0.3,0.4
+```
+
+Change a parameter for specific PFTs:
+```bash
+tools/param_utils/set_paramfile -i paramfile.nc -o output.nc -p needleleaf_evergreen_temperate_tree,c4_grass medlynintercept=99.9,100.1 medlynslope=2.99,1.99
+```
+
+Set a parameter to the fill value:
+```bash
+tools/param_utils/set_paramfile -i paramfile.nc -o output.nc -p needleleaf_evergreen_temperate_tree,c4_grass fleafcn=nan,nan
+```
diff --git a/python/ctsm/args_utils.py b/python/ctsm/args_utils.py
index 612e3f09a1..ed0e934d32 100644
--- a/python/ctsm/args_utils.py
+++ b/python/ctsm/args_utils.py
@@ -46,3 +46,12 @@ def plon_type(plon):
             "ERROR: Longitude should be between 0 and 360 or -180 and 180."
         )
     return plon_float
+
+
+def comma_separated_list(value):
+    """
+    Helper function for argparse to split comma-separated strings into a list.
+    """
+    if value is None:
+        return None
+    return [v.strip() for v in value.split(",")]
diff --git a/python/ctsm/netcdf_utils.py b/python/ctsm/netcdf_utils.py
new file mode 100644
index 0000000000..025c7f73b6
--- /dev/null
+++ b/python/ctsm/netcdf_utils.py
@@ -0,0 +1,141 @@
+"""
+Helper functions for working with netCDF files
+"""
+
+import numpy as np
+import xarray as xr
+from netCDF4 import Dataset  # pylint: disable=no-name-in-module
+
+
+def _is_dtype_nan_capable(ndarray: np.ndarray):
+    """
+    Given a numpy array, return True if it's capable of taking a NaN
+    """
+    try:
+        np.isnan(ndarray)
+        return True
+    except TypeError:
+        return False
+
+
+def _are_dicts_identical_nansequal(dict0: dict, dict1: dict, keys_to_ignore=None):
+    """
+    Compare two dictionaries, considering NaNs to be equal. Don't be strict here about types; if
+    they can be coerced to comparable types and then they match, return True.
+    """
+    # pylint: disable=too-many-return-statements
+
+    if keys_to_ignore is None:
+        keys_to_ignore = []
+    keys_to_ignore = np.array(keys_to_ignore)
+
+    if len(dict0) != len(dict1):
+        return False
+    for key, value0 in dict0.items():
+        if key in keys_to_ignore:
+            continue
+        if key not in dict1:
+            return False
+        value1 = dict1[key]
+
+        # Coerce to numpy arrays to simplify comparison code
+        value0 = np.array(value0)
+        value1 = np.array(value1)
+
+        # Compare, only asking to check equal NaNs if both are capable of taking NaN values
+        both_are_nan_capable = _is_dtype_nan_capable(value0) and _is_dtype_nan_capable(value1)
+        if not np.array_equal(value0, value1, equal_nan=both_are_nan_capable):
+            return False
+
+    return True
+
+
+def get_netcdf_format(file_path):
+    """
+    Get format of netCDF file
+    """
+    with Dataset(file_path, "r") as netcdf_file:
+        netcdf_format = netcdf_file.data_model
+    return netcdf_format
+
+
+def _is_dataarray_metadata_identical(da0: xr.DataArray, da1: xr.DataArray, keys_to_ignore=None):
+    """
+    Check whether two DataArrays have identical-enough metadata
+    """
+
+    # Check data type
+    if da0.dtype != da1.dtype:
+        return False
+
+    # Check encoding
+    if not _are_dicts_identical_nansequal(
+        da0.encoding, da1.encoding, keys_to_ignore=keys_to_ignore
+    ):
+        return False
+
+    # Check attributes
+    if not _are_dicts_identical_nansequal(da0.attrs, da1.attrs):
+        return False
+
+    # Check name
+    if da0.name != da1.name:
+        return False
+
+    # Check dims
+    if da0.dims != da1.dims:
+        return False
+
+    return True
+
+
+def _is_dataarray_data_identical(da0: xr.DataArray, da1: xr.DataArray):
+    """
+    Check whether two DataArrays have identical data
+    """
+    # pylint: disable=too-many-return-statements
+
+    # Check sizes
+    if da0.sizes != da1.sizes:
+        return False
+
+    # Check coordinates
+    if bool(da0.coords) or bool(da1.coords):
+        if not bool(da0.coords) or not bool(da1.coords):
+            return False
+        if not da0.coords.equals(da1.coords):
+            return False
+
+    # Check values ("The array's data converted to numpy.ndarray")
+    if not np.array_equal(da0.values, da1.values):
+        # Try-except to avoid TypeError from putting NaN-incapable dtypes through
+        # np.array_equal(..., equal_nan=True)
+        try:
+            if not np.array_equal(da0.values, da1.values, equal_nan=True):
+                return False
+        except TypeError:
+            return False
+
+    # Check data ("The DataArray's data as an array. The underlying array type (e.g. dask, sparse,
+    # pint) is preserved.")
+    da0_data_type = type(da0.data)
+    if not isinstance(da1.data, da0_data_type):
+        return False
+    if not isinstance(da0.data, np.ndarray):
+        raise NotImplementedError(f"Add support for comparing two objects of type {da0_data_type}")
+
+    return True
+
+
+def are_xr_dataarrays_identical(da0: xr.DataArray, da1: xr.DataArray, keys_to_ignore=None):
+    """
+    Comprehensively check whether two DataArrays are identical
+    """
+    if not _is_dataarray_metadata_identical(da0, da1, keys_to_ignore=keys_to_ignore):
+        return False
+
+    if not _is_dataarray_data_identical(da0, da1):
+        return False
+
+    # Fallback to however xarray defines equality, in case we missed something above
+    return da0.equals(da1)
diff --git a/python/ctsm/param_utils/__init__.py b/python/ctsm/param_utils/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/python/ctsm/param_utils/paramfile_shared.py b/python/ctsm/param_utils/paramfile_shared.py
new file mode 100644
index 0000000000..b9860be0db
--- /dev/null
+++ b/python/ctsm/param_utils/paramfile_shared.py
@@ -0,0 +1,134 @@
+"""
+Functions etc. shared among parameter file utilities
+"""
+
+import argparse
+import xarray as xr
+
+from ctsm.netcdf_utils import are_xr_dataarrays_identical
+
+PFTNAME_VAR = "pftname"
+
+
+def are_paramfile_dataarrays_identical(da0: xr.DataArray, da1: xr.DataArray):
+    """
+    Check whether parameter DataArrays are identical enough, ignoring some metadata
+    """
+    return are_xr_dataarrays_identical(da0, da1, keys_to_ignore=["source", "original_shape"])
+
+
+def check_pfts_in_paramfile(selected_pfts, ds):
+    """
+    Check that the given PFTs are present in the parameter file.
+
+    Parameters
+    ----------
+    selected_pfts : list of str
+        List of PFT names to check.
+    ds : xarray.Dataset
+        The parameter file dataset.
+
+    Returns
+    -------
+    list of str
+        List of all PFT names in the file.
+
+    Raises
+    ------
+    KeyError
+        If any selected PFT is not found in the file, or if PFTNAME_VAR is missing.
+    """
+    if PFTNAME_VAR not in ds:
+        raise KeyError(f"paramfile missing variable: {PFTNAME_VAR}")
+    pft_names = get_pft_names(ds)
+    pfts_not_in_file = []
+    for pft in selected_pfts:
+        if pft not in pft_names:
+            pfts_not_in_file += [pft]
+    if pfts_not_in_file:
+        raise KeyError(f"PFT(s) not found in parameter file: {', '.join(pfts_not_in_file)}")
+
+    return pft_names
+
+
+def get_pft_names(ds):
+    """
+    Get the list of PFT names from the parameter file dataset.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        The parameter file dataset.
+
+    Returns
+    -------
+    list of str
+        List of PFT names.
+    """
+    pft_names = [pft.decode().strip() for pft in ds[PFTNAME_VAR].values]
+    return pft_names
+
+
+def get_selected_pft_indices(selected_pfts, pft_names):
+    """
+    Get indices of selected PFTs in the list of all PFT names.
+
+    Parameters
+    ----------
+    selected_pfts : list of str
+        List of PFT names to select.
+    pft_names : list of str
+        List of all PFT names.
+
+    Returns
+    -------
+    list of int
+        Indices of selected PFTs.
+    """
+    indices = [i for i, name in enumerate(pft_names) if name in selected_pfts]
+    return indices
+
+
+def open_paramfile(file_in, mask_and_scale=False):
+    """
+    Open a parameter file as an xarray.Dataset.
+
+    Parameters
+    ----------
+    file_in : str
+        Path to the input netCDF file.
+    mask_and_scale : bool, optional
+        Whether to apply mask and scale (default: False).
+
+    Returns
+    -------
+    xarray.Dataset
+        The opened dataset.
+    """
+    return xr.open_dataset(file_in, decode_timedelta=False, mask_and_scale=mask_and_scale)
+
+
+def paramfile_parser_setup(description):
+    """
+    Set up an argument parser for parameter file utilities.
+
+    Parameters
+    ----------
+    description : str
+        Description for the argument parser.
+
+    Returns
+    -------
+    tuple
+        (parser, pft_flags) where parser is an ArgumentParser and pft_flags is a list of flags for
+        PFT argument.
+    """
+    parser = argparse.ArgumentParser(
+        description=description, formatter_class=argparse.RawTextHelpFormatter
+    )
+    parser.add_argument("-i", "--input", required=True, help="Input netCDF file")
+
+    # Flags that can be used for the PFT argument
+    pft_flags = ["-p", "--pft"]
+
+    return parser, pft_flags
diff --git a/python/ctsm/param_utils/query_paramfile.py b/python/ctsm/param_utils/query_paramfile.py
new file mode 100644
index 0000000000..6103010827
--- /dev/null
+++ b/python/ctsm/param_utils/query_paramfile.py
@@ -0,0 +1,110 @@
+"""
+Query parameters in a CTSM paramfile.
+
+This script allows users to print the values of one or more parameters from a CTSM parameter file
+(netCDF format). Users can specify variables and optionally filter output by PFTs. The script
+handles multi-dimensional variables and provides formatted output for PFT-specific parameters.
+"""
+
+from ctsm.args_utils import comma_separated_list
+from ctsm.param_utils.paramfile_shared import paramfile_parser_setup
+from ctsm.param_utils.paramfile_shared import open_paramfile
+from ctsm.param_utils.paramfile_shared import PFTNAME_VAR, check_pfts_in_paramfile
+from ctsm.param_utils.paramfile_shared import get_selected_pft_indices, get_pft_names
+
+
+def get_arguments():
+    """
+    Parse command-line arguments for querying variables from a netCDF file.
+
+    Returns
+    -------
+    argparse.Namespace
+        Parsed arguments with attributes:
+            - input: Path to the netCDF file
+            - variables: List of variable names to query
+            - pft: Optional list of PFT names to print
+    """
+    parser, pft_flags = paramfile_parser_setup(
+        "Print values of one or more parameters from a CTSM paramfile."
+    )
+    parser.add_argument(
+        "variables",
+        help="Names of variables to query (space-separated)",
+        nargs="*",
+    )
+    parser.add_argument(
+        *pft_flags,
+        help="Comma-separated list of PFT names to print (only applies to PFT-specific variables)",
+        type=comma_separated_list,
+    )
+    args = parser.parse_args()
+    return args
+
+
+def print_values(ds, var, selected_pfts, pft_names):
+    """
+    Print the values of a variable from the dataset, optionally filtered by PFTs.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        The opened netCDF dataset.
+    var : str
+        Variable name to print.
+    selected_pfts : list or None
+        List of selected PFT names to print, or None to print all.
+    pft_names : list or None
+        List of all PFT names in the file.
+    """
+    data = ds[var].values
+    if list(ds[var].dims) == ["pft"]:
+        print(var + ":")
+        indices = range(len(pft_names))
+        if selected_pfts is not None:
+            indices = get_selected_pft_indices(selected_pfts, pft_names)
+            max_name_len = max(len(pft_names[i]) for i in indices) if indices else 0
+        else:
+            max_name_len = max(len(name) for name in pft_names)
+        for p in indices:
+            print(f"   {pft_names[p]:<{max_name_len}}: {data[p]}")
+    elif ds[var].ndim > 1:
+        print(f"{var}:")
+        print(data)
+    else:
+        print(f"{var}: {data}")
+
+
+def main():
+    """
+    Main entry point for query_paramfile.
+
+    Parses arguments, opens the netCDF file, and prints requested variable values,
+    optionally filtered by PFTs.
+    """
+    args = get_arguments()
+
+    ds = open_paramfile(args.input, mask_and_scale=True)
+
+    # If user didn't specify variables, print all
+    if not args.variables:
+        args.variables = ds.variables
+
+    selected_pfts = args.pft
+    pft_names = None
+    if selected_pfts:
+        pft_names = check_pfts_in_paramfile(selected_pfts, ds)
+    elif PFTNAME_VAR in ds.coords:
+        pft_names = get_pft_names(ds)
+
+    for var in args.variables:
+        if var in ds.variables:
+            print_values(ds, var, selected_pfts, pft_names)
+        else:
+            print(f"Variable '{var}' not found in {args.input}")
+
+    ds.close()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/python/ctsm/param_utils/set_paramfile.py b/python/ctsm/param_utils/set_paramfile.py
new file mode 100644
index 0000000000..2ceb52069b
--- /dev/null
+++ b/python/ctsm/param_utils/set_paramfile.py
@@ -0,0 +1,381 @@
+"""
+Tool for changing parameters on CTSM paramfile.
+
+This script allows users to modify one or more parameters in a CTSM parameter file (netCDF format).
+It supports selecting specific PFTs, dropping other PFTs, and changing parameter values, including
+setting values to fill (missing) values. The script ensures safe file handling and provides
+detailed error checking for argument validity and parameter changes.
+"""
+
+import os
+import sys
+from datetime import datetime
+import numpy as np
+import xarray as xr
+
+from ctsm.args_utils import comma_separated_list
+from ctsm.netcdf_utils import get_netcdf_format
+from ctsm.param_utils.paramfile_shared import paramfile_parser_setup, open_paramfile
+from ctsm.param_utils.paramfile_shared import check_pfts_in_paramfile, get_selected_pft_indices
+from ctsm.param_utils.paramfile_shared import PFTNAME_VAR
+
+
+def check_arguments(args):
+    """
+    Validate command-line arguments for set_paramfile.
+
+    Checks for existence of input file, prevents overwriting output files,
+    and ensures logical consistency of PFT-related options.
+    """
+    if not os.path.exists(args.input):
+        raise FileNotFoundError(args.input)
+
+    # Avoid potentially overwriting canonical files
+    if os.path.exists(args.output):
+        raise FileExistsError(args.output)
+
+    # --drop-other-pfts makes no sense without --pfts
+    if args.drop_other_pfts and not args.pft:
+        raise RuntimeError("--drop-other-pfts makes no sense without -p/--pft")
+
+
+def get_arguments():
+    """
+    Parse command-line arguments for setting variables on a netCDF file.
+
+    Returns
+    -------
+    argparse.Namespace
+        Parsed arguments with attributes:
+            - input: Path to the input netCDF file
+            - output: Path to the output netCDF file
+            - pft: Optional list of PFT names whose values you want to change
+            - drop_other_pfts: Boolean flag to drop PFTs not specified
+            - param_changes: List of parameter changes to apply
+    """
+    parser, pft_flags = paramfile_parser_setup(
+        "Change values of one or more parameters in a CTSM paramfile."
+    )
+
+    parser.add_argument(
+        "-o", "--output", required=True, help="Output netCDF file. Must not already exist."
+    )
+
+    # TODO: Add mutually-exclusive --exclude-pfts argument for PFTs you DON'T want to include
+    parser.add_argument(
+        *pft_flags,
+        help="Comma-separated list of PFTs to include (only applies to PFT-specific variables)",
+        type=comma_separated_list,
+    )
+
+    parser.add_argument(
+        "--drop-other-pfts",
+        help=f"Do not include PFTs other than the ones given in {'/'.join(pft_flags)}",
+        action="store_true",
+    )
+
+    parser.add_argument(
+        "param_changes",
+        help=(
+            "Parameter changes to apply. Use nan to set to the fill value. E.g.:\n"
+            "   param1=new_value1 pftparam=pft1_val,nan,... param3=nan"
+        ),
+        nargs="*",
+    )
+
+    args = parser.parse_args()
+    check_arguments(args)
+
+    return args
+
+
+def is_integer(obj):
+    """
+    Determine if an object is an integer or a numpy array of integer dtype.
+
+    Parameters
+    ----------
+    obj : object
+        Object to check.
+
+    Returns
+    -------
+    bool
+        True if obj is an integer or numpy array of integer dtype, False otherwise.
+    """
+    if isinstance(obj, np.ndarray):
+        obj_type = obj.dtype
+    else:
+        obj_type = type(obj)
+    return np.issubdtype(obj_type, np.integer)
+
+
+def check_correct_ndims(da, new_value, throw_error=False):
+    """
+    Check that the new value for a parameter has the correct number of dimensions.
+
+    Parameters
+    ----------
+    da : xarray.DataArray
+        The parameter DataArray to check against.
+    new_value : array-like
+        The new value to assign.
+    throw_error : bool, optional
+        If True, raise an error on mismatch.
+
+    Returns
+    -------
+    bool
+        True if dimensions match, False otherwise.
+    """
+    expected = da.ndim
+    actual = np.array(new_value).ndim
+    is_ndim_correct = expected == actual
+    if throw_error and not is_ndim_correct:
+        raise RuntimeError(f"Incorrect N dims: Expected {expected}, got {actual}")
+    return is_ndim_correct
+
+
+def drop_other_pfts(selected_pfts, ds):
+    """
+    Drop PFTs from the dataset that are not in the selected list.
+
+    Parameters
+    ----------
+    selected_pfts : list of str
+        List of PFT names to retain.
+    ds : xarray.Dataset
+        The parameter file dataset.
+
+    Returns
+    -------
+    xarray.Dataset
+        Dataset containing only the selected PFTs.
+    """
+    pft_names = check_pfts_in_paramfile(selected_pfts, ds)
+    indices = get_selected_pft_indices(selected_pfts, pft_names)
+    ds = ds.isel({"pft": indices})
+    return ds
+
+
+def _add_cmd_to_history(ds):
+    """
+    Prepend the calling command and timestamp to the netCDF history attribute.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        Dataset to update.
+
+    Returns
+    -------
+    xarray.Dataset
+        Dataset with updated history attribute.
+    """
+    if "history" not in ds.attrs:
+        ds.attrs["history"] = ""
+    cmd_items = [f"'{x}'" if " " in x else x for x in sys.argv]
+    datetime_str = datetime.now().strftime("%Y-%m-%d %H:%M:%S")
+    ds.attrs["history"] = f"{datetime_str}: {' '.join(cmd_items)}\n{ds.attrs['history']}"
+    return ds
+
+
+def save_paramfile(ds_out: xr.Dataset, output_path, *, nc_format="NETCDF3_CLASSIC"):
+    """
+    Save an xarray Dataset to a netCDF parameter file.
+
+    Parameters
+    ----------
+    ds_out : xarray.Dataset
+        Dataset to save.
+    output_path : str
+        Path to output netCDF file.
+    nc_format : str, optional
+        NetCDF format to use (default: "NETCDF3_CLASSIC").
+    """
+
+    # We don't want to add _FillValue to parameters that didn't already have one. This dict will
+    # track such parameters and be passed to .to_netcdf(..., encoding=encoding)
+    encoding = {}
+    for var in ds_out:
+        if "_FillValue" not in ds_out[var].encoding:
+            encoding[var] = {"_FillValue": None}
+
+    ds_out = _add_cmd_to_history(ds_out)
+
+    ds_out.to_netcdf(output_path, format=nc_format, encoding=encoding)
+
+
+def _replace_nans_with_fill(var_encoding, new_value, *, chg=None):
+    """
+    Replace NaNs in the new parameter value array with the fill value.
+
+    Parameters
+    ----------
+    var_encoding : dict
+        Encoding of an xarray DataArray.
+    new_value : numpy.ndarray
+        Array of new values.
+    chg : str, optional (keyword-only)
+        Change string from command line, used for error messages.
+
+    Returns
+    -------
+    numpy.ndarray
+        Array with NaNs replaced by fill value.
+    """
+    if any(np.isnan(np.atleast_1d(new_value))):
+        # TODO: Add code to add fill value to parameters without it
+        if "_FillValue" not in var_encoding:
+            msg = "Can't set parameter to fill value if it doesn't already have one"
+            if chg is not None:
+                msg += f": {chg}"
+            raise NotImplementedError(msg)
+        fill_value = var_encoding["_FillValue"]
+        new_value[np.isnan(new_value)] = fill_value
+
+    return new_value
+
+
+def _convert_to_output_dtype(ds_out, var, new_value, *, chg=None):
+    """
+    Convert new_value array to the output variable's data type.
+
+    Parameters
+    ----------
+    ds_out : xarray.Dataset
+        Output dataset.
+    var : str
+        Variable name.
+    new_value : numpy.ndarray
+        Array of new values.
+    chg : str, optional (keyword-only)
+        Change string from command line, used for error messages.
+
+    Returns
+    -------
+    numpy.ndarray
+        Array converted to output dtype.
+    """
+    try:
+        new_value = new_value.astype(type(ds_out[var].dtype))
+    except ValueError as e:
+        msg = str(e)
+        if "invalid literal for int() with base 10" in msg:
+            # Throw a nicer error message including the entire requested change
+            err_msg = "Invalid assignment to an integer parameter"
+            if chg is not None:
+                err_msg += f": {chg}"
+            raise ValueError(err_msg) from e
+        raise e
+
+    return new_value
+
+
+def apply_new_value_to_parameter(args, ds_out, var, new_value, var_encoding, *, chg=None):
+    """
+    Apply a new value to a parameter in the output dataset, handling PFT selection, dimension
+    checks, fill value replacement, and assignment.
+
+    Parameters
+    ----------
+    args : argparse.Namespace
+        Parsed command-line arguments.
+    ds_out : xarray.Dataset
+        Output dataset to modify.
+    var : str
+        Name of the variable to change.
+    new_value : numpy.ndarray
+        Array of new values to assign.
+    var_encoding : dict
+        Encoding dictionary for the variable, used for fill value replacement.
+    chg : str, optional (keyword-only)
+        Change string from command line, used for error messages.
+
+    Returns
+    -------
+    xarray.Dataset
+        Modified output dataset with the new parameter value applied.
+    """
+    # Are we acting on just some PFTs? If so, we'll need some stuff.
+    just_some_pfts = PFTNAME_VAR in ds_out[var].coords and args.pft
+    # pylint is probably wrong with the possibly-used-before-assignment warning, but do this
+    # here just to placate it. Make it an invalid index so we get an error if we try to use
+    # it.
+    indices = -1
+    if just_some_pfts:
+        pft_names = check_pfts_in_paramfile(args.pft, ds_out)
+        indices = get_selected_pft_indices(args.pft, pft_names)
+
+    # Check that correct number of dimensions were given for new values. Special handling needed
+    # if we're just acting on one PFT.
+    da_to_check = ds_out[var]
+    if just_some_pfts and len(args.pft) == 1:
+        da_to_check = da_to_check.isel(pft=indices).squeeze()
+    check_correct_ndims(da_to_check, new_value, throw_error=True)
+
+    # Handle the situation where we're only changing values for some PFTs but keeping the others
+    if just_some_pfts and not args.drop_other_pfts:
+        tmp = ds_out[var].values.copy()
+        tmp[indices] = new_value
+        new_value = tmp
+
+    # Ensure that any NaNs are replaced with the fill value
+    new_value = _replace_nans_with_fill(var_encoding, new_value, chg=chg)
+
+    # This can be needed if, e.g., you're selecting and changing just one PFT
+    if ds_out[var].values.ndim > 0 and new_value.ndim == 0:
+        new_value = np.atleast_1d(new_value)
+
+    ds_out[var].values = new_value
+    return ds_out
+
+
+def main():
+    """
+    Main entry point for set_paramfile.
+
+    Parses arguments, opens the input netCDF file, applies requested changes,
+    and saves the modified dataset to a new netCDF file.
+    """
+    args = get_arguments()
+
+    ds_in_masked_scaled = open_paramfile(args.input, mask_and_scale=True)
+
+    ds_in = open_paramfile(args.input)
+    ds_out = ds_in.copy()
+
+    # If --drop-other-pfts was given, drop PFTs not in args.pft
+    if args.drop_other_pfts:
+        ds_out = drop_other_pfts(args.pft, ds_out)
+
+    # Apply parameter changes, if any
+    for chg in args.param_changes:
+        var, new_value = chg.split("=")
+
+        # TODO: Add handling of multi-dimensional parameters
+        if ds_out[var].ndim > 1:
+            raise NotImplementedError("Can't yet change multi-dimensional parameters")
+
+        # Split at commas, if any, and convert to numpy array
+        new_value = np.array(new_value.split(",")).squeeze()
+
+        # TODO: Add code to set integer variables to their missing value. This is harder than it
+        # sounds.
+        if np.any(np.char.lower(new_value) == "nan") and is_integer(ds_out[var].values):
+            raise NotImplementedError(f"Can't set integer parameter to fill value: {chg}")
+
+        # Convert to the output data type
+        new_value = _convert_to_output_dtype(ds_out, var, new_value, chg=chg)
+
+        # Extract some information
+        var_encoding = ds_in_masked_scaled[var].encoding
+
+        # Apply new value to parameter
+        ds_out = apply_new_value_to_parameter(args, ds_out, var, new_value, var_encoding, chg=chg)
+
+    save_paramfile(ds_out, args.output, nc_format=get_netcdf_format(args.input))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/python/ctsm/run_sys_tests.py b/python/ctsm/run_sys_tests.py
index 8108668246..88fad06169 100644
--- a/python/ctsm/run_sys_tests.py
+++ b/python/ctsm/run_sys_tests.py
@@ -749,6 +749,14 @@ def _check_py_env(test_attributes):
         except ModuleNotFoundError as err:
             raise ModuleNotFoundError("modify_fsurdat" + err_msg) from err
 
+    # Check requirements for using set_paramfile Python module, if needed
+    set_paramfile_users = ["SETPARAMFILE"]
+    if any(any(u in t for u in set_paramfile_users) for t in test_attributes):
+        try:
+            import ctsm.param_utils.set_paramfile
+        except ModuleNotFoundError as err:
+            raise ModuleNotFoundError("set_paramfile" + err_msg) from err
+
     # Check requirements for RXCROPMATURITY, if needed
     if any("RXCROPMATURITY" in t for t in test_attributes):
         try:
diff --git a/python/ctsm/test/test_sys_query_paramfile.py b/python/ctsm/test/test_sys_query_paramfile.py
new file mode 100755
index 0000000000..c01e7bc5e8
--- /dev/null
+++ b/python/ctsm/test/test_sys_query_paramfile.py
@@ -0,0 +1,147 @@
+#!/usr/bin/env python3
+
+"""System tests for query_paramfile"""
+
+import unittest
+import os
+import sys
+import io
+from contextlib import redirect_stdout
+import tempfile
+import shutil
+import xarray as xr
+
+from ctsm import unit_testing
+
+from ctsm.param_utils import query_paramfile as qp
+
+# Allow names that pylint doesn't like, because otherwise I find it hard
+# to make readable unit test names
+# pylint: disable=invalid-name
+
+PARAMFILE = os.path.join(
+    os.path.dirname(__file__), "testinputs", "ctsm5.3.041.Nfix_params.v13.c250221_upplim250.nc"
+)
+
+
+class TestSysQueryParamfile(unittest.TestCase):
+    """System tests of query_paramfile"""
+
+    def setUp(self):
+        self.orig_argv = sys.argv
+        self.tempdir = tempfile.mkdtemp()
+
+    def tearDown(self):
+        sys.argv = self.orig_argv
+        shutil.rmtree(self.tempdir, ignore_errors=True)
+
+    def test_query_paramfile_scalar_nopfts(self):
+        """Test that print_values works with scalar parameter and no PFTs specified"""
+
+        sys.argv = ["get_arguments", "-i", PARAMFILE, "phenology_soil_depth"]
+
+        f = io.StringIO()
+        with redirect_stdout(f):
+            qp.main()
+        out = f.getvalue()
+        self.assertEqual("phenology_soil_depth: 0.08\n", out)
+
+    def test_query_paramfile_scalar_ignorepfts(self):
+        """Test that print_values works with scalar parameter and PFTs specified (ignored)"""
+
+        sys.argv = ["get_arguments", "-i", PARAMFILE, "phenology_soil_depth", "--pft", "c3_crop"]
+
+        f = io.StringIO()
+        with redirect_stdout(f):
+            qp.main()
+        out = f.getvalue()
+        self.assertEqual("phenology_soil_depth: 0.08\n", out)
+
+    def test_query_paramfile_pft_selectnone(self):
+        """Test that print_values works with PFT-dim parameter and no PFTs specified"""
+
+        sys.argv = ["get_arguments", "-i", PARAMFILE, "rswf_min"]
+
+        f = io.StringIO()
+        with redirect_stdout(f):
+            qp.main()
+        out = f.getvalue()
+        self.assertRegex(
+            out,
+            (
+                r"rswf_min:\n"
+                r"\s+not_vegetated\s*: 0\.25\n"
+                r"\s+needleleaf_evergreen_temperate_tree\s*: 0\.25\n.*"
+            ),
+        )
+
+    def test_query_paramfile_pft_select2(self):
+        """Test that print_values works with PFT-dim parameter and two PFTs specified"""
+
+        sys.argv = [
+            "get_arguments",
+            "-i",
+            PARAMFILE,
+            "--pft",
+            "not_vegetated,needleleaf_evergreen_temperate_tree",
+            "rswf_min",
+        ]
+
+        f = io.StringIO()
+        with redirect_stdout(f):
+            qp.main()
+        out = f.getvalue()
+        self.assertRegex(
+            out,
+            (
+                r"rswf_min:\n"
+                r"\s+not_vegetated\s*: 0\.25\n"
+                r"\s+needleleaf_evergreen_temperate_tree\s*: 0\.25\n"
+            ),
+        )
+
+    def test_query_paramfile_no_variables_fake(self):
+        """
+        Test that print_values prints every variable when no variables are given. Use a small fake
+        paramfile so we can check that what gets printed is what we expect.
+        """
+
+        fake_da1 = xr.DataArray(data=[1, 2, 3])
+        fake_da2 = xr.DataArray(data=[4, 5, 6])
+        fake_ds = xr.Dataset(data_vars={"fake1": fake_da1, "fake2": fake_da2})
+        fake_nc_path = os.path.join(self.tempdir, "fake_da.nc")
+        fake_ds.to_netcdf(fake_nc_path)
+
+        sys.argv = [
+            "query_paramfile",
+            "-i",
+            fake_nc_path,
+        ]
+
+        f = io.StringIO()
+        with redirect_stdout(f):
+            qp.main()
+        out = f.getvalue()
+        self.assertRegex(out, (r"fake1: \[1 2 3\]\nfake2: \[4 5 6\]\n"))
+
+    def test_query_paramfile_no_variables_real(self):
+        """
+        Test that query_paramfile doesn't error when trying to print every variable from a real
+        paramfile. Don't actually check that it matches what we expect; that's done in
+        test_query_paramfile_no_variables_fake.
+        """
+
+        sys.argv = [
+            "query_paramfile",
+            "-i",
+            PARAMFILE,
+        ]
+
+        f = io.StringIO()
+        with redirect_stdout(f):
+            qp.main()
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_sys_set_paramfile.py b/python/ctsm/test/test_sys_set_paramfile.py
new file mode 100755
index 0000000000..7d82a32302
--- /dev/null
+++ b/python/ctsm/test/test_sys_set_paramfile.py
@@ -0,0 +1,738 @@
+#!/usr/bin/env python3
+
+"""System tests for set_paramfile"""
+
+import unittest
+import os
+import sys
+import shutil
+import tempfile
+import numpy as np
+import xarray as xr
+
+from ctsm import unit_testing
+
+from ctsm.netcdf_utils import get_netcdf_format
+from ctsm.param_utils import set_paramfile as sp
+from ctsm.param_utils.paramfile_shared import open_paramfile, are_paramfile_dataarrays_identical
+from ctsm.param_utils.paramfile_shared import check_pfts_in_paramfile, get_selected_pft_indices
+
+# Allow names that pylint doesn't like, because otherwise I find it hard
+# to make readable unit test names
+# pylint: disable=invalid-name
+
+
+PARAMFILE = os.path.join(
+    os.path.dirname(__file__), "testinputs", "ctsm5.3.041.Nfix_params.v13.c250221_upplim250.nc"
+)
+
+
+class TestSysSetParamfile(unittest.TestCase):
+    """System tests of set_paramfile"""
+
+    # pylint: disable=too-many-public-methods
+
+    def setUp(self):
+        self.orig_argv = sys.argv
+        self.tempdir = tempfile.mkdtemp()
+
+    def tearDown(self):
+        sys.argv = self.orig_argv
+        shutil.rmtree(self.tempdir, ignore_errors=True)
+
+    def test_set_paramfile_copyfile(self):
+        """Test that set_paramfile can straight-up copy to a new file"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        sys.argv = ["set_paramfile", "-i", PARAMFILE, "-o", output_path]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+        ds_in = open_paramfile(PARAMFILE)
+        ds_out = open_paramfile(output_path)
+
+        # Check that contents are functionally identical
+        self.assertEqual(ds_in, ds_out)
+
+        # Check that both are the same kind of netCDF
+        self.assertEqual(get_netcdf_format(PARAMFILE), get_netcdf_format(output_path))
+
+    def test_set_paramfile_copy_without_adding_fillvalue(self):
+        """Test that set_paramfile can copy to a new file without adding _FillValue"""
+        input_path = os.path.join(self.tempdir, "input.nc")
+        output_path = os.path.join(self.tempdir, "output.nc")
+        param_name = "param0"
+
+        # Save a test paramfile without _FillValue
+        da = xr.DataArray(data=np.float64([1, 2, 3]))
+        da.encoding["_FillValue"] = None
+        ds = xr.Dataset(data_vars={param_name: da})
+        ds.to_netcdf(input_path, encoding={param_name: {"_FillValue": None}})
+        ds_in = open_paramfile(input_path)
+        self.assertFalse("_FillValue" in ds_in[param_name].encoding)
+        self.assertFalse("_FillValue" in ds_in[param_name].attrs)
+
+        # Use set_paramfile to copy to new file
+        sys.argv = ["set_paramfile", "-i", input_path, "-o", output_path]
+        sp.main()
+
+        # Check that _FillValue wasn't added
+        ds_out = open_paramfile(output_path)
+        self.assertFalse("_FillValue" in ds_out[param_name].encoding)
+        self.assertFalse("_FillValue" in ds_out[param_name].attrs)
+
+    def test_set_paramfile_extractpfts(self):
+        """Test that set_paramfile can copy to a new file with only some requested PFTs"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        pfts_to_include = ["not_vegetated", "needleleaf_evergreen_temperate_tree"]
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "-p",
+            ",".join(pfts_to_include),
+            "--drop-other-pfts",
+        ]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+        ds_in = open_paramfile(PARAMFILE)
+        ds_out = open_paramfile(output_path)
+
+        # Check that included variables/coords match
+        for var in ds_in.variables:
+            actual = ds_out[var]
+            if sp.PFTNAME_VAR in ds_in[var].coords:
+                expected = ds_in[var].isel(pft=[0, 1])
+            else:
+                expected = ds_in[var]
+            self.assertTrue(are_paramfile_dataarrays_identical(expected, actual))
+
+    def test_set_paramfile_changeparams_scalar_errors_given_list(self):
+        """Test that set_paramfile errors if given a list for a scalar parameter"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "a_coef=0.87,0.91",
+        ]
+        with self.assertRaisesRegex(RuntimeError, "Incorrect N dims"):
+            sp.main()
+
+    def test_set_paramfile_changeparam_1d_errors_given_scalar(self):
+        """Test that set_paramfile errors if given a scalar for a 1-d parameter"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "xl=0.724",
+        ]
+        with self.assertRaisesRegex(RuntimeError, "Incorrect N dims"):
+            sp.main()
+
+    def test_set_paramfile_changeparams_scalar_double(self):
+        """Test that set_paramfile can copy to a new file with some scalar double params changed"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "a_coef=0.87",
+            "bgc_cn_s2=87",
+        ]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+        ds_in = open_paramfile(PARAMFILE)
+        ds_out = open_paramfile(output_path)
+
+        for var in ds_in.variables:
+            # Check that all variables/coords are equal except the ones we changed, which should be
+            # set to what we asked
+            if var == "a_coef":
+                self.assertTrue(ds_in[var].values == 0.13)
+                self.assertTrue(ds_out[var].values == 0.87)
+            elif var == "bgc_cn_s2":
+                self.assertTrue(ds_in[var].values == 11)
+                self.assertTrue(ds_out[var].values == 87)
+            else:
+                self.assertTrue(are_paramfile_dataarrays_identical(ds_in[var], ds_out[var]))
+
+            # Check that data type hasn't changed
+            self.assertTrue(ds_in[var].dtype == ds_out[var].dtype)
+
+            # Check that fill value hasn't changed
+            if "_FillValue" in ds_in[var].encoding:
+                fv_in = ds_in[var].encoding["_FillValue"]
+                fv_out = ds_out[var].encoding["_FillValue"]
+                if isinstance(fv_in, bytes):
+                    self.assertTrue(isinstance(fv_out, bytes))
+                    self.assertEqual(fv_in, fv_out)
+                else:
+                    self.assertEqual(np.isnan(fv_in), np.isnan(fv_out))
+                    if not np.isnan(fv_in):
+                        self.assertEqual(fv_in, fv_out)
+
+    def test_set_paramfile_changeparams_1d_double(self):
+        """
+        Test that set_paramfile can copy to a new file with a 1-d double param changed (not PFT-
+        dimensioned)
+        """
+        output_path = os.path.join(self.tempdir, "output.nc")
+        this_var = "mimics_fmet"
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            f"{this_var}=0.1,0.2,0.3,0.4",
+        ]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+        ds_in = open_paramfile(PARAMFILE)
+        ds_out = open_paramfile(output_path)
+
+        for var in ds_in.variables:
+            # Check that all variables/coords are equal except the ones we changed, which should be
+            # set to what we asked
+            if var == this_var:
+                self.assertTrue(
+                    np.array_equal(ds_in[var].values, np.array([0.75, 0.85, 0.013, 40]))
+                )
+                self.assertTrue(np.array_equal(ds_out[var].values, np.array([0.1, 0.2, 0.3, 0.4])))
+            else:
+                self.assertTrue(are_paramfile_dataarrays_identical(ds_in[var], ds_out[var]))
+
+            # Check that data type hasn't changed
+            self.assertTrue(ds_in[var].dtype == ds_out[var].dtype)
+
+            # Check that fill value hasn't changed
+            if "_FillValue" in ds_in[var].encoding:
+                fv_in = ds_in[var].encoding["_FillValue"]
+                fv_out = ds_out[var].encoding["_FillValue"]
+                if isinstance(fv_in, bytes):
+                    self.assertTrue(isinstance(fv_out, bytes))
+                    self.assertEqual(fv_in, fv_out)
+                else:
+                    self.assertEqual(np.isnan(fv_in), np.isnan(fv_out))
+                    if not np.isnan(fv_in):
+                        self.assertEqual(fv_in, fv_out)
+
+    def test_set_paramfile_changeparams_scalar_int(self):
+        """Test that set_paramfile can copy to a new file with a scalar integer param changed"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        this_var = "upplim_destruct_metamorph"
+        new_value = 1987
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            f"upplim_destruct_metamorph={new_value}",
+        ]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+        ds_in = open_paramfile(PARAMFILE)
+        ds_out = open_paramfile(output_path)
+
+        # Check that the variable in question is actually an integer to begin with
+        self.assertTrue(sp.is_integer(ds_in[this_var].values))
+        # Also check that it actually differs from our new value
+        self.assertTrue(ds_in[this_var].values != new_value)
+
+        for var in ds_in.variables:
+            # Check that all variables/coords are equal except the one we changed, which should be
+            # set to what we asked
+            if var == this_var:
+                self.assertTrue(ds_out[var].values == new_value)
+            else:
+                self.assertTrue(ds_in[var].equals(ds_out[var]))
+
+            # Check that data type hasn't changed
+            self.assertTrue(ds_in[var].dtype == ds_out[var].dtype)
+
+    def test_set_paramfile_extractpfts_changeparam_dbl(self):
+        """
+        Test that set_paramfile can (1) copy to a new file with only some requested PFTs and (2)
+        change the values of double parameters of those PFTs
+        """
+        output_path = os.path.join(self.tempdir, "output.nc")
+        pfts_to_include = ["not_vegetated", "needleleaf_evergreen_temperate_tree"]
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "-p",
+            ",".join(pfts_to_include),
+            "--drop-other-pfts",
+            "xl=0.724,0.87",
+        ]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+        ds_in = open_paramfile(PARAMFILE)
+        ds_out = open_paramfile(output_path)
+
+        # Check that included variables/coords match as expected
+        for var in ds_in.variables:
+            if var == "xl":
+                self.assertTrue(np.array_equal(np.array([0.724, 0.87]), ds_out[var].values))
+            elif sp.PFTNAME_VAR in ds_in[var].coords:
+                self.assertTrue(ds_in[var].isel(pft=[0, 1]).equals(ds_out[var]))
+            else:
+                self.assertTrue(ds_in[var].equals(ds_out[var]))
+
+    def test_set_paramfile_changeparam_dbl_onlysomepfts(self):
+        """
+        Test that set_paramfile can (1) copy to a new file with only some requested PFTs and (2)
+        change the values of double parameters of those PFTs
+        """
+        output_path = os.path.join(self.tempdir, "output.nc")
+        pfts_to_include = ["not_vegetated", "needleleaf_evergreen_temperate_tree"]
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "-p",
+            ",".join(pfts_to_include),
+            "xl=0.724,0.87",
+        ]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+        ds_in = open_paramfile(PARAMFILE)
+        ds_out = open_paramfile(output_path)
+
+        # Check that included variables/coords match as expected
+        for var in ds_in.variables:
+            if var == "xl":
+
+                # Changed values (first 2)
+                this_slice = slice(0, 2)
+                expected = np.array([0.724, 0.87])
+                result = ds_out[var].isel(pft=this_slice).values
+                self.assertTrue(np.array_equal(expected, result))
+
+                # Preserved values (everything but the first 2)
+                this_slice = slice(2, None)
+                expected = ds_in["xl"].isel(pft=this_slice)
+                result = ds_out["xl"].isel(pft=this_slice)
+                self.assertTrue(are_paramfile_dataarrays_identical(expected, result))
+            else:
+                self.assertTrue(are_paramfile_dataarrays_identical(ds_in[var], ds_out[var]))
+
+    def test_set_paramfile_extractpfts_changeparam_int(self):
+        """
+        Test that set_paramfile can (1) copy to a new file with only some requested PFTs and (2)
+        change the values of integer parameters of those PFTs
+        """
+        output_path = os.path.join(self.tempdir, "output.nc")
+        pfts_to_include = ["not_vegetated", "needleleaf_evergreen_temperate_tree"]
+        this_var = "max_NH_planting_date"
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "-p",
+            ",".join(pfts_to_include),
+            "--drop-other-pfts",
+            f"{this_var}=1986,1987",
+        ]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+        ds_in = open_paramfile(PARAMFILE)
+        ds_out = open_paramfile(output_path)
+
+        # Check that the variable in question is actually an integer to begin with
+        self.assertTrue(sp.is_integer(ds_in[this_var].values))
+        # Also check that it actually differs from our new values
+        self.assertTrue(ds_in[this_var].values[0] != 1986)
+        self.assertTrue(ds_in[this_var].values[1] != 1987)
+
+        # Check that included variables/coords match as expected
+        for var in ds_in.variables:
+            if var == this_var:
+                self.assertTrue(np.array_equal(np.array([1986, 1987]), ds_out[var].values))
+            elif sp.PFTNAME_VAR in ds_in[var].coords:
+                self.assertTrue(
+                    are_paramfile_dataarrays_identical(ds_in[var].isel(pft=[0, 1]), ds_out[var])
+                )
+            else:
+                self.assertTrue(are_paramfile_dataarrays_identical(ds_in[var], ds_out[var]))
+
+    def test_set_paramfile_fill_value_scalar_double_nan(self):
+        """
+        Test that setting scalar double to fill value writes a literal NaN if that's the _FillValue
+        """
+        # Create paramfile with a double variable with fill value NaN
+        input_path = os.path.join(self.tempdir, "input.nc")
+        ds = open_paramfile(PARAMFILE, mask_and_scale=True)
+        new_param_name = "new_param_abc123"
+        ds[new_param_name] = xr.DataArray(data=np.array(3.14))
+        ds[new_param_name].encoding["_FillValue"] = np.nan
+        self.assertTrue(new_param_name in ds)
+        ds.to_netcdf(input_path)
+
+        # Check that its fill value is NaN
+        ds_in = open_paramfile(input_path, mask_and_scale=True)
+        self.assertTrue("_FillValue" in ds_in[new_param_name].encoding)
+        self.assertTrue(np.isnan(ds_in[new_param_name].encoding["_FillValue"]))
+
+        # Ask to set it to the FillValue
+        output_path = os.path.join(self.tempdir, "output.nc")
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            input_path,
+            "-o",
+            output_path,
+            f"{new_param_name}=nan",
+        ]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+
+        # Ensure it wrote a literal NaN
+        ds_out = open_paramfile(output_path, mask_and_scale=False)
+        self.assertTrue(np.isnan(ds_out[new_param_name]))
+
+        # Ensure it preserved NaN FillValue
+        ds_out = open_paramfile(output_path, mask_and_scale=True)
+        self.assertTrue(np.isnan(ds_out[new_param_name].encoding["_FillValue"]))
+
+    def test_set_paramfile_fill_value_scalar_double_real(self):
+        """
+        Test that setting scalar double to fill value does NOT write NaN if that's not the
+        _FillValue
+        """
+        # Create paramfile with a double variable with fill value -999.9
+        input_path = os.path.join(self.tempdir, "input.nc")
+        ds = open_paramfile(PARAMFILE, mask_and_scale=True)
+        new_param_name = "new_param_abc123"
+        ds[new_param_name] = xr.DataArray(data=np.array(3.14))
+        fill_value = -999.9
+        ds[new_param_name].encoding["_FillValue"] = fill_value
+        self.assertTrue(new_param_name in ds)
+        ds.to_netcdf(input_path)
+
+        # Check that its fill value is what we asked for
+        ds_in = open_paramfile(input_path, mask_and_scale=True)
+        self.assertTrue("_FillValue" in ds_in[new_param_name].encoding)
+        self.assertEqual(fill_value, ds_in[new_param_name].encoding["_FillValue"])
+
+        # Ask to set it to the FillValue
+        output_path = os.path.join(self.tempdir, "output.nc")
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            input_path,
+            "-o",
+            output_path,
+            f"{new_param_name}=nan",
+        ]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+
+        # Ensure it preserved FillValue
+        ds_out = open_paramfile(output_path, mask_and_scale=True)
+        self.assertEqual(fill_value, ds_out[new_param_name].encoding["_FillValue"])
+        self.assertTrue(np.isnan(ds_out[new_param_name].values))
+
+        # Ensure it wrote the FillValue and not a literal NaN
+        ds_out = open_paramfile(output_path, mask_and_scale=False)
+        self.assertEqual(fill_value, ds_out[new_param_name].values)
+
+    def test_set_paramfile_setparams_scalar_double_tonan_with_nancaps(self):
+        """Test setting scalar double to NaN using 'NaN'"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        this_var = "a_coef"
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            f"{this_var}=NaN",
+        ]
+
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+        ds_out = open_paramfile(output_path, mask_and_scale=True)
+        self.assertTrue(np.isnan(ds_out[this_var]))
+
+    def test_set_paramfile_setparams_pft_double_tonan_with_nan(self):
+        """Test setting PFT-dimensioned double to NaN using 'nan'"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        pfts_to_include = ["not_vegetated", "needleleaf_evergreen_temperate_tree"]
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "-p",
+            ",".join(pfts_to_include),
+            "--drop-other-pfts",
+            "xl=nan,nan",
+            "planting_temp=nan,nan",
+        ]
+
+        # Check that planting_temp is already nan and xl isn't
+        ds_in = open_paramfile(PARAMFILE, mask_and_scale=True)
+        self.assertTrue(all(np.isnan(ds_in["planting_temp"].isel(pft=[0, 1]))))
+        self.assertTrue(not any(np.isnan(ds_in["xl"].isel(pft=[0, 1]))))
+
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+        ds_out = open_paramfile(output_path, mask_and_scale=True)
+        self.assertTrue(all(np.isnan(ds_out["xl"])))
+        self.assertTrue(all(np.isnan(ds_out["planting_temp"])))
+
+    def test_set_paramfile_setparams_nan_but_no_fillvalue(self):
+        """Test that NotImplementedError is given if trying to set NaN but param has no FillValue"""
+
+        # Create paramfile with a double variable without fill value
+        input_path = os.path.join(self.tempdir, "input.nc")
+        ds = open_paramfile(PARAMFILE, mask_and_scale=True)
+        new_param_name = "new_param_abc123"
+        ds[new_param_name] = xr.DataArray(data=np.array(3.14))
+        ds[new_param_name].encoding["_FillValue"] = None
+        self.assertTrue(new_param_name in ds)
+        ds.to_netcdf(input_path)
+
+        # Check that it doesn't have fill value
+        ds_in = open_paramfile(input_path, mask_and_scale=True)
+        self.assertFalse("_FillValue" in ds_in[new_param_name].encoding)
+
+        # Call set_paramfile, trying to set it to NaN
+        output_path = os.path.join(self.tempdir, "output.nc")
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            input_path,
+            "-o",
+            output_path,
+            f"{new_param_name}=nan",
+        ]
+        with self.assertRaisesRegex(
+            NotImplementedError, "Can't set parameter to fill value if it doesn't already have one:"
+        ):
+            sp.main()
+
+    def test_set_paramfile_setparams_scalar_int_tonan_with_nan(self):
+        """Test that NotImplementedError is given if trying to set NaN for an integer"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        this_var = "upplim_destruct_metamorph"
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            f"{this_var}=nan",
+        ]
+
+        with self.assertRaisesRegex(
+            NotImplementedError, "Can't set integer parameter to fill value:"
+        ):
+            sp.main()
+
+    # TODO: Test changing param when extracting just one PFT
+
+    # TODO: Test changing PFT name
+
+    def test_set_paramfile_changeparam_multidim_errors(self):
+        """
+        Test that set_paramfile errors if requesting change of a multi-dimensional parameter. This
+        test will obviously need to be replaced once that functionality is added.
+        """
+        output_path = os.path.join(self.tempdir, "output.nc")
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "mimics_till_decompk_multipliers=dummy",
+        ]
+
+        with self.assertRaises(NotImplementedError):
+            sp.main()
+
+    def test_set_paramfile_setparams_just_one_pft(self):
+        """Test changing just one PFT's value of something without dropping others"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        pft_to_include = "needleleaf_deciduous_boreal_tree"
+        this_var = "rswf_max"
+        new_value = 0.7
+
+        # Ensure it wasn't new_value before
+        ds_in = open_paramfile(PARAMFILE)
+        pft_names = check_pfts_in_paramfile([pft_to_include], ds_in)
+        pft_index = get_selected_pft_indices([pft_to_include], pft_names)[0]
+        self.assertFalse(ds_in[this_var].values[pft_index] == new_value)
+
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "-p",
+            pft_to_include,
+            f"{this_var}={new_value}",
+        ]
+        sp.main()
+
+        ds_out = open_paramfile(output_path)
+        for i, value in enumerate(ds_out[this_var]):
+            if i == pft_index:
+                self.assertTrue(value == new_value)
+            else:
+                self.assertTrue(value == ds_in[this_var].values[i])
+
+    def test_set_paramfile_setparams_just_one_pft_dropothers_noset(self):
+        """Test dropping all but one PFT without changing any parameters"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        pft_to_include = "needleleaf_deciduous_boreal_tree"
+
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "-p",
+            pft_to_include,
+            "--drop-other-pfts",
+        ]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+
+        # Check that the file is just what you get if you drop all but the one PFT
+        ds_in = open_paramfile(PARAMFILE)
+        ds_in_1pft = sp.drop_other_pfts([pft_to_include], ds_in)
+        ds_out = open_paramfile(output_path)
+        self.assertTrue(set(ds_in_1pft.variables) == set(ds_out.variables))
+        for var in ds_in_1pft:
+            self.assertTrue(are_paramfile_dataarrays_identical(ds_in_1pft[var], ds_out[var]))
+        self.assertTrue(ds_in_1pft.equals(ds_out))
+        self.assertEqual(ds_in_1pft.sizes["pft"], 1)
+        self.assertEqual(ds_out.sizes["pft"], 1)
+
+    def test_set_paramfile_setparams_just_one_pft_dropothers_doset(self):
+        """Test dropping all but one PFT, changing one parameter"""
+        output_path = os.path.join(self.tempdir, "output.nc")
+        pft_to_include = "needleleaf_deciduous_boreal_tree"
+        this_var = "rswf_max"
+        new_value = 0.7
+
+        # Ensure it wasn't new_value before
+        ds_in = open_paramfile(PARAMFILE)
+        pft_names = check_pfts_in_paramfile([pft_to_include], ds_in)
+        pft_index = get_selected_pft_indices([pft_to_include], pft_names)[0]
+        self.assertFalse(ds_in[this_var].values[pft_index] == new_value)
+
+        sys.argv = [
+            "set_paramfile",
+            "-i",
+            PARAMFILE,
+            "-o",
+            output_path,
+            "-p",
+            pft_to_include,
+            "--drop-other-pfts",
+            f"{this_var}={new_value}",
+        ]
+        sp.main()
+        self.assertTrue(os.path.exists(output_path))
+
+        # Check that all variables match except for the one we changed
+        ds_out = open_paramfile(output_path)
+        ds_in_1pft = sp.drop_other_pfts([pft_to_include], ds_in)
+        for var in ds_in.variables:
+            da_in = ds_in_1pft[var]
+            da_out = ds_out[var]
+            if var == this_var:
+                self.assertFalse(are_paramfile_dataarrays_identical(da_in, da_out))
+            else:
+                self.assertTrue(are_paramfile_dataarrays_identical(da_in, da_out))
+
+    def test_set_paramfile_int_errors_given_float_point0(self):
+        """
+        Test that set_paramfile errors if given float value ending in .0 for an int field
+        """
+        input_path = os.path.join(self.tempdir, "input.nc")
+        output_path = os.path.join(self.tempdir, "output.nc")
+        param_name = "param0"
+
+        # Save a test paramfile with an int field
+        da = xr.DataArray(data=np.int32(3))
+        ds = xr.Dataset(data_vars={param_name: da})
+        ds.to_netcdf(input_path, encoding={param_name: {"_FillValue": None}})
+        ds_in = open_paramfile(input_path)
+        self.assertTrue(sp.is_integer(ds_in[param_name].values))
+
+        # Try giving it a value ending in .0
+        sys.argv = ["set_paramfile", "-i", input_path, "-o", output_path, f"{param_name}=4.0"]
+        with self.assertRaisesRegex(ValueError, "Invalid assignment to an integer parameter:"):
+            sp.main()
+
+    def test_set_paramfile_int_errors_given_float_point1(self):
+        """
+        Test that set_paramfile errors if given float value NOT ending in .0 for an int field
+        """
+        input_path = os.path.join(self.tempdir, "input.nc")
+        output_path = os.path.join(self.tempdir, "output.nc")
+        param_name = "param0"
+
+        # Save a test paramfile with an int field
+        da = xr.DataArray(data=np.int32(3))
+        ds = xr.Dataset(data_vars={param_name: da})
+        ds.to_netcdf(input_path, encoding={param_name: {"_FillValue": None}})
+        ds_in = open_paramfile(input_path)
+        self.assertTrue(sp.is_integer(ds_in[param_name].values))
+
+        # Try giving it a value ending in .1
+        sys.argv = ["set_paramfile", "-i", input_path, "-o", output_path, f"{param_name}=4.1"]
+        with self.assertRaisesRegex(ValueError, "Invalid assignment to an integer parameter:"):
+            sp.main()
+
+    def test_set_paramfile_double_ok_given_int(self):
+        """
+        Test that set_paramfile works if given int value for a double field
+        """
+        input_path = os.path.join(self.tempdir, "input.nc")
+        output_path = os.path.join(self.tempdir, "output.nc")
+        param_name = "param0"
+
+        # Save a test paramfile with a double field
+        da = xr.DataArray(data=np.float32(3.14))
+        ds = xr.Dataset(data_vars={param_name: da})
+        ds.to_netcdf(input_path, encoding={param_name: {"_FillValue": None}})
+        ds_in = open_paramfile(input_path)
+        self.assertFalse(sp.is_integer(ds_in[param_name].values))
+
+        # Give it an integer
+        sys.argv = ["set_paramfile", "-i", input_path, "-o", output_path, f"{param_name}=4"]
+        sp.main()
+
+        # Check that it's still a double after saving
+        ds_out = open_paramfile(output_path)
+        self.assertFalse(sp.is_integer(ds_out[param_name].values))
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_unit_args_utils.py b/python/ctsm/test/test_unit_args_utils.py
index 548e7d9389..b4b101dda6 100755
--- a/python/ctsm/test/test_unit_args_utils.py
+++ b/python/ctsm/test/test_unit_args_utils.py
@@ -17,6 +17,7 @@
 
 # pylint: disable=wrong-import-position
 from ctsm.args_utils import plon_type, plat_type
+from ctsm.args_utils import comma_separated_list
 from ctsm import unit_testing
 
 # pylint: disable=invalid-name
@@ -118,6 +119,48 @@ def test_platType_outOfBounds_negative(self):
             _ = plat_type(-91)
 
 
+class TestArgsCommaSeparatedList(unittest.TestCase):
+    """
+    Test comma_separated_list argparse helper
+    """
+
+    def setUp(self):
+        self.orig_argv = sys.argv
+
+    def tearDown(self):
+        sys.argv = self.orig_argv
+
+    def test_comma_separated_list_1(self):
+        """
+        Test comma_separated_list with one item in list
+        """
+        parser = argparse.ArgumentParser()
+        parser.add_argument("--list-arg", type=comma_separated_list)
+        sys.argv = ["scriptname", "--list-arg", "abc"]
+        args = parser.parse_args()
+        self.assertEqual(["abc"], args.list_arg)
+
+    def test_comma_separated_list_1plus(self):
+        """
+        Test comma_separated_list with one item in list but a comma too
+        """
+        parser = argparse.ArgumentParser()
+        parser.add_argument("--list-arg", type=comma_separated_list)
+        sys.argv = ["scriptname", "--list-arg", "abc,"]
+        args = parser.parse_args()
+        self.assertEqual(["abc", ""], args.list_arg)
+
+    def test_comma_separated_list_2(self):
+        """
+        Test comma_separated_list with two items in list
+        """
+        parser = argparse.ArgumentParser()
+        parser.add_argument("--list-arg", type=comma_separated_list)
+        sys.argv = ["scriptname", "--list-arg", "abc,def"]
+        args = parser.parse_args()
+        self.assertEqual(["abc", "def"], args.list_arg)
+
+
 if __name__ == "__main__":
     unit_testing.setup_for_tests()
     unittest.main()
diff --git a/python/ctsm/test/test_unit_netcdf_utils.py b/python/ctsm/test/test_unit_netcdf_utils.py
new file mode 100755
index 0000000000..50fc2d2c81
--- /dev/null
+++ b/python/ctsm/test/test_unit_netcdf_utils.py
@@ -0,0 +1,390 @@
+#!/usr/bin/env python3
+"""
+Unit tests for netcdf_utils.py functions
+"""
+
+import os
+import sys
+import unittest
+import shutil
+import tempfile
+import numpy as np
+import pandas as pd
+import xarray as xr
+
+# -- add python/ctsm  to path (needed if we want to run the test stand-alone)
+_CTSM_PYTHON = os.path.join(os.path.dirname(os.path.realpath(__file__)), os.pardir, os.pardir)
+sys.path.insert(1, _CTSM_PYTHON)
+
+# pylint: disable=wrong-import-position
+import ctsm.netcdf_utils as nu
+from ctsm import unit_testing
+
+# pylint: disable=invalid-name
+# pylint: disable=protected-access
+
+
+class TestUnitGetNetcdfFormat(unittest.TestCase):
+    """
+    Unit tests for get_netcdf_format
+    """
+
+    def setUp(self):
+        self.tempdir = tempfile.mkdtemp()
+        self.outfile = os.path.join(self.tempdir, "file.nc")
+        da = xr.DataArray(data=[1, 2, 3])
+        self.ds = xr.Dataset(data_vars={"var": da})
+
+    def tearDown(self):
+        shutil.rmtree(self.tempdir, ignore_errors=True)
+
+    def test_get_netcdf_format_classic(self):
+        """
+        Test that get_netcdf_format() gets "classic" format right
+        """
+        nc_format = "NETCDF3_CLASSIC"
+        self.ds.to_netcdf(self.outfile, format=nc_format)
+        self.assertEqual(nu.get_netcdf_format(self.outfile), nc_format)
+
+    def test_get_netcdf_format_netcdf4(self):
+        """
+        Test that get_netcdf_format() gets "netCDF4" format right
+        """
+        nc_format = "NETCDF4"
+        self.ds.to_netcdf(self.outfile, format=nc_format)
+        self.assertEqual(nu.get_netcdf_format(self.outfile), nc_format)
+
+
+class TestUnitAreXrDataArraysIdentical(unittest.TestCase):
+    """
+    Unit tests for are_xr_dataarrays_identical
+    """
+
+    # pylint: disable=too-many-public-methods
+
+    def test_are_xr_dataarrays_identical_data_types_match(self):
+        """Should be true if dtypes match"""
+        da0 = xr.DataArray(int(1))
+        da1 = xr.DataArray(int(1))
+        self.assertTrue(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_data_types_differ(self):
+        """Should be false if dtypes differ"""
+        da0 = xr.DataArray(int(1))
+        da1 = xr.DataArray(np.float64(1))
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_data_types_differ_precision(self):
+        """Should be false if dtypes differ only in their precision"""
+        da0 = xr.DataArray(np.float64(1))
+        da1 = xr.DataArray(np.float32(1))
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_encodings_match(self):
+        """Should be true if encodings are the exact same"""
+        da0 = xr.DataArray(int(1))
+        da1 = xr.DataArray(int(1))
+        da0.encoding["dummy0"] = -999
+        da1.encoding["dummy0"] = -999
+        da0.encoding["dummy1"] = -999
+        da1.encoding["dummy1"] = -999
+        self.assertTrue(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_encodings_match_except_order(self):
+        """Should be true if encodings are the same as long as you don't care about order"""
+        da0 = xr.DataArray(int(1))
+        da1 = xr.DataArray(int(1))
+        da0.encoding["dummy0"] = -999
+        da0.encoding["dummy1"] = -999
+        da1.encoding["dummy1"] = -999
+        da1.encoding["dummy0"] = -999
+        self.assertTrue(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_encodings_differ_number(self):
+        """Should be false if encodings have a different number of keys"""
+        da0 = xr.DataArray(int(1))
+        da1 = xr.DataArray(int(1))
+        da0.encoding["dummy0"] = -999
+        da1.encoding["dummy0"] = -999
+        da1.encoding["dummy1"] = -999
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_encodings_differ(self):
+        """Should be false if encodings have the same number and order of keys but not values"""
+        da0 = xr.DataArray(int(1))
+        da1 = xr.DataArray(int(1))
+        da0.encoding["dummy0"] = -999
+        da1.encoding["dummy0"] = -999
+        da0.encoding["dummy1"] = -998
+        da1.encoding["dummy1"] = -999
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_attributes_match(self):
+        """Should be true if attributes are the exact same"""
+        da0 = xr.DataArray(int(1))
+        da1 = xr.DataArray(int(1))
+        da0.attrs["dummy0"] = -999
+        da1.attrs["dummy0"] = -999
+        da0.attrs["dummy1"] = -999
+        da1.attrs["dummy1"] = -999
+        self.assertTrue(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_attributes_match_nan(self):
+        """Should be true if attributes are the exact same and NaN"""
+        da0 = xr.DataArray(1.0)
+        da1 = xr.DataArray(1.0)
+        da0.attrs["_FillValue"] = np.float64(np.nan)
+        da1.attrs["_FillValue"] = np.float64(np.nan)
+        self.assertTrue(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_attributes_match_except_order(self):
+        """Should be true if attributes are the same as long as you don't care about order"""
+        da0 = xr.DataArray(int(1))
+        da1 = xr.DataArray(int(1))
+        da0.attrs["dummy0"] = -999
+        da0.attrs["dummy1"] = -999
+        da1.attrs["dummy1"] = -999
+        da1.attrs["dummy0"] = -999
+        self.assertTrue(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_attributes_differ_number(self):
+        """Should be false if attributes have a different number of keys"""
+        da0 = xr.DataArray(int(1))
+        da1 = xr.DataArray(int(1))
+        da0.attrs["dummy0"] = -999
+        da1.attrs["dummy0"] = -999
+        da1.attrs["dummy1"] = -999
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_attributes_differ(self):
+        """Should be false if attributes have the same number and order of keys but not values"""
+        da0 = xr.DataArray(int(1))
+        da1 = xr.DataArray(int(1))
+        da0.attrs["dummy0"] = -999
+        da1.attrs["dummy0"] = -999
+        da0.attrs["dummy1"] = -998
+        da1.attrs["dummy1"] = -999
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_values_ndarrays_match(self):
+        """Should be true if values match and they're both numpy arrays under the hood"""
+        da0 = xr.DataArray(np.array(int(1)))
+        da1 = xr.DataArray(np.array(int(1)))
+        self.assertTrue(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_values_ndarrays_differ(self):
+        """Should be false if values differ and they're both numpy arrays under the hood"""
+        da0 = xr.DataArray(np.array(int(1)))
+        da1 = xr.DataArray(np.array(int(2)))
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_values_ndarrays_differ_nantypeerror(self):
+        """
+        Should be false if values differ, they're both numpy arrays under the hood, and they can't
+        be coerced to a type capable of NaN
+        """
+        da0 = xr.DataArray(np.array(["a", "b", "c"]))
+        da1 = xr.DataArray(np.array(["d", "e", "f"]))
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_coords_match(self):
+        """Should be true if coordinates match"""
+        time = pd.date_range("2000-01-01", periods=3)
+        da0 = xr.DataArray(
+            dims=["time"],
+            coords={"time": time},
+        )
+        da1 = xr.DataArray(
+            dims=["time"],
+            coords={"time": time},
+        )
+        self.assertTrue(nu.are_xr_dataarrays_identical(da0, da1))
+
+        # To make sure we don't need a separate "test_are_xr_dataarrays_identical_indexes_match"
+        self.assertEqual(len(da0.indexes), len(da1.indexes))
+        for key in da0.indexes:
+            self.assertTrue(da0.indexes[key].equals(da1.indexes.get(key)))
+
+    def test_are_xr_dataarrays_identical_coords_onemissing(self):
+        """Should be false if only one has coords"""
+        data = [1, 2, 3]
+        time = pd.date_range("2000-01-01", periods=len(data))
+        da0 = xr.DataArray(
+            data=data,
+            dims=["time"],
+            coords={"time": time},
+        )
+        da1 = xr.DataArray(
+            data=data,
+            dims=["time"],
+        )
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_coords_differ(self):
+        """Should be false if coordinates differ"""
+        time0 = pd.date_range("2000-01-01", periods=3)
+        da0 = xr.DataArray(
+            dims=["time"],
+            coords={"time": time0},
+        )
+        time1 = pd.date_range("1987-01-01", periods=3)
+        da1 = xr.DataArray(
+            dims=["time"],
+            coords={"time": time1},
+        )
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_dims_match(self):
+        """Should be true if dimensions match"""
+        da0 = xr.DataArray(data=[1], dims=["dim"])
+        da1 = xr.DataArray(data=[1], dims=["dim"])
+        self.assertTrue(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_dims_differ(self):
+        """Should be false if dimensions differ"""
+        da0 = xr.DataArray(data=[1], dims=["dim0"])
+        da1 = xr.DataArray(data=[1], dims=["dim1"])
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_names_match(self):
+        """Should be true if names match"""
+        da0 = xr.DataArray(data=[1], name="da_name")
+        da1 = xr.DataArray(data=[1], name="da_name")
+        self.assertTrue(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_names_differ(self):
+        """Should be false if names differ"""
+        da0 = xr.DataArray(data=[1], name="da0_name")
+        da1 = xr.DataArray(data=[1], name="da1_name")
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    def test_are_xr_dataarrays_identical_sizes_differ(self):
+        """Should be false if sizes differ"""
+        da0 = xr.DataArray(data=[1])
+        da1 = xr.DataArray(data=[1, 1])
+        self.assertFalse(nu.are_xr_dataarrays_identical(da0, da1))
+
+    # Waiting on dask, sparse, or pint to be in ctsm_pylib:
+    # TODO: False if data types don't match
+    # TODO: NotImplementedError if data types match but aren't np.array
+
+    # Waiting on dask to be in ctsm_pylib:
+    # TODO: True if the only difference is chunked or not
+    # TODO: True if the only difference is chunk sizes
+
+
+class TestUnitAreDictsIdenticalNansEqual(unittest.TestCase):
+    """
+    Unit tests for _are_dicts_identical_nansequal
+    """
+
+    def test_are_dicts_identical_nansequal_yes_noignore_nonan(self):
+        """
+        Test two identical dicts with no keys being ignored and no nans
+        """
+        dict0 = {"a": 1, "b": 2}
+        dict1 = {"a": 1, "b": 2}
+        self.assertTrue(nu._are_dicts_identical_nansequal(dict0, dict1))
+
+    def test_are_dicts_identical_nansequal_yes_ignorestr_nonan(self):
+        """
+        Test two identical dicts with a key being ignored (as str) and no nans
+        """
+        dict0 = {"a": 1, "b": 2}
+        dict1 = {"a": 1, "b": 3}
+        self.assertTrue(nu._are_dicts_identical_nansequal(dict0, dict1, keys_to_ignore="b"))
+
+    def test_are_dicts_identical_nansequal_yes_ignorelist_nonan(self):
+        """
+        Test two identical dicts with a key being ignored (as list) and no nans
+        """
+        dict0 = {"a": 1, "b": 2}
+        dict1 = {"a": 1, "b": 3}
+        self.assertTrue(nu._are_dicts_identical_nansequal(dict0, dict1, keys_to_ignore=["b"]))
+
+    def test_are_dicts_identical_nansequal_yes_ignorestr_nan(self):
+        """
+        Test two identical dicts with a key being ignored (as str) and one key with matching nan
+        values
+        """
+        dict0 = {"a": np.nan, "b": 2}
+        dict1 = {"a": np.nan, "b": 3}
+        self.assertTrue(nu._are_dicts_identical_nansequal(dict0, dict1, keys_to_ignore="b"))
+
+    def test_are_dicts_identical_nansequal_no_nonan(self):
+        """
+        Test two different dicts with no nans
+        """
+        dict0 = {"a": 1, "b": 2}
+        dict1 = {"a": 1, "b": 4}
+        self.assertFalse(nu._are_dicts_identical_nansequal(dict0, dict1))
+
+    def test_are_dicts_identical_nansequal_no_nansmatch(self):
+        """
+        Test two different dicts with matching nans
+        """
+        dict0 = {"a": np.nan, "b": 2}
+        dict1 = {"a": np.nan, "b": 4}
+        self.assertFalse(nu._are_dicts_identical_nansequal(dict0, dict1))
+
+    def test_are_dicts_identical_nansequal_no_nansdiffer(self):
+        """
+        Test two different dicts that only differ in that one has a nan where another doesn't
+        """
+        dict0 = {"a": 1, "b": 2}
+        dict1 = {"a": 1, "b": np.nan}
+        self.assertFalse(nu._are_dicts_identical_nansequal(dict0, dict1))
+
+    def test_are_dicts_identical_nansequal_no_keysdiffer(self):
+        """
+        Test two different dicts that have different keys
+        """
+        dict0 = {"a": 1, "b": 2}
+        dict1 = {"a": 1, "c": 2}
+        self.assertFalse(nu._are_dicts_identical_nansequal(dict0, dict1))
+
+    def test_are_dicts_identical_nansequal_no_lengthsdiffer(self):
+        """
+        Test two different dicts that have different lengths
+        """
+        dict0 = {"a": 1, "b": 2}
+        dict1 = {"a": 1}
+        self.assertFalse(nu._are_dicts_identical_nansequal(dict0, dict1))
+
+    def test_are_dicts_identical_nansequal_no_both_nparrays(self):
+        """
+        Test two different dicts that have differing numpy arrays for one value
+        """
+        dict0 = {"a": 1, "b": np.array([1, 2])}
+        dict1 = {"a": 1, "b": np.array([1, 3])}
+        self.assertFalse(nu._are_dicts_identical_nansequal(dict0, dict1))
+
+    def test_are_dicts_identical_nansequal_no_differ_nparrays(self):
+        """
+        Test two dicts where one has a value that's a numpy array and the other doesn't, but they're
+        identical if you coerce them both to numpy arrays
+        """
+        dict0 = {"a": 1, "b": np.array([1, 2])}
+        dict1 = {"a": 1, "b": [1, 2]}
+        self.assertTrue(nu._are_dicts_identical_nansequal(dict0, dict1))
+
+    def test_are_dicts_identical_nansequal_yes_both_nparrays(self):
+        """
+        Test two different dicts that have identical numpy arrays for one value
+        """
+        dict0 = {"a": 1, "b": np.array([1, 2])}
+        dict1 = {"a": 1, "b": np.array([1, 2])}
+        self.assertTrue(nu._are_dicts_identical_nansequal(dict0, dict1))
+
+    def test_are_dicts_identical_nansequal_yes_both_nparrays_str(self):
+        """
+        Test two different dicts that have identical numpy arrays of strings for one value
+        """
+        dict0 = {"a": 1, "b": np.array(["1", "2"])}
+        dict1 = {"a": 1, "b": np.array(["1", "2"])}
+        self.assertTrue(nu._are_dicts_identical_nansequal(dict0, dict1))
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_unit_query_paramfile.py b/python/ctsm/test/test_unit_query_paramfile.py
new file mode 100755
index 0000000000..b796253ace
--- /dev/null
+++ b/python/ctsm/test/test_unit_query_paramfile.py
@@ -0,0 +1,113 @@
+#!/usr/bin/env python3
+
+"""Unit tests for query_paramfile"""
+
+import unittest
+import sys
+import io
+from contextlib import redirect_stdout
+import xarray as xr
+
+from ctsm import unit_testing
+
+from ctsm.param_utils import query_paramfile as qp
+from ctsm.param_utils.paramfile_shared import PFTNAME_VAR
+
+# Allow names that pylint doesn't like, because otherwise I find it hard
+# to make readable unit test names
+# pylint: disable=invalid-name
+
+
+def _setup_pft_parameter_ds():
+    """
+    Set up a parameter Dataset with a PFT-dimensioned parameter
+    """
+    pft_dimname = "pft"
+    pft_names_list = ["pft0", "pft1"]
+    pft_names_da = xr.DataArray(
+        name=PFTNAME_VAR,
+        data=pft_names_list,
+        dims=[pft_dimname],
+        coords={pft_dimname: pft_names_list},
+    )
+    var_name = "pft_param"
+    pft_param_da = xr.DataArray(
+        data=[1986.0325, 1987.0724], dims=[pft_dimname], coords={PFTNAME_VAR: pft_names_da}
+    )
+    ds = xr.Dataset(data_vars={var_name: pft_param_da, pft_dimname: pft_names_da})
+    return ds, var_name, pft_names_list
+
+
+class TestUnitQueryParamfile(unittest.TestCase):
+    """Unit tests of query_paramfile"""
+
+    def setUp(self):
+        self.orig_argv = sys.argv
+
+    def tearDown(self):
+        sys.argv = self.orig_argv
+
+    def test_query_paramfile_args_short(self):
+        """Test that all arguments can be set correctly with shortnames"""
+        input_path = "/path/to/input.nc"
+        sys.argv = ["get_arguments", "-i", input_path, "-p", "pft1,pft2", "var1", "var2"]
+        args = qp.get_arguments()
+        self.assertEqual(input_path, args.input)
+        self.assertEqual(["pft1", "pft2"], args.pft)
+        self.assertEqual(["var1", "var2"], args.variables)
+
+    def test_query_paramfile_args_long(self):
+        """Test that all arguments can be set correctly with longnames"""
+        input_path = "/path/to/input.nc"
+        sys.argv = ["get_arguments", "--input", input_path, "--pft", "pft1,pft2", "var1", "var2"]
+        args = qp.get_arguments()
+        self.assertEqual(input_path, args.input)
+        self.assertEqual(["pft1", "pft2"], args.pft)
+        self.assertEqual(["var1", "var2"], args.variables)
+
+    def test_query_paramfile_print_scalar(self):
+        """Test that print_values works with a scalar parameter"""
+        scalar_da = xr.DataArray(data=1987.0724)
+        var_name = "scalar_param"
+        ds = xr.Dataset(data_vars={var_name: scalar_da})
+
+        f = io.StringIO()
+        with redirect_stdout(f):
+            qp.print_values(ds, var_name, selected_pfts=None, pft_names=None)
+        out = f.getvalue()
+        self.assertEqual("scalar_param: 1987.0724\n", out)
+
+    def test_query_paramfile_print_pfts_selectnone(self):
+        """Test that print_values works with PFT-dimensioned parameter, selecting no PFTs"""
+        ds, var_name, pft_names_list = _setup_pft_parameter_ds()
+
+        f = io.StringIO()
+        with redirect_stdout(f):
+            qp.print_values(ds, var_name, selected_pfts=None, pft_names=pft_names_list)
+        out = f.getvalue()
+        self.assertEqual("pft_param:\n   pft0: 1986.0325\n   pft1: 1987.0724\n", out)
+
+    def test_query_paramfile_print_pfts_selectall(self):
+        """Test that print_values works with PFT-dimensioned parameter, selecting all PFTs"""
+        ds, var_name, pft_names_list = _setup_pft_parameter_ds()
+
+        f = io.StringIO()
+        with redirect_stdout(f):
+            qp.print_values(ds, var_name, selected_pfts=pft_names_list, pft_names=pft_names_list)
+        out = f.getvalue()
+        self.assertEqual("pft_param:\n   pft0: 1986.0325\n   pft1: 1987.0724\n", out)
+
+    def test_query_paramfile_print_pfts_selectone(self):
+        """Test that print_values works with PFT-dimensioned parameter, selecting one PFT"""
+        ds, var_name, pft_names_list = _setup_pft_parameter_ds()
+
+        f = io.StringIO()
+        with redirect_stdout(f):
+            qp.print_values(ds, var_name, selected_pfts=["pft0"], pft_names=pft_names_list)
+        out = f.getvalue()
+        self.assertEqual("pft_param:\n   pft0: 1986.0325\n", out)
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_unit_set_paramfile.py b/python/ctsm/test/test_unit_set_paramfile.py
new file mode 100755
index 0000000000..6b091f8967
--- /dev/null
+++ b/python/ctsm/test/test_unit_set_paramfile.py
@@ -0,0 +1,281 @@
+#!/usr/bin/env python3
+
+"""Unit tests for set_paramfile"""
+
+import unittest
+import os
+import sys
+import tempfile
+import shutil
+import numpy as np
+import xarray as xr
+
+from ctsm import unit_testing
+
+from ctsm.netcdf_utils import get_netcdf_format
+from ctsm.param_utils import set_paramfile as sp
+from ctsm.param_utils.paramfile_shared import open_paramfile, are_paramfile_dataarrays_identical
+
+# Allow names that pylint doesn't like, because otherwise I find it hard
+# to make readable unit test names
+# pylint: disable=invalid-name
+
+
+PARAMFILE = os.path.join(
+    os.path.dirname(__file__), "testinputs", "ctsm5.3.041.Nfix_params.v13.c250221_upplim250.nc"
+)
+
+
+def save_paramfile_with_integer_that_has_fillvalue(tempdir):
+    """
+    Convenience function for creating a parameter file that has an integer parameter with a fill
+    value
+    """
+    input_path = os.path.join(tempdir, "input.nc")
+    new_param_name = "new_param_abc123"
+    fill_value = -999
+
+    # Construct the Dataset
+    data = np.array([1, 2, fill_value, 4, 5], dtype=np.int32)
+    da = xr.DataArray(data)
+    da.encoding["_FillValue"] = fill_value
+    ds = xr.Dataset(data_vars={new_param_name: da})
+
+    # Check it
+    assert new_param_name in ds
+    assert sp.is_integer(ds[new_param_name].values)
+
+    # Save it
+    sp.save_paramfile(ds, input_path)
+
+    return input_path, new_param_name, fill_value
+
+
+class TestUnitCheckCorrectNdims(unittest.TestCase):
+    """Unit tests of check_correct_ndims"""
+
+    def test_checkcorrectndims_0d_int(self):
+        """Check True when given a standard int for a 0d parameter"""
+        da = xr.DataArray(data=1)
+        self.assertTrue(sp.check_correct_ndims(da, 1))
+
+    def test_checkcorrectndims_0d_int_np(self):
+        """Check True when given a numpy int for a 0d parameter"""
+        da = xr.DataArray(data=1)
+        self.assertTrue(sp.check_correct_ndims(da, np.int32(1)))
+
+    def test_checkcorrectndims_1d_int(self):
+        """Check False when given a standard int for a 0d parameter"""
+        da = xr.DataArray(data=[1, 2])
+        self.assertFalse(sp.check_correct_ndims(da, 1))
+
+    def test_checkcorrectndims_1d_int_np(self):
+        """Check False when given a numpy int for a 0d parameter"""
+        da = xr.DataArray(data=[1, 2])
+        self.assertFalse(sp.check_correct_ndims(da, np.int32(1)))
+
+    def test_checkcorrectndims_0d_list(self):
+        """Check False when given a list for a 0d parameter"""
+        da = xr.DataArray(data=1)
+        self.assertFalse(sp.check_correct_ndims(da, [1, 2]))
+
+    def test_checkcorrectndims_0d_nparray(self):
+        """Check False when given a numpy array for a 0d parameter"""
+        da = xr.DataArray(data=1)
+        self.assertFalse(sp.check_correct_ndims(da, np.array([1, 2])))
+
+    def test_checkcorrectndims_0d_nparray_error(self):
+        """Check for error when given a numpy array for a 0d parameter and requesting throw_error"""
+        da = xr.DataArray(data=1)
+        with self.assertRaisesRegex(RuntimeError, "Incorrect N dims: Expected 0, got 1"):
+            sp.check_correct_ndims(da, np.array([1, 2]), throw_error=True)
+
+    def test_checkcorrectndims_1d_list(self):
+        """Check True when given a list for a 1d parameter"""
+        da = xr.DataArray(data=[1, 2])
+        self.assertTrue(sp.check_correct_ndims(da, [1, 2]))
+
+    def test_checkcorrectndims_1d_nparray(self):
+        """Check True when given a numpy array for a 1d parameter"""
+        da = xr.DataArray(data=[1, 2])
+        self.assertTrue(sp.check_correct_ndims(da, np.array([1, 2])))
+
+
+class TestUnitIsInteger(unittest.TestCase):
+    """Unit tests of is_integer"""
+
+    def test_isinteger_obj_int(self):
+        """Check True if given an object of type int"""
+        self.assertTrue(sp.is_integer(int(1)))
+
+    def test_isinteger_obj_int_np(self):
+        """Check True if given an object of a numpy integer type"""
+        self.assertTrue(sp.is_integer(np.int32(1)))
+
+    def test_isinteger_obj_int_array0d_np(self):
+        """Check True if given an object of a numpy array with integer dtype"""
+        self.assertTrue(sp.is_integer(np.array(1, dtype=np.int32)))
+
+    def test_isinteger_obj_int_0d_xr(self):
+        """Check True if given a numpy scalar integer object via xarray"""
+        da = xr.DataArray(data=1)
+        self.assertTrue(sp.is_integer(da.values))
+
+    def test_isinteger_obj_int_1d_xr(self):
+        """Check True if given a numpy 1-d integer object via xarray"""
+        da = xr.DataArray(data=[1])
+        self.assertTrue(sp.is_integer(da.values))
+        da = xr.DataArray(data=[1, 2])
+        self.assertTrue(sp.is_integer(da.values))
+
+    def test_isinteger_obj_float(self):
+        """Check False if given an object of type float"""
+        self.assertFalse(sp.is_integer(float(3.14)))
+
+    def test_isinteger_obj_float_np(self):
+        """Check False if given an object of a numpy float type"""
+        self.assertFalse(sp.is_integer(np.float32(3.14)))
+
+    def test_isinteger_type_int(self):
+        """Check False if given a type"""
+        self.assertFalse(sp.is_integer(int))
+
+
+class TestUnitSetParamfile(unittest.TestCase):
+    """Unit tests of set_paramfile"""
+
+    def setUp(self):
+        self.orig_argv = sys.argv
+
+    def tearDown(self):
+        sys.argv = self.orig_argv
+
+    def test_set_paramfile_args_short(self):
+        """Test that all arguments can be set correctly with shortnames"""
+        output_path = "/path/to/output.nc"
+        sys.argv = [
+            "get_arguments",
+            "-i",
+            PARAMFILE,
+            "-p",
+            "pft1,pft2",
+            "-o",
+            output_path,
+            "param1=new_value1",
+            "param2=new_value2",
+        ]
+        args = sp.get_arguments()
+        self.assertEqual(PARAMFILE, args.input)
+        self.assertEqual(["pft1", "pft2"], args.pft)
+        self.assertEqual(output_path, args.output)
+        self.assertEqual(["param1=new_value1", "param2=new_value2"], args.param_changes)
+
+    def test_set_paramfile_args_long(self):
+        """Test that all arguments can be set correctly with longnames"""
+        output_path = "/path/to/output.nc"
+        sys.argv = [
+            "get_arguments",
+            "--input",
+            PARAMFILE,
+            "--pft",
+            "pft1,pft2",
+            "--drop-other-pfts",
+            "--output",
+            output_path,
+            "param1=new_value1",
+            "param2=new_value2",
+        ]
+        args = sp.get_arguments()
+        self.assertEqual(PARAMFILE, args.input)
+        self.assertEqual(["pft1", "pft2"], args.pft)
+        self.assertEqual(output_path, args.output)
+        self.assertEqual(["param1=new_value1", "param2=new_value2"], args.param_changes)
+
+    def test_set_paramfile_error_missing_input(self):
+        """Test that it errors if input file doesn't exist"""
+        output_path = "/path/to/output.nc"
+        sys.argv = [
+            "get_arguments",
+            "--input",
+            "nwuefweirbdfdiurbe",
+            "--pft",
+            "pft1,pft2",
+            "--output",
+            output_path,
+        ]
+        with self.assertRaises(FileNotFoundError):
+            sp.get_arguments()
+
+    def test_set_paramfile_error_existing_output(self):
+        """Test that it errors if output file already exists"""
+        sys.argv = [
+            "get_arguments",
+            "--input",
+            PARAMFILE,
+            "--pft",
+            "pft1,pft2",
+            "--output",
+            PARAMFILE,
+        ]
+        with self.assertRaises(FileExistsError):
+            sp.get_arguments()
+
+    def test_set_paramfile_error_dropotherpfts_without_pft(self):
+        """Test that it errors if given --drop-other-pfts without --pft"""
+        output_path = "/path/to/output.nc"
+        sys.argv = [
+            "get_arguments",
+            "--input",
+            PARAMFILE,
+            "--drop-other-pfts",
+            "--output",
+            output_path,
+        ]
+        with self.assertRaises(RuntimeError):
+            sp.get_arguments()
+
+
+class TestUnitSaveParamfile(unittest.TestCase):
+    """Unit tests of save_paramfile"""
+
+    def setUp(self):
+        self.tempdir = tempfile.mkdtemp()
+        self.output_path = os.path.join(self.tempdir, "output.nc")
+
+    def tearDown(self):
+        shutil.rmtree(self.tempdir, ignore_errors=True)
+
+    def test_save_paramfile(self):
+        """Test that save_paramfile can save our usual test file to a new file without changes"""
+        input_path = PARAMFILE
+        ds_in = open_paramfile(input_path)
+        sp.save_paramfile(ds_in, self.output_path, nc_format=get_netcdf_format(input_path))
+        ds_out = open_paramfile(self.output_path)
+        self.assertTrue(ds_out.equals(ds_in))
+        self.assertTrue(set(ds_in.variables) == set(ds_out.variables))
+        for var in ds_in:
+            self.assertTrue(are_paramfile_dataarrays_identical(ds_in[var], ds_out[var]))
+
+    def test_save_paramfile_integer_with_fillvalue(self):
+        """Test that save_paramfile can successfully save an integer parameter with a fill value"""
+
+        # Create paramfile with a integer variable with fill value -999
+        input_path, new_param_name, fill_value = save_paramfile_with_integer_that_has_fillvalue(
+            self.tempdir
+        )
+
+        # Read it, checking that its fill value is what we asked for. Note: We need to mask, because
+        # otherwise the fill value won't be read.
+        ds = xr.open_dataset(input_path, mask_and_scale=True)
+        self.assertTrue("_FillValue" in ds[new_param_name].encoding)
+        self.assertEqual(fill_value, ds[new_param_name].encoding["_FillValue"])
+
+        # Check that the saved variable is an integer type. Note: We need to NOT mask, because
+        # masking converts fill values to NaN, which forces conversion to float.
+        ds = xr.open_dataset(input_path, mask_and_scale=False)
+        self.assertTrue(sp.is_integer(ds[new_param_name].values))
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/testinputs/ctsm5.3.041.Nfix_params.v13.c250221_upplim250.nc b/python/ctsm/test/testinputs/ctsm5.3.041.Nfix_params.v13.c250221_upplim250.nc
new file mode 100644
index 0000000000..ab47c530f2
--- /dev/null
+++ b/python/ctsm/test/testinputs/ctsm5.3.041.Nfix_params.v13.c250221_upplim250.nc
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:539679e78631b252a2a9e00fe73024cbc9fdf455e496921e94483750c6cf610c
+size 200764
diff --git a/tools/param_utils/query_paramfile b/tools/param_utils/query_paramfile
new file mode 100755
index 0000000000..6bde5f09ec
--- /dev/null
+++ b/tools/param_utils/query_paramfile
@@ -0,0 +1,19 @@
+#!/usr/bin/env python3
+"""
+For description and instructions, please see README.
+"""
+
+import os
+import sys
+
+_CTSM_PYTHON = os.path.join(os.path.dirname(os.path.realpath(__file__)),
+                            os.pardir,
+                            os.pardir,
+                            'python')
+sys.path.insert(1, _CTSM_PYTHON)
+
+from ctsm.param_utils.query_paramfile import main
+
+if __name__ == "__main__":
+    main()
+
diff --git a/tools/param_utils/set_paramfile b/tools/param_utils/set_paramfile
new file mode 100755
index 0000000000..cebdb81d30
--- /dev/null
+++ b/tools/param_utils/set_paramfile
@@ -0,0 +1,19 @@
+#!/usr/bin/env python3
+"""
+For description and instructions, please see README.
+"""
+
+import os
+import sys
+
+_CTSM_PYTHON = os.path.join(os.path.dirname(os.path.realpath(__file__)),
+                            os.pardir,
+                            os.pardir,
+                            'python')
+sys.path.insert(1, _CTSM_PYTHON)
+
+from ctsm.param_utils.set_paramfile import main
+
+if __name__ == "__main__":
+    main()
+