diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs
index 3035923261..f86a330408 100644
--- a/.git-blame-ignore-revs
+++ b/.git-blame-ignore-revs
@@ -2,6 +2,7 @@
 b429b63824e09f82e95d2982f14311cbbd8e4a37
 d229b5c6689efc4c2a6cef077515c4ccd5c18ff6
 4cd83cb3ee6d85eb909403487abf5eeaf4d98911
+d229b5c6689efc4c2a6cef077515c4ccd5c18ff6
 0aa2957c1f8603c63fa30b11295c06cfddff44a5
 2cdb380febb274478e84cd90945aee93f29fa2e6
 e44dc469439e02e9ee582dab274d890ebdfab104
diff --git a/cime_config/usermods_dirs/f09_37x288pt_PanBoreal/shell_commands b/cime_config/usermods_dirs/f09_37x288pt_PanBoreal/shell_commands
new file mode 100644
index 0000000000..1176d22871
--- /dev/null
+++ b/cime_config/usermods_dirs/f09_37x288pt_PanBoreal/shell_commands
@@ -0,0 +1,8 @@
+
+# Change below line if you move the subset data directory
+
+./xmlchange CLM_USRDAT_DIR='$DIN_LOC_ROOT/lnd/clm2/regional_datasets/f09_38x288pt_PanBoreal'
+
+./xmlchange ATM_DOMAIN_MESH='$CLM_USRDAT_DIR/domain.lnd.fv0.9x1.25_gx1v7_f09_38x288pt_PanBoreal_c230524_ESMF_UNSTRUCTURED_MESH.nc'
+./xmlchange LND_DOMAIN_MESH='$CLM_USRDAT_DIR/domain.lnd.fv0.9x1.25_gx1v7_f09_38x288pt_PanBoreal_c230524_ESMF_UNSTRUCTURED_MESH.nc'
+./xmlchange MASK_MESH='$CLM_USRDAT_DIR/domain.lnd.fv0.9x1.25_gx1v7_f09_38x288pt_PanBoreal_c230524_ESMF_UNSTRUCTURED_MESH.nc'
diff --git a/cime_config/usermods_dirs/f09_37x288pt_PanBoreal/user_nl_clm b/cime_config/usermods_dirs/f09_37x288pt_PanBoreal/user_nl_clm
new file mode 100644
index 0000000000..025aa390e7
--- /dev/null
+++ b/cime_config/usermods_dirs/f09_37x288pt_PanBoreal/user_nl_clm
@@ -0,0 +1,22 @@
+!----------------------------------------------------------------------------------
+! Users should add all user specific namelist changes below in the form of 
+! namelist_var = new_namelist_value 
+!
+! EXCEPTIONS: 
+! Set use_cndv           by the compset you use and the CLM_BLDNML_OPTS -dynamic_vegetation setting
+! Set use_vichydro       by the compset you use and the CLM_BLDNML_OPTS -vichydro           setting
+! Set use_cn             by the compset you use and CLM_BLDNML_OPTS -bgc  setting
+! Set use_crop           by the compset you use and CLM_BLDNML_OPTS -crop setting
+! Set spinup_state       by the CLM_BLDNML_OPTS -bgc_spinup      setting
+! Set co2_ppmv           with CCSM_CO2_PPMV                      option
+! Set fatmlndfrc         with LND_DOMAIN_PATH/LND_DOMAIN_FILE    options
+! Set finidat            with RUN_REFCASE/RUN_REFDATE/RUN_REFTOD options for hybrid or branch cases
+!                        (includes $inst_string for multi-ensemble cases)
+!                        or with CLM_FORCE_COLDSTART to do a cold start
+!                        or set it with an explicit filename here.
+! Set maxpatch_glc       with GLC_NEC                            option
+! Set glc_do_dynglacier  with GLC_TWO_WAY_COUPLING               env variable
+!----------------------------------------------------------------------------------
+
+
+fsurdat = '$CLM_USRDAT_DIR/surfdata_f09_38x288pt_PanBoreal_hist_16pfts_Irrig_CMIP6_simyr2000_c230523.nc'
diff --git a/doc/design/python_script_user_interface.rst b/doc/design/python_script_user_interface.rst
index 3ad6a4d2cf..87c1b3e8fc 100644
--- a/doc/design/python_script_user_interface.rst
+++ b/doc/design/python_script_user_interface.rst
@@ -81,3 +81,10 @@ More verbose output should go in ``logger.info`` or ``logger.debug`` statements
       * e.g. You might want to output a ``logging.debug`` statement for every variable in a file you are editing.
 
 Near the top of each python module where logging is used, there should be a line, ``logger = logging.getLogger(__name__)``. Then logging statements should be done using statements like ``logger.info(...)``, *not* ``logging.info(...)``: this allows more contextual information in logging output.
+
+====================================================
+ Considerations on inclusion of python packages
+====================================================
+
+Since, this is somewhat an implementation detail the discussion for this is in ``../../python/README.python_pkgs.rst``. The python
+packages used is somewhat both an important part of the user interface, the tool design, and an implementation detail.
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 6c00f5a686..894e9b6c66 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -45,8 +45,8 @@
 # The suffix(es) of source filenames.
 # You can specify multiple suffix as a list of string:
 #
-# source_suffix = ['.rst', '.md']
-source_suffix = '.rst'
+source_suffix = ['.rst', '.md']
+#source_suffix = '.rst'
 
 # The master toctree document.
 master_doc = 'index'
diff --git a/doc/source/how-to-make-mesh.md b/doc/source/how-to-make-mesh.md
new file mode 100644
index 0000000000..0620598e05
--- /dev/null
+++ b/doc/source/how-to-make-mesh.md
@@ -0,0 +1,98 @@
+# Creating an ESMF mesh file from a netCDF file
+
+This gist includes instructions for creating and visualizing a mesh file from a netcdf file with valid 1D or 2D lats and lons coordinates.
+
+* **ESMF Mesh file** aka **Unstructured Grid File Format** is a netcdf file (format) that includes the information about the grids coordinates and their connectivity to each other. 
+
+Additional information about ESMF mesh files are available [here](https://earthsystemmodeling.org/docs/release/ESMF_8_0_1/ESMF_refdoc/node3.html#SECTION03028200000000000000).
+
+------
+
+In this example, we will use `./mesh_maker.py` which uses `mesh_type.py` to create a mesh file and visualize it.
+
+1- First clone my fork and branch that includes these capabilities:
+``` Shell
+git clone https://github.com/negin513/ctsm.git ctsm_mesh
+cd ctsm_mesh
+
+git checkout subset_mesh_dask
+```
+
+2- Next run mesh_maker.py for a netcdf file:
+
+```
+cd tools/site_and_regional
+```
+Check all the avaialble options:
+
+```
+./mesh_maker.py --help
+```
+
+The output shows all available options for this script:
+```
+|------------------------------------------------------------------|
+|---------------------  Instructions  -----------------------------|
+|------------------------------------------------------------------|
+This script creates ESMF unstructured GRID (mesh file) from a netcdf
+file with valid lats and lons. Provided lats and lons can be 1D or 2D.
+
+For example for running WRF-CTSM cases, the user can create a mesh
+file for their domain :
+    ./mesh_maker.py --input wrfinput_d01 --output my_region
+        --lat XLAT --lon XLONG --verbose
+
+optional arguments:
+  -h, --help        show this help message and exit
+  --input INPUT     Netcdf input file for creating ESMF mesh.
+  --output OUTPUT   Name of the ESMF mesh created.
+  --outdir OUT_DIR  Output directory (only if name of output mesh is not
+                    defined)
+  --lat LAT_NAME    Name of latitude varibale on netcdf input file. If none
+                    given, looks to find variables that include 'lat'.
+  --lon LON_NAME    Name of latitude varibale on netcdf input file. If none
+                    given, looks to find variables that include 'lon'.
+  --mask MASK_NAME  Name of mask varibale on netcdf input file. If none given,
+                    create a fake mask with values of 1.
+  --area AREA_NAME  Name of area variable on netcdf input file. If none given,
+                    ESMF calculates element areas automatically.
+  --overwrite       If meshfile exists, overwrite the meshfile.
+  -v, --verbose     Increase output verbosity
+  
+  ```
+  
+Let's create a mesh file from a netcdf file with 1D lats and lons. On the sample files provided 1D lat and long coordinates are saved on `lsmlat` and `lsmlon` variables. 
+
+```
+./mesh_maker.py --input /glade/scratch/negins/example_files/surfdata_4x5_hist_78pfts_CMIP6_simyr1850_275.0-330.0_-40-15_c220705.nc --output test_mesh_1d.nc --lat lsmlat --lon lsmlon --overwrite
+```
+`--verbose` option also provide additional information for debugging. 
+
+This script will create regional and global mesh plots. For example for the above files, the plos are:
+test_mesh_1d_regional.png
+![image](https://user-images.githubusercontent.com/17344536/200441736-972a8136-5c05-4bc9-9bca-b498d972914a.png)
+
+
+test_mesh_1d_global.png
+
+![image](https://user-images.githubusercontent.com/17344536/200441753-d06e95d1-d85b-4216-9c23-d11ba89a31e4.png)
+
+
+
+------
+  ## Creating Mesh files for a WRF domain:
+For running WRF-CTSM cases, we need to create ESMF mesh files for the WRF domain. We can create mesh file from wrfinput (wrf initial condition files). wrfinput has 2D coordinate information on `XLAT` and `XLONG` variable.
+
+For example, let's create a mesh file from a WRF input file for WRF-CTSM run.  
+ ```
+./mesh_maker.py --input /glade/scratch/negins/example_files/wrfinput_d01  --output test_mesh_wrf.nc --lat XLAT --lon XLONG --overwrite
+```
+
+This produce mesh files for running for our WRF domain. 
+
+Here is how the regional plot looks like for this mesh file:
+
+ ![image](https://user-images.githubusercontent.com/17344536/200442002-1ee5595c-9252-4934-a07c-2f6ad86aff1b.png)
+
+
+ 
\ No newline at end of file
diff --git a/py_env_create b/py_env_create
index 4b3612cfda..9a9b8d940c 100755
--- a/py_env_create
+++ b/py_env_create
@@ -20,6 +20,7 @@ if [ $error != 0 ]; then
    echo "For notes on installing on a user system see: https://docs.conda.io/projects/conda/en/latest/user-guide/install/index.html"
    echo "Error code was $error"
    cat condahelp.txt
+   rm condahelp.txt
    exit -1
 fi
 rm condahelp.txt
@@ -44,7 +45,7 @@ usage() {
   echo "[-v|--verbose]  "
   echo "     Run with verbose mode for the install so you see the progress bar"
   echo "[-f|--file <file>]  "
-  echo "     Conda environment file to use (can be a text format or YAML format)"
+  echo "     Conda environment requirements text file to use (text format) in addition to the others"
   echo "     Assumed to be under the directory: $condadir"
   echo "     Default is: $condafile"
   echo "[--option <option>]  "
@@ -80,6 +81,9 @@ while [ $# -gt 0 ]; do
    shift
 done
 
+#
+# Error checking options and setup
+#
 if [ ! -f $condadir/$condafile ]; then
    echo "$condadir/$condafile does NOT exist"
    echo "Use the --file option with a valid filename"
@@ -89,11 +93,13 @@ fi
 echo "Use conda to install the python environment needed to run the CTSM python tools in the conda environment: $ctsm_python"
 echo "Using the file: $condadir/$condafile"
 
+#
 # Check if the environment already exists, if it does continue, if not create it
+#
 conda list -n $ctsm_python >& /dev/null
 if [ $? != 0 ]; then
    echo "Create $ctsm_python"
-   cmd="conda create --force -n $ctsm_python -q"
+   cmd="conda create --force --name $ctsm_python --quiet"
    echo "$cmd"
    $cmd
    if [ $? != 0 ]; then
@@ -103,12 +109,20 @@ if [ $? != 0 ]; then
 else
    echo "$ctsm_python environment already exists"
 fi
-echo "Install $ctsm_python this can take a long time, be patient...."
-verbosity="-q"
+#
+# Install the environemnt
+#
+echo "Install $ctsm_python this can take a long time (12 to 20 minutes is expected), be patient...."
+echo "    ...."
+echo "    ...."
+echo "    ...."
+verbosity="--quiet"
+loglevel="ERROR"
 if [ "$verbose" == "Yes" ]; then
-  verbosity="-v"
+  verbosity="--verbose"
+  loglevel="INFO"
 fi
-cmd="conda install --yes $verbosity -c conda-forge -n $ctsm_python --file $condadir/$condafile $option"
+cmd="conda install --yes $verbosity --channel conda-forge --name $ctsm_python $condadir/$condafile $option"
 echo "$cmd"
 $cmd
 if [ $? != 0 ]; then
@@ -117,6 +131,9 @@ if [ $? != 0 ]; then
    echo "Change the file and try again"
    exit -2
 fi
+#
+# Report on success
+#
 echo "Successfully installed the $ctsm_python python environment"
 echo
 echo "activate the environment by doing the following..."
diff --git a/python/README.python_pkgs.rst b/python/README.python_pkgs.rst
new file mode 100644
index 0000000000..e5a9af742e
--- /dev/null
+++ b/python/README.python_pkgs.rst
@@ -0,0 +1,33 @@
+.. sectnum::
+
+.. contents::
+
+=====================================================================================
+ Requirements to consider for python packages to be included for CTSM python tools
+=====================================================================================
+
+Requirements for CTSM python tools:
+
+Any added dependencies should be discussed and approved by the CTSM-software team. Criteria for evaluation of third-party packages
+include:
+- How much value does this package provide to us beyond what we could get without it?
+- How difficult it is to install this package (including: how much does this complicate / slow down the creation of the conda
+environment)?
+- How stable is this package?
+- How well maintained is this package?
+- Are there other packages that are more stable or better maintained that would provide nearly the same level of value?
+- Tools that require extra packages should be done in a "contrib" type area out of the main part of tools (this would apply for
+advanced plotting capability for example)
+- We need to be able to reproduce working conda environments minimally on our test machines (currently cheyenne and izumi), but also
+on any machine that we run CTSM on. If there is a machine that we can run CTSM on that we can't build the conda environments or run
+the tools on -- that needs to be fixed.
+- We need to tell the user how long to expect the conda environment to load, and give them options if the conda load is taking too
+long
+- Conda environments need to build robustly even for users who don't have ctsm_pylib loaded in their conda environment
+- Currently dask will NOT be something we require for any of the main CTSM tools
+- Currently we won't use conda-lock
+- We specify the black version exactly so that black will function identically for all users
+- We specify the pylint version exactly because pylint is finicky with version and we need it to work identically for all developers
+- We might remove the need for python packages that aren't providing enough utility
+Packages where specific versions seem to be required will have the version requirements in a least the >= form if not an exact
+version
diff --git a/python/conda_env_ctsm_py.txt b/python/conda_env_ctsm_py.txt
index abccf74773..6cdd06e633 100644
--- a/python/conda_env_ctsm_py.txt
+++ b/python/conda_env_ctsm_py.txt
@@ -22,3 +22,5 @@ xesmf
 numba=0.55.2  # Avoid 0.56 until numpy>=1.20.  This is the minimum for xesmf
 pylint=2.8.3
 black=22.3.0  # NOTE: The version here needs to be coordinated with the black github action under ../.github/workflows
+cartopy=0.20.2
+matplotlib=3.3.2
diff --git a/python/conda_env_ctsm_py_latest.txt b/python/conda_env_ctsm_py_latest.txt
index cc442503c4..456dccf34d 100644
--- a/python/conda_env_ctsm_py_latest.txt
+++ b/python/conda_env_ctsm_py_latest.txt
@@ -12,3 +12,5 @@ xarray>=2022.3.0
 pylint>=2.8.3,<2.9.0  # Once, you move off of 2.8.3, make lint shows 2 errors, at 2.11.1 there are 2 more errors beyond that
                       # By, 2.14.4 there are 75 errors
 black>=22.6.0       # NOTE: The version here needs to be coordinated with the black github action under ../.github/workflows
+cartopy=0.20.2
+matplotlib=3.3.2
\ No newline at end of file
diff --git a/python/ctsm/.pylintrc b/python/ctsm/.pylintrc
index ceff04c7d8..db78949288 100644
--- a/python/ctsm/.pylintrc
+++ b/python/ctsm/.pylintrc
@@ -433,6 +433,12 @@ good-names=i,
            k,
            ex,
            Run,
+           ds,
+           m,
+           n,
+           l1,
+           l2,
+           ax,
            _,
 # --- default list is above here, our own list is below here ---
 # Allow logger as a global name in each module, because this seems to follow general recommended convention:
diff --git a/python/ctsm/mesh_maker.py b/python/ctsm/mesh_maker.py
new file mode 100644
index 0000000000..84526f3217
--- /dev/null
+++ b/python/ctsm/mesh_maker.py
@@ -0,0 +1,313 @@
+#!/usr/bin/env python3
+"""
+|------------------------------------------------------------------|
+|---------------------  Instructions  -----------------------------|
+|------------------------------------------------------------------|
+This script creates ESMF unstructured GRID (mesh file) from a netcdf
+file with valid lats and lons. Provided lats and lons can be 1D or 2D.
+
+For example for running WRF-CTSM cases, the user can create a mesh
+file for their domain :
+    ./mesh_maker --input wrfinput_d01 --output my_region \
+        --lat XLAT --lon XLONG --verbose
+
+"""
+import os
+import logging
+import argparse
+import textwrap
+from datetime import datetime
+import xarray as xr
+import numpy as np
+
+from ctsm.site_and_regional.mesh_type import MeshType
+from ctsm.utils import abort
+from ctsm.ctsm_logging import (
+    setup_logging_pre_config,
+    add_logging_args,
+    process_logging_args,
+)
+
+
+def get_parser():
+    """
+    Get the parser object for mesh_maker script.
+
+    Returns:
+        parser (ArgumentParser):
+            ArgumentParser which includes all the parser information.
+    """
+    parser = argparse.ArgumentParser(
+        description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
+    )
+
+    parser.print_usage = parser.print_help
+
+    parser.add_argument(
+        "--input",
+        help="Netcdf input file for creating ESMF mesh.",
+        action="store",
+        dest="input",
+        required=True,
+    )
+    parser.add_argument(
+        "--output",
+        help="Name of the ESMF mesh created.",
+        action="store",
+        dest="output",
+        required=False,
+    )
+
+    parser.add_argument(
+        "--outdir",
+        help="Output directory (only if name of output mesh is not defined)",
+        action="store",
+        dest="out_dir",
+        type=str,
+    )
+
+    parser.add_argument(
+        "--lat",
+        help="Name of latitude variable on netcdf input file.",
+        action="store",
+        dest="lat_name",
+        type=str,
+        required=True,
+    )
+    parser.add_argument(
+        "--lon",
+        help="Name of latitude variable on netcdf input file.",
+        action="store",
+        dest="lon_name",
+        type=str,
+        required=True,
+    )
+    parser.add_argument(
+        "--mask",
+        help="Name of mask variable on netcdf input file."
+        + " If none given, create a fake mask with values of 1.",
+        action="store",
+        dest="mask_name",
+        type=str,
+        required=False,
+    )
+    parser.add_argument(
+        "--area",
+        help="Name of area variable on netcdf input file."
+        + " If none given, ESMF calculates element areas automatically. ",
+        action="store",
+        dest="area_name",
+        type=str,
+        required=False,
+    )
+    parser.add_argument(
+        "--overwrite",
+        help="If meshfile exists, overwrite the meshfile.",
+        action="store_true",
+        dest="overwrite",
+        required=False,
+    )
+
+    add_logging_args(parser)
+    return parser
+
+
+def process_and_check_args(args):
+    """Process and check the arguments"""
+    if args.output and args.out_dir:
+        logging.error(" Both --outdir and --output cannot be provided at the same time.")
+        err_msg = textwrap.dedent(
+            """
+                \n ------------------------------------
+                \n You have provided both --outdir and --output.
+                \n Please provide only one of these options to proceed:
+                \n --outdir : directory to save mesh file. mesh file name automatically created.
+                \n --output : Absolute or relative path of the ESMF mesh file created.\n
+                """
+        )
+        abort(err_msg)
+
+    # -- no file name and output path:
+    if not args.output and not args.out_dir:
+        args.out_dir = os.path.join(os.getcwd(), "meshes")
+
+    if not args.output:
+        # -- make output path if does not exist.
+        if not os.path.isdir(args.out_dir):
+            os.mkdir(args.out_dir)
+
+        today = datetime.today()
+        today_string = today.strftime("%y%m%d")
+        args.output = os.path.join(
+            args.out_dir,
+            os.path.splitext(args.input)[0]
+            + "_ESMF_UNSTRUCTURED_MESH"
+            + "_c"
+            + today_string
+            + ".nc",
+        )
+
+    # -- exit if mesh_out exists and --overwrite is not specified.
+    if os.path.exists(args.output):
+        if args.overwrite:
+            os.remove(args.output)
+        else:
+            err_msg = (
+                "output meshfile exists, please choose --overwrite to overwrite the mesh file."
+            )
+            abort(err_msg)
+
+    if not os.path.isfile(args.input):
+        err_msg = textwrap.dedent(
+            """\
+                \n ------------------------------------
+                \n Input file not found. Please make sure to provide the
+                \n full path of Netcdf input file for making the mesh.
+                \n ------------------------------------
+                """
+        )
+        abort(err_msg)
+
+    return args
+
+
+def check_input_file(args, ds):
+    """Check that the input file has the variables expected"""
+    if args.lat_name not in ds.coords and args.lat_name not in ds.variables:
+        err_msg = "Input file does not have variable named " + args.lat_name
+        abort(err_msg)
+
+    else:
+        logging.debug(
+            "- %s exist in the provided netcdf file with dimension of %s.",
+            args.lat_name,
+            len(ds[args.lat_name].dims).__str__(),
+        )
+
+    if args.lon_name not in ds.coords and args.lon_name not in ds.variables:
+        err_msg = "Input file does not have variable named " + args.lon_name
+        abort(err_msg)
+    else:
+        logging.debug(
+            "- %s exist in the provided netcdf file with dimension of %s.",
+            args.lon_name,
+            len(ds[args.lon_name].dims).__str__(),
+        )
+    if args.mask_name is not None:
+        if args.mask_name not in ds.variables:
+            err_msg = "Input file does not have mask variable named " + args.mask_name
+            abort(err_msg)
+
+    if args.area_name is not None:
+        if args.area_name not in ds.variables:
+            err_msg = "Input file does not have area variable named " + args.area_name
+            abort(err_msg)
+        if "units" not in ds[args.area_name].attrs:
+            err_msg = "Units attribute is NOT on the area variable"
+            abort(err_msg)
+        areaunits = ["radians^2", "radian2", "radians2", "radian^2"]
+        if not any(name in ds[args.area_name].attrs["units"] for name in areaunits):
+            err_msg = "Area does NOT have the correct units of radians^2 but has " + str(
+                ds[args.area_name].attrs["units"]
+            )
+            abort(err_msg)
+    if len(ds[args.lon_name]) == 1:
+        err_msg = "No need to create a mesh file for a single point grid."
+        abort(err_msg)
+
+
+def main():
+    """Main function to create a mesh file from another file"""
+
+    setup_logging_pre_config()
+    parser = get_parser()
+    args = parser.parse_args()
+
+    # --------------------------------- #
+    # process logging args (i.e. debug and verbose)
+    process_logging_args(args)
+
+    args = process_and_check_args(args)
+
+    nc_file = args.input
+    mesh_out = args.output
+    mask_name = args.mask_name
+    area_name = args.area_name
+
+    ds = xr.open_dataset(nc_file, mask_and_scale=False, decode_times=False).transpose()
+
+    check_input_file(args, ds)
+
+    lats = ds[args.lat_name].astype(np.float32)
+    lons = ds[args.lon_name].astype(np.float32)
+
+    if (len(lats.dims) > 2) or (len(lons.dims) > 2):
+        time_dims = [dim for dim in lats.dims if "time" in dim.lower()]
+        if time_dims:
+            logging.debug("- time dimension found on lat %s", str(time_dims))
+            logging.debug("- removing time dimensions from lats and lons. ")
+            lats = lats[:, :, 0]
+            lons = lons[:, :, 0]
+        else:
+            err_msg = (
+                "latitude or longitude has more than 2 dimensions and "
+                + "the third dimension cannot be detected as time."
+            )
+            abort(err_msg)
+
+    logging.info("Creating mesh file from : %s", nc_file)
+    logging.info("Writing mesh file to    : %s", mesh_out)
+
+    if mask_name is not None:
+        mask = ds[mask_name].astype(np.float32)
+        if mask.max() > 1.0 or mask.min() < 0.0:
+            abort("Mask variable is not within 0 to 1")
+    else:
+        mask = None
+
+    if area_name is not None:
+        area = np.array(ds[area_name].astype(np.float32))
+    else:
+        area = None
+
+    this_mesh = MeshType(lats, lons, mask=mask, area=area)
+    this_mesh.calculate_corners()
+    this_mesh.calculate_nodes()
+    this_mesh.create_esmf(mesh_out)
+
+
+def read_main():
+    """Main function to read a mesh file and output it again"""
+
+    setup_logging_pre_config()
+    parser = get_parser()
+    args = parser.parse_args()
+
+    # --------------------------------- #
+    # process logging args (i.e. debug and verbose)
+    process_logging_args(args)
+
+    args = process_and_check_args(args)
+
+    nc_file = args.input
+    mesh_out = args.output
+
+    ds = xr.open_dataset(nc_file, mask_and_scale=False, decode_times=False).transpose()
+
+    lon0 = np.array([120.0])
+    lat0 = np.array([45.0])
+    x_dim = "lon"
+    y_dim = "lat"
+    lons = xr.DataArray(lon0, name="lon", dims=x_dim, coords={x_dim: lon0})
+    lats = xr.DataArray(lat0, name="lat", dims=y_dim, coords={y_dim: lat0})
+
+    logging.info("Reading  mesh file from : %s", nc_file)
+    logging.info("Writing mesh file to    : %s", mesh_out)
+
+    this_mesh = MeshType(lats, lons)
+    this_mesh.read_file(ds)
+    this_mesh.create_esmf(mesh_out)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/python/ctsm/mesh_plotter.py b/python/ctsm/mesh_plotter.py
new file mode 100644
index 0000000000..7bcbae6fb4
--- /dev/null
+++ b/python/ctsm/mesh_plotter.py
@@ -0,0 +1,158 @@
+#!/usr/bin/env python3
+"""
+|------------------------------------------------------------------|
+|---------------------  Instructions  -----------------------------|
+|------------------------------------------------------------------|
+This script plots an ESMF unstructured GRID (mesh file).
+
+"""
+import os
+import logging
+import argparse
+import textwrap
+from datetime import datetime
+import xarray as xr
+import numpy as np
+
+from ctsm.site_and_regional.mesh_plot_type import MeshPlotType
+from ctsm.utils import abort
+from ctsm.ctsm_logging import (
+    setup_logging_pre_config,
+    add_logging_args,
+    process_logging_args,
+)
+
+
+def get_parser():
+    """
+    Get the parser object for mesh_plotter script.
+
+    Returns:
+        parser (ArgumentParser):
+            ArgumentParser which includes all the parser information.
+    """
+    parser = argparse.ArgumentParser(
+        description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
+    )
+
+    parser.print_usage = parser.print_help
+
+    parser.add_argument(
+        "--input",
+        help="Netcdf ESMF Mesh input file to plot.",
+        action="store",
+        dest="input",
+        required=True,
+    )
+    parser.add_argument(
+        "--output",
+        help="Prefix for the names of the plot files to create.",
+        action="store",
+        dest="output",
+        required=False,
+    )
+
+    parser.add_argument(
+        "--outdir",
+        help="Output directory for plots",
+        action="store",
+        dest="out_dir",
+        type=str,
+    )
+
+    parser.add_argument(
+        "--overwrite",
+        help="If plots  xists, overwrite them.",
+        action="store_true",
+        dest="overwrite",
+        required=False,
+    )
+
+    add_logging_args(parser)
+    return parser
+
+
+def process_and_check_args(args):
+    """Process and check the arguments"""
+    if args.output and args.out_dir:
+        logging.error(" Both --outdir and --output cannot be provided at the same time.")
+        err_msg = textwrap.dedent(
+            """
+                \n ------------------------------------
+                \n You have provided both --outdir and --output.
+                \n Please provide only one of these options to proceed:
+                \n --outdir : directory to save mesh file plots, file names assumed.
+                \n --output : Absolute or relative path of the prefix for plot filenames.\n
+                """
+        )
+        abort(err_msg)
+
+    # -- no file name and output path:
+    if not args.output and not args.out_dir:
+        args.out_dir = os.path.join(os.getcwd(), "meshes")
+
+    if not args.output:
+        # -- make output path if does not exist.
+        if not os.path.isdir(args.out_dir):
+            os.mkdir(args.out_dir)
+
+        today = datetime.today()
+        today_string = today.strftime("%y%m%d")
+        args.output = os.path.join(
+            args.out_dir,
+            os.path.splitext(args.input)[0] + "_c" + today_string,
+        )
+
+    if not os.path.isfile(args.input):
+        err_msg = textwrap.dedent(
+            """\
+                \n ------------------------------------
+                \n Input file not found. Please make sure to provide the
+                \n full path of Netcdf input mesh file
+                \n ------------------------------------
+                """
+        )
+        abort(err_msg)
+
+    return args
+
+
+def main():
+    """Main function to read a mesh file and plot it"""
+
+    setup_logging_pre_config()
+    parser = get_parser()
+    args = parser.parse_args()
+
+    # --------------------------------- #
+    # process logging args (i.e. debug and verbose)
+    process_logging_args(args)
+
+    args = process_and_check_args(args)
+
+    nc_file = args.input
+    mesh_out = args.output
+
+    ds = xr.open_dataset(nc_file, mask_and_scale=False, decode_times=False).transpose()
+
+    lon0 = np.array([120.0])
+    lat0 = np.array([45.0])
+    x_dim = "lon"
+    y_dim = "lat"
+    lons = xr.DataArray(lon0, name="lon", dims=x_dim, coords={x_dim: lon0})
+    lats = xr.DataArray(lat0, name="lat", dims=y_dim, coords={y_dim: lat0})
+
+    logging.info("Reading  mesh file from : %s", nc_file)
+
+    this_mesh = MeshPlotType(lats, lons)
+    this_mesh.read_file(ds)
+
+    plot_regional = os.path.splitext(mesh_out)[0] + "_regional" + ".png"
+
+    plot_global = os.path.splitext(mesh_out)[0] + "_global" + ".png"
+
+    this_mesh.make_mesh_plot(plot_regional, plot_global)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/python/ctsm/site_and_regional/mesh_plot_type.py b/python/ctsm/site_and_regional/mesh_plot_type.py
new file mode 100644
index 0000000000..872c0e8101
--- /dev/null
+++ b/python/ctsm/site_and_regional/mesh_plot_type.py
@@ -0,0 +1,136 @@
+"""
+This module extends the mesh type for advanced plotting
+"""
+import logging
+
+import numpy as np
+
+# -- libraries for plotting mesh (make_mesh_plot)
+# Turn import error off in case you are using a
+# conda environment that doesn't include these.
+# pylint: disable=import-error
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
+
+from ctsm.site_and_regional.mesh_type import MeshType
+
+logger = logging.getLogger(__name__)
+
+
+class MeshPlotType(MeshType):
+    """
+    Extend mesh type with some advanced plotting capability
+    """
+
+    def make_mesh_plot(self, plot_regional, plot_global):
+        """
+        Create plots for the ESMF mesh file
+
+        Parameters
+        ----------
+        plot_regional : str
+            The path to write the ESMF meshfile regional plot
+        plot_global : str
+            The path to write the ESMF meshfile global plot
+        """
+
+        self.mesh_plot(plot_regional, regional=True)
+        self.mesh_plot(plot_global, regional=False)
+
+    def mesh_plot(self, plot_file, regional):
+        """Make a plot of a mesh file in either a regional or global grid"""
+        # -- regional settings
+        if regional:
+            plt.figure(num=None, figsize=(15, 13), facecolor="w", edgecolor="k")
+            # pylint: disable=abstract-class-instantiated
+            ax = plt.axes(projection=ccrs.PlateCarree())
+            plot_type = "regional"
+            line_width = 1
+            marker = "x"
+            marker_size = 1
+        # global settings
+        else:
+            fig = plt.figure(num=None, figsize=(15, 10), facecolor="w", edgecolor="k")
+            # pylint: disable=abstract-class-instantiated
+            ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())
+            plot_type = "global"
+            line_width = 0.5
+            marker = "o"
+            marker_size = None
+
+        ax.add_feature(cfeature.COASTLINE, edgecolor="black")
+        ax.add_feature(cfeature.BORDERS, edgecolor="black")
+        ax.add_feature(cfeature.OCEAN)
+        ax.add_feature(cfeature.BORDERS)
+        ax.add_feature(cfeature.LAND, edgecolor="black")
+        ax.add_feature(cfeature.LAKES, edgecolor="black")
+        ax.add_feature(cfeature.RIVERS)
+        ax.gridlines(
+            color="black",
+            linestyle="dotted",
+            draw_labels=True,
+        )
+        if not regional:
+            ax.set_global()
+
+        # -- plot corner coordinates
+        # clats, clons = self.node_coords.T.compute()
+        # elem_conn_vals = self.elem_conn.compute()
+        clats, clons = self.node_coords.T
+        elem_conn_vals = self.elem_conn
+        element_counts = elem_conn_vals.shape[0]
+        # Scale marker and line size for a large number of points
+        if element_counts > 60000:
+            marker_size = 0.01
+            line_width = line_width * 0.05
+
+        for index in range(element_counts):
+            conns = [int(x) - 1 for x in elem_conn_vals[index]]
+
+            lat_corners = clats[conns]
+            lon_corners = clons[conns]
+            poly_corners = np.zeros((len(lat_corners), 2))
+            poly_corners[:, 1] = lon_corners
+            poly_corners[:, 0] = lat_corners
+            # pylint: disable=abstract-class-instantiated
+            poly = mpatches.Polygon(
+                poly_corners,
+                closed=True,
+                ec="black",
+                transform=ccrs.PlateCarree(),
+                zorder=10,
+                facecolor="none",
+                linewidth=line_width,
+            )
+            ax.add_patch(poly)
+
+        # -- plot center coordinates
+        # clon, clat = self.center_coords.T.compute()
+        clon, clat = self.center_coords.T
+
+        # pylint: disable=abstract-class-instantiated
+        ax.scatter(
+            clon,
+            clat,
+            color="tomato",
+            marker=marker,
+            s=marker_size,
+            transform=ccrs.PlateCarree(),
+            zorder=11,
+        )
+        if regional:
+            lc_colors = {
+                "Corner Coordinates": "black",  # value=0
+                "Center Coordinates": "tomato",  # value=1
+            }
+            labels, handles = zip(
+                *[(k, mpatches.Rectangle((0, 0), 1, 1, facecolor=v)) for k, v in lc_colors.items()]
+            )
+
+            ax.legend(handles, labels)
+
+        plt.savefig(plot_file, bbox_inches="tight")
+
+        logger.info("Successfully created %s plots for ESMF Mesh file : %s", plot_type, plot_file)
diff --git a/python/ctsm/site_and_regional/mesh_type.py b/python/ctsm/site_and_regional/mesh_type.py
new file mode 100644
index 0000000000..be785e745d
--- /dev/null
+++ b/python/ctsm/site_and_regional/mesh_type.py
@@ -0,0 +1,580 @@
+"""
+This module includes the definition and functions for defining a Grid or Mesh.
+This enables creating ESMF mesh file (unstructured grid file)for valid 1D or 2D lats and lons.
+"""
+import logging
+import argparse
+import datetime
+
+import numpy as np
+import xarray as xr
+
+# import dask.array as da
+# import dask.dataframe as dd
+import pandas
+
+from ctsm.utils import abort
+
+logger = logging.getLogger(__name__)
+
+
+class MeshType:
+    """
+    An object for describing mesh or grid files.
+    """
+
+    # pylint: disable=too-many-instance-attributes
+    def __init__(self, center_lats, center_lons, mask=None, mesh_name=None, area=None):
+        """
+        Construct a mesh object
+
+        Attributes
+        ----------
+        center_lats : xarray DataArray
+            array of latitudes (either 1D or 2D)
+        center_lons : xarray DataArray
+            array longitudes (either 1D or 2D)
+        mesh_name : str or None, optional
+            Name of the mesh
+        mask : np array or None
+            numpy array that include landmask
+        area : numpy.ndarray or None
+            Array containing element areas for the ESMF mesh file
+            If None, ESMF calculates element areas internally.
+
+        Methods
+        -------
+        check_lat_lon_dims:
+            check if dimensions of latitude and longitude is valid.
+
+        create_artificial_mask:
+            create an artificial mask of ones if mask does not exits.
+
+        create_2d_coords:
+            if 1d coords is provided, this method creates 2d coords from 1d coords.
+
+        calculate_corners:
+            calculate corner coordinates for each polygon in the grid
+
+        calculate_node_coords:
+            extract node coordiantes for nodeCoords variable on mesh file
+
+        calculate_elem_conn
+            calculate element connectivity for elementConn variable on mesh file
+
+        create_esmf
+            write mesh file to netcdf file.
+        """
+
+        self.mesh_name = mesh_name
+        self.center_lats = center_lats
+        self.center_lons = center_lons
+
+        # -- dims of lat and lon (1d or 2d)
+        self.lat_dims = len(self.center_lats.dims)
+        self.lon_dims = len(self.center_lons.dims)
+        self.check_lat_lon_dims()
+
+        if mask is None:
+            self.create_artificial_mask()
+        else:
+            # self.mask = da.from_array(np.array(mask.astype(np.int8)))
+            self.mask = np.array(mask.astype(np.int8))
+
+        self.area = area
+
+        self.unit = "degrees"
+
+        #
+        # variables that will be set to proper values later...
+        #
+        # These are only used for a conversion
+        self.center_lat2d = 1
+        self.center_lon2d = 1
+        self.corner_lats = 1
+        self.corner_lons = 1
+        # These are normally calculated
+        self.node_coords = 1
+        self.elem_conn = 1
+        self.center_coords = 1
+        # This one is only used when read from a file
+        self.num_elem_conn = 1
+
+    def read_file(self, xrds_in):
+        """
+        Read an input mesh file
+        """
+        # Check that corners are set
+        if self.are_corners_set():
+            abort(
+                "Corners have been set in the mesh object,"
+                + " read_file can only be used from an empty mesh_type object"
+            )
+
+        if not isinstance(xrds_in, xr.Dataset):
+            abort("Input file is not a X-Array DataSet type")
+
+        expected_var_list = ["nodeCoords", "elementConn", "numElementConn", "centerCoords"]
+        for var in expected_var_list:
+            if not var in xrds_in.variables:
+                abort("Variable expected to be on an ESMF mesh file is NOT found: " + var)
+
+        # Check for optional fields
+        if "elementMask" in xrds_in.variables:
+            if max(xrds_in["elementMask"]) > 1.0 or min(xrds_in["elementMask"]) < 0.0:
+                abort("elementMask variable is not within 0 to 1")
+            mask = xrds_in["elementMask"].astype(np.float32)
+            self.mask = np.array(mask.astype(np.int8))
+        else:
+            self.create_artificial_mask()
+
+        if "elementArea" in xrds_in.variables:
+            if "units" in xrds_in["elementArea"].attrs:
+                self.area = np.array(xrds_in["elementArea"])
+                if xrds_in["elementArea"].attrs["units"] != "radians^2":
+                    abort("elementArea is on the mesh file, but without the correct units")
+
+        #
+        # Get data
+        #
+        self.center_coords = xrds_in["centerCoords"].T
+        num_element_count = self.center_coords.shape[0]
+        num_coord_dim = self.center_coords.shape[1]
+        self.elem_conn = xrds_in["elementConn"].T
+        if self.elem_conn.shape[0] != num_element_count:
+            abort(
+                "size of elementConn is not consistent is "
+                + str(self.elem_conn.shape[0])
+                + " expected "
+                + str(num_element_count)
+            )
+        self.num_elem_conn = xrds_in["numElementConn"]
+        if self.num_elem_conn.shape[0] != num_element_count:
+            abort(
+                "size of elementConn is not consistent is "
+                + str(self.num_elem_conn.shape[0])
+                + " expected "
+                + str(num_element_count)
+            )
+        self.node_coords = xrds_in["nodeCoords"].T
+        if self.node_coords.shape[1] != num_coord_dim:
+            abort(
+                "size of node_coords is not consistent is "
+                + str(self.node_coords.shape[1])
+                + " expected "
+                + str(num_coord_dim)
+            )
+        if max(self.num_elem_conn) > self.elem_conn.shape[1]:
+            abort(
+                "max size of num_elem_conn is greater than dimension"
+                + str(max(self.num_elem_conn))
+                + " expected "
+                + str(self.elem_conn.shape[1])
+            )
+
+    def check_lat_lon_dims(self):
+        """
+        Check latitude and longitude dimensions to make sure they are valid.
+
+        -------------
+        Raises:
+            Error (ArgumentTypeError):
+                If the provided latitude has dimension >2.
+            Error (ArgumentTypeError):
+                If the provided longitude has dimension >2.
+        """
+
+        if self.lat_dims not in [1, 2]:
+            err_msg = """
+                Unrecognized grid! \n
+                The dimension of latitude should be either 1 or 2 but it is {}.
+                """.format(
+                self.lat_dims
+            )
+            raise argparse.ArgumentTypeError(err_msg)
+
+        if self.lon_dims not in [1, 2]:
+            err_msg = """
+                Unrecognized grid! \n
+                The dimension of longitude should be either 1 or 2 but it is {}.
+                """.format(
+                self.lon_dims
+            )
+            raise argparse.ArgumentTypeError(err_msg)
+
+    def create_artificial_mask(self):
+        """
+        create an artificial mask of 1 (i.e. all ocean) if no land mask is provided.
+        """
+
+        logger.info("Creating an artificial mask for this region...")
+
+        if self.lat_dims == 1:
+            # -- 1D mask (lat x lon)
+            lats_size = self.center_lats.size
+            lons_size = self.center_lons.size
+            mask = np.ones([lons_size, lats_size], dtype=np.int8)
+        elif self.lat_dims == 2:
+            # -- 2D mask
+            mask = np.ones(self.center_lats.shape, dtype=np.int8)
+        # mask_da = da.from_array(mask)
+        mask_da = mask
+        self.mask = mask_da
+
+    def create_2d_coords(self):
+        """
+        Create 2d center points for our mesh
+        and convert them to Dask Array.
+        """
+        self.center_lats = self.center_lats.astype(np.float64, copy=False)
+        self.center_lons = self.center_lons.astype(np.float64, copy=False)
+
+        if self.lat_dims == 1:
+            # -- 1D lats and lons
+            lats_size = self.center_lats.size
+            lons_size = self.center_lons.size
+
+            # -- convert center points from 1d to 2d
+            self.center_lat2d = np.broadcast_to(self.center_lats[:], (lons_size, lats_size))
+            self.center_lon2d = np.broadcast_to(self.center_lons[:], (lons_size, lats_size))
+        elif self.lat_dims == 2:
+            # -- 2D lats and lons
+            dims = self.center_lons.shape
+
+            # -- convert 2D lats and lons to number x and y
+            lons_size = dims[0]
+            lats_size = dims[1]
+
+            # -- convert to dask array
+            # self.center_lat2d = da.from_array(np.array(self.center_lats))
+            # self.center_lon2d = da.from_array(np.array(self.center_lons))
+            self.center_lat2d = np.array(self.center_lats)
+            self.center_lon2d = np.array(self.center_lons)
+
+    def calculate_corners(self, unit="degrees"):
+        """
+        calculate corner coordinates by averaging adjacent cells
+
+        Parameters
+        ----------
+        unit : {'degrees', 'radians'}, optional
+            The unit of corner coordinates.
+        """
+
+        self.create_2d_coords()
+        # -- pad center_lats for calculating edge gridpoints
+        # -- otherwise we cannot calculate the corner coords
+        # -- for the edge rows/columns.
+
+        padded_lat2d = np.pad(self.center_lat2d, (1, 1), mode="reflect", reflect_type="odd")
+
+        # -- pad center_lons for calculating edge grids
+        padded_lon2d = np.pad(self.center_lon2d, (1, 1), mode="reflect", reflect_type="odd")
+
+        # -- calculate corner lats for each grid
+        north_east = (
+            padded_lat2d[1:-1, 1:-1]
+            + padded_lat2d[0:-2, 1:-1]
+            + padded_lat2d[1:-1, 2:]
+            + padded_lat2d[0:-2, 2:]
+        ) / 4.0
+        north_west = (
+            padded_lat2d[1:-1, 1:-1]
+            + padded_lat2d[0:-2, 1:-1]
+            + padded_lat2d[1:-1, 0:-2]
+            + padded_lat2d[0:-2, 0:-2]
+        ) / 4.0
+        south_west = (
+            padded_lat2d[1:-1, 1:-1]
+            + padded_lat2d[1:-1, 0:-2]
+            + padded_lat2d[2:, 1:-1]
+            + padded_lat2d[2:, 0:-2]
+        ) / 4.0
+        south_east = (
+            padded_lat2d[1:-1, 1:-1]
+            + padded_lat2d[1:-1, 2:]
+            + padded_lat2d[2:, 1:-1]
+            + padded_lat2d[2:, 2:]
+        ) / 4.0
+
+        # -- order counter-clockwise
+        self.corner_lats = np.stack(
+            [
+                north_west.T.reshape((-1,)).T,
+                south_west.T.reshape((-1,)).T,
+                south_east.T.reshape((-1,)).T,
+                north_east.T.reshape((-1,)).T,
+            ],
+            axis=1,
+        )
+        # Corners don't go beyond the poles
+        indx = np.argwhere(self.corner_lats > 90.0)
+        self.corner_lats[indx] = 90.0
+        indx2 = np.argwhere(self.corner_lats < -90.0)
+        self.corner_lats[indx2] = 90.0
+
+        # -- calculate corner lons for each grid
+        north_east = (
+            padded_lon2d[1:-1, 1:-1]
+            + padded_lon2d[0:-2, 1:-1]
+            + padded_lon2d[1:-1, 2:]
+            + padded_lon2d[0:-2, 2:]
+        ) / 4.0
+        north_west = (
+            padded_lon2d[1:-1, 1:-1]
+            + padded_lon2d[0:-2, 1:-1]
+            + padded_lon2d[1:-1, 0:-2]
+            + padded_lon2d[0:-2, 0:-2]
+        ) / 4.0
+        south_west = (
+            padded_lon2d[1:-1, 1:-1]
+            + padded_lon2d[1:-1, 0:-2]
+            + padded_lon2d[2:, 1:-1]
+            + padded_lon2d[2:, 0:-2]
+        ) / 4.0
+        south_east = (
+            padded_lon2d[1:-1, 1:-1]
+            + padded_lon2d[1:-1, 2:]
+            + padded_lon2d[2:, 1:-1]
+            + padded_lon2d[2:, 2:]
+        ) / 4.0
+
+        # -- order counter-clockwise
+        self.corner_lons = np.stack(
+            [
+                north_west.T.reshape((-1,)).T,
+                south_west.T.reshape((-1,)).T,
+                south_east.T.reshape((-1,)).T,
+                north_east.T.reshape((-1,)).T,
+            ],
+            axis=1,
+        )
+        # Longitudes should stay within 0 to 360
+        if np.any(self.corner_lons > 360.0):
+            abort("Corners have longitudes greater than 360")
+        if np.any(self.corner_lons < 0.0):
+            logger.warning(
+                "Corners have longitudes less than zero -- %s %s",
+                "this will work, but may result in answers",
+                "different to roundoff to corners within 0 to 360",
+            )
+
+        self.unit = unit
+
+    def calculate_node_coords(self):
+        """
+        Calculates coordinates of each node (for 'nodeCoords' in ESMF mesh).
+        In ESMF mesh, 'nodeCoords' is a two-dimensional array
+        with dimension ('nodeCount','coordDim')
+        """
+        # Check that corners are set
+        if not self.are_corners_set():
+            abort("Corners have not been set in the mesh object")
+
+        # -- create an array of corner pairs
+        corner_pairs = np.stack(
+            [self.corner_lons.T.reshape((-1,)).T, self.corner_lats.T.reshape((-1,)).T],
+            axis=1,
+        )
+
+        # -- convert to float32 to find duplicates
+
+        # -- remove coordinates that are shared between the elements
+        # node_coords = dd.from_dask_array(corner_pairs).drop_duplicates().values
+        # node_coords.compute_chunk_sizes()
+        node_coords = pandas.DataFrame(data=corner_pairs).drop_duplicates().values
+
+        # -- check size of unique coordinate pairs
+        dims = self.mask.shape
+        nlon = dims[0]
+        nlat = dims[1]
+        elem_conn_size = (nlon * nlat) + nlon
+        elem_conn_size = elem_conn_size + nlat + 1
+
+        self.node_coords = node_coords
+
+        # -- error check to avoid issues later
+        if self.node_coords.shape[0] != elem_conn_size:
+            logger.error(
+                "The size of unique coordinate pairs is %d but expected size is %d!",
+                self.node_coords.shape[0],
+                elem_conn_size,
+            )
+            logger.error("Including a pole point in your grid is the usual reason for this problem")
+            abort("Expected size for element connections is wrong")
+
+    def are_corners_set(self):
+        """Check if the corners have been set for this object"""
+        if isinstance(self.corner_lons, int):
+            return False
+        if isinstance(self.corner_lats, int):
+            return False
+        return True
+
+    def calculate_elem_conn(self):
+        """
+        Calculate element connectivity (for 'elementConn' in ESMF mesh).
+        In ESMF mesh, 'elementConn' describes how the nodes are connected together.
+        """
+        # corners = dd.concat(
+        #    [
+        #        dd.from_dask_array(corner)
+        #        for corner in [
+        #            self.corner_lons.T.reshape((-1,)).T,
+        #            self.corner_lats.T.reshape((-1,)).T,
+        #        ]
+        #    ],
+        #    axis=1,
+        # )
+        # corners.columns = ["lon", "lat"]
+
+        # elem_conn = corners.compute().groupby(["lon", "lat"], sort=False).ngroup() + 1
+        # elem_conn = da.from_array(elem_conn.to_numpy())
+
+        ## -- reshape to write to ESMF
+        # self.elem_conn = elem_conn.T.reshape((4, -1)).T
+
+        # Check that corners are set
+        if not self.are_corners_set():
+            abort("Corners have not been set in the mesh object")
+
+        corners = pandas.concat(
+            [
+                pandas.DataFrame(corner)
+                for corner in [
+                    self.corner_lons.T.reshape((-1,)).T,
+                    self.corner_lats.T.reshape((-1,)).T,
+                ]
+            ],
+            axis=1,
+        )
+        corners.columns = ["lon", "lat"]
+
+        elem_conn = corners.groupby(["lon", "lat"], sort=False).ngroup() + 1
+        elem_conn = elem_conn.to_numpy()
+
+        # -- reshape to write to ESMF
+        self.elem_conn = elem_conn.T.reshape((4, -1)).T
+
+    def are_nodes_set(self):
+        """Check if the nodes have been set for this object"""
+        if isinstance(self.node_coords, int):
+            return False
+        if isinstance(self.center_coords, int):
+            return False
+        if isinstance(self.elem_conn, int):
+            return False
+        return True
+
+    def calculate_nodes(self):
+        """
+        Calculate the node coordinates, element connections and center coordinates
+        """
+        # Check that corners are set
+        if not self.are_corners_set():
+            abort("Corners have not been set in the mesh object")
+
+        # -- calculate node coordinates
+        self.calculate_node_coords()
+
+        # -- calculate element connections
+        self.calculate_elem_conn()
+
+        self.center_coords = np.stack(
+            [
+                self.center_lon2d.T.reshape((-1,)).T,
+                self.center_lat2d.T.reshape((-1,)).T,
+            ],
+            axis=1,
+        )
+
+    def create_esmf(self, mesh_fname):
+        """
+        Create an ESMF mesh file for the mesh
+
+        Parameters
+        ----------
+        mesh_fname : str
+            The path to write the ESMF meshfile
+
+        """
+        if not self.are_nodes_set():
+            abort("Nodes  have not been set in the mesh object")
+
+        # create output Xarray dataset
+        ds_out = xr.Dataset()
+
+        # origGridDims is only output if a mesh file was converted
+        if isinstance(self.center_lon2d, np.ndarray):
+            ds_out["origGridDims"] = xr.DataArray(
+                np.array(self.center_lon2d.shape, dtype=np.int32), dims=("origGridRank")
+            )
+
+        ds_out["nodeCoords"] = xr.DataArray(
+            self.node_coords, dims=("nodeCount", "coordDim"), attrs={"units": self.unit}
+        )
+        ds_out["elementConn"] = xr.DataArray(
+            self.elem_conn,
+            dims=("elementCount", "maxNodePElement"),
+            attrs={
+                "long_name": "Node indices that define the element connectivity",
+                "_FillValue": -1,
+            },
+        )
+        ds_out.elementConn.encoding = {"dtype": np.int32}
+
+        if isinstance(self.center_lon2d, np.ndarray):
+            ds_out["numElementConn"] = xr.DataArray(
+                4 * np.ones(self.center_lon2d.size, dtype=np.int32),
+                dims=("elementCount"),
+                attrs={"long_name": "Number of nodes per element"},
+            )
+        else:
+            ds_out["numElementConn"] = xr.DataArray(
+                self.num_elem_conn,
+                dims=("elementCount"),
+                attrs={"long_name": "Number of nodes per element"},
+            )
+
+        ds_out["centerCoords"] = xr.DataArray(
+            self.center_coords,
+            dims=("elementCount", "coordDim"),
+            attrs={"units": self.unit},
+        )
+
+        # -- add mask
+        ds_out["elementMask"] = xr.DataArray(
+            self.mask.T.reshape((-1,)).T,
+            dims=("elementCount"),
+            attrs={"units": "unitless", "_FillValue": -9999.0},
+        )
+        ds_out.elementMask.encoding = {"dtype": np.int32}
+
+        # -- add area if provided
+        if self.area is not None:
+            ds_out["elementArea"] = xr.DataArray(
+                self.area.T.reshape((-1,)).T,
+                dims=("elementCount"),
+                attrs={"units": "radians^2", "long_name": "area weights"},
+            )
+
+        # -- force no '_FillValue' if not specified (default Nan)
+        for var in ds_out.variables:
+            if "_FillValue" not in ds_out[var].encoding:
+                ds_out[var].encoding["_FillValue"] = None
+
+        # -- add global attributes
+        if self.mesh_name:
+            ds_out.attrs["title"] = "ESMF unstructured grid file  " + self.mesh_name
+        else:
+            ds_out.attrs["title"] = "ESMF unstructured grid file  "
+        ds_out.attrs["gridType"] = "unstructured mesh"
+        ds_out.attrs["version"] = "0.9"
+        ds_out.attrs["conventions"] = "ESMFMESH"
+        ds_out.attrs["date_created"] = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
+
+        # -- write Xarray dataset to file
+        if mesh_fname is not None:
+            logger.info("Writing ESMF Mesh file to : %s", mesh_fname)
+            ds_out.to_netcdf(mesh_fname)
+            logger.info("Successfully created ESMF Mesh file : %s", mesh_fname)
diff --git a/python/ctsm/site_and_regional/regional_case.py b/python/ctsm/site_and_regional/regional_case.py
index 2a2dc66bce..2f44b564ea 100644
--- a/python/ctsm/site_and_regional/regional_case.py
+++ b/python/ctsm/site_and_regional/regional_case.py
@@ -5,13 +5,18 @@
 # -- Import Python Standard Libraries
 import logging
 import os
+import argparse
+from datetime import datetime
 
 # -- 3rd party libraries
 import numpy as np
+import xarray as xr
 
 # -- import local classes for this script
 from ctsm.site_and_regional.base_case import BaseCase, USRDAT_DIR
+from ctsm.site_and_regional.mesh_type import MeshType
 from ctsm.utils import add_tag_to_filename
+from ctsm.utils import abort
 
 logger = logging.getLogger(__name__)
 
@@ -35,6 +40,8 @@ class RegionalCase(BaseCase):
         Region's name
     create_domain : bool
         flag for creating domain file
+    create_mesh : bool
+        flag for creating mesh file
     create_surfdata : bool
         flag for creating surface dataset
     create_landuse : bool
@@ -54,17 +61,34 @@ class RegionalCase(BaseCase):
         region's name or a combination of bounds of this
         region lat1-lat2_lon1-lon2
 
+    check_region_bounds
+        Check for the regional bounds
+
+    check_region_lons
+        Check for the regional lons
+
+    check_region_lats
+        Check for the regional lats
+
     create_domain_at_reg
         Create domain file at this region
 
     create_surfdata_at_reg
         Create surface dataset at this region
 
+    extract_mesh_at_reg
+        Extract mesh from the domain dataset created by create_domain_at_reg
+
     create_landuse_at_reg
         Create landuse file at this region
 
+    write_shell_commands(namelist)
+        write out xml commands to a file for usermods (i.e. shell_commands) for regional settings.
     """
 
+    # pylint: disable=too-many-instance-attributes
+    # the ones we have are useful
+
     def __init__(
         self,
         lat1,
@@ -77,6 +101,7 @@ def __init__(
         create_landuse,
         create_datm,
         create_user_mods,
+        create_mesh,
         out_dir,
         overwrite,
     ):
@@ -96,8 +121,13 @@ def __init__(
         self.lon1 = lon1
         self.lon2 = lon2
         self.reg_name = reg_name
+        self.create_mesh = create_mesh
+        self.mesh = None
         self.out_dir = out_dir
+        self.check_region_bounds()
         self.create_tag()
+        self.ni = None
+        self.nj = None
 
     def create_tag(self):
         """
@@ -112,6 +142,49 @@ def create_tag(self):
                 str(self.lon1), str(self.lon2), str(self.lat1), str(self.lat2)
             )
 
+    def check_region_bounds(self):
+        """
+        Check for the regional bounds
+        """
+        self.check_region_lons()
+        self.check_region_lats()
+
+    def check_region_lons(self):
+        """
+        Check for the regional lon bounds
+        """
+        if self.lon1 >= self.lon2:
+            err_msg = """
+            \n
+            ERROR: lon1 is bigger than lon2.
+            lon1 points to the westernmost longitude of the region. {}
+            lon2 points to the easternmost longitude of the region. {}
+            Please make sure lon1 is smaller than lon2.
+
+            Please note that if longitude in -180-0, the code automatically
+            convert it to 0-360.
+            """.format(
+                self.lon1, self.lon2
+            )
+            raise argparse.ArgumentTypeError(err_msg)
+
+    def check_region_lats(self):
+        """
+        Check for the regional lat bound
+        """
+        if self.lat1 >= self.lat2:
+            err_msg = """
+            \n
+            ERROR: lat1 is bigger than lat2.
+            lat1 points to the westernmost longitude of the region. {}
+            lat2 points to the easternmost longitude of the region. {}
+            Please make sure lat1 is smaller than lat2.
+
+            """.format(
+                self.lat1, self.lat2
+            )
+            raise argparse.ArgumentTypeError(err_msg)
+
     def create_domain_at_reg(self, indir, file):
         """
         Create domain file for this RegionalCase class.
@@ -145,6 +218,16 @@ def create_domain_at_reg(self, indir, file):
         f_in.close()
         f_out.close()
 
+        if self.create_mesh:
+            mesh_out = os.path.join(
+                self.out_dir,
+                os.path.splitext(fdomain_out)[0] + "_ESMF_UNSTRUCTURED_MESH.nc",
+            )
+            self.mesh = mesh_out
+            logger.info("creating mesh file from domain file: %s", wfile)
+            ds = xr.open_dataset(wfile, mask_and_scale=False, decode_times=False).transpose()
+            self.extract_mesh_at_reg(ds)
+
     def create_surfdata_at_reg(self, indir, file, user_mods_dir, specify_fsurf_out):
         """
         Create surface data file for this RegionalCase class.
@@ -189,6 +272,30 @@ def create_surfdata_at_reg(self, indir, file, user_mods_dir, specify_fsurf_out):
                 line = "fsurdat = '${}'".format(os.path.join(USRDAT_DIR, fsurf_out))
                 self.write_to_file(line, nl_clm)
 
+    def extract_mesh_at_reg(self, ds_in):
+        """
+        Create Mesh from Surface dataset netcdf file.
+        """
+        logger.info("Creating meshfile for  at region: %s", self.tag)
+
+        lat_name = "yc"
+        lon_name = "xc"
+
+        lats = ds_in[lat_name].astype(np.float32)
+        lons = ds_in[lon_name].astype(np.float32)
+
+        self.ni = len(lats.ni)
+        self.nj = len(lats.nj)
+
+        mask = ds_in["mask"].astype(np.float32)
+        if mask.max() > 1.0 or mask.min() < 0.0:
+            abort("Mask variable is not within 0 to 1")
+
+        this_mesh = MeshType(lats, lons, mask=mask)
+        this_mesh.calculate_corners()
+        this_mesh.calculate_nodes()
+        this_mesh.create_esmf(self.mesh)
+
     def create_landuse_at_reg(self, indir, file, user_mods_dir):
         """
         Create land use data file for this RegionalCase class.
@@ -228,3 +335,224 @@ def create_landuse_at_reg(self, indir, file, user_mods_dir):
                 # line = "landuse = '${}'".format(os.path.join(USRDAT_DIR, fluse_out))
                 line = "flanduse_timeseries = '${}'".format(os.path.join(USRDAT_DIR, fluse_out))
                 self.write_to_file(line, nl_clm)
+
+    def create_mesh_at_reg(self, mesh_dir, mesh_surf):
+        """
+        Create a mesh subsetted for the RegionalCase class.
+        """
+        logger.info("----------------------------------------------------------------------")
+        logger.info("Subsetting mesh file for region: %s", self.tag)
+
+        today = datetime.today()
+        today_string = today.strftime("%y%m%d")
+
+        mesh_in = os.path.join(mesh_dir, mesh_surf)
+        mesh_out = os.path.join(
+            self.out_dir,
+            os.path.splitext(mesh_surf)[0] + "_" + self.tag + "_c" + today_string + ".nc",
+        )
+
+        logger.info("mesh_in  :  %s", mesh_in)
+        logger.info("mesh_out :  %s", mesh_out)
+
+        self.mesh = mesh_out
+
+        node_coords, subset_element, subset_node, conn_dict = self.subset_mesh_at_reg(mesh_in)
+
+        f_in = xr.open_dataset(mesh_in)
+        self.write_mesh(f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out)
+
+    def subset_mesh_at_reg(self, mesh_in):
+        """
+        This function subsets the mesh based on lat and lon bounds given by RegionalCase class.
+        """
+        f_in = xr.open_dataset(mesh_in)
+        elem_count = len(f_in["elementCount"])
+        elem_conn = f_in["elementConn"]
+        num_elem_conn = f_in["numElementConn"]
+        node_count = len(f_in["nodeCount"])
+        node_coords = f_in["nodeCoords"]
+
+        subset_element = []
+        cnt = 0
+
+        for n in range(elem_count):
+            endx = elem_conn[n, : num_elem_conn[n].values].values
+            endx[
+                :,
+            ] -= 1  # convert to zero based index
+            endx = [int(xi) for xi in endx]
+
+            nlon = node_coords[endx, 0].values
+            nlat = node_coords[endx, 1].values
+
+            l1 = np.logical_or(nlon <= self.lon1, nlon >= self.lon2)
+            l2 = np.logical_or(nlat <= self.lat1, nlat >= self.lat2)
+
+            if np.any(np.logical_or(l1, l2)):
+                pass
+            else:
+                subset_element.append(n)
+                cnt += 1
+
+        subset_node = []
+        conn_dict = {}
+        cnt = 1
+        for n in range(node_count):
+            nlon = node_coords[n, 0].values
+            nlat = node_coords[n, 1].values
+
+            l1 = np.logical_or(nlon <= self.lon1, nlon >= self.lon2)
+            l2 = np.logical_or(nlat <= self.lat1, nlat >= self.lat2)
+
+            if np.logical_or(l1, l2):
+                conn_dict[n + 1] = -9999
+            else:
+                subset_node.append(n)
+                conn_dict[n + 1] = cnt
+                cnt += 1
+
+            # -- reverse logic
+            # l1 = np.logical_and(nlon >= self.lon1,nlon <= self.lon2)
+            # l2 = np.logical_and(nlat >= self.lat1,nlat <= self.lat2)
+            # if np.any(l1) and np.any(l2):
+            #    subset_node.append(n)
+            #    conn_dict[n+1] = cnt
+            #    cnt+=1
+            # else:
+            #    conn_dict[n+1] = -9999
+
+        return node_coords, subset_element, subset_node, conn_dict
+
+    @staticmethod
+    def write_mesh(f_in, node_coords, subset_element, subset_node, conn_dict, mesh_out):
+        # pylint: disable=unused-argument
+        """
+        This function writes out the subsetted mesh file.
+        """
+        corner_pairs = f_in.variables["nodeCoords"][
+            subset_node,
+        ]
+        variables = f_in.variables
+        global_attributes = f_in.attrs
+
+        max_node_dim = len(f_in["maxNodePElement"])
+
+        elem_count = len(subset_element)
+        elem_conn_out = np.empty(shape=[elem_count, max_node_dim])
+        elem_conn_index = f_in.variables["elementConn"][
+            subset_element,
+        ]
+
+        for n in range(elem_count):
+            for m in range(max_node_dim):
+                ndx = int(elem_conn_index[n, m])
+                elem_conn_out[n, m] = conn_dict[ndx]
+
+        num_elem_conn_out = np.empty(
+            shape=[
+                elem_count,
+            ]
+        )
+        num_elem_conn_out[:] = f_in.variables["numElementConn"][
+            subset_element,
+        ]
+
+        center_coords_out = np.empty(shape=[elem_count, 2])
+        center_coords_out[:, :] = f_in.variables["centerCoords"][subset_element, :]
+
+        if "elementMask" in variables:
+            elem_mask_out = np.empty(
+                shape=[
+                    elem_count,
+                ]
+            )
+            elem_mask_out[:] = f_in.variables["elementMask"][
+                subset_element,
+            ]
+
+        if "elementArea" in variables:
+            elem_area_out = np.empty(
+                shape=[
+                    elem_count,
+                ]
+            )
+            elem_area_out[:] = f_in.variables["elementArea"][
+                subset_element,
+            ]
+
+        # -- create output dataset
+        f_out = xr.Dataset()
+
+        f_out["nodeCoords"] = xr.DataArray(
+            corner_pairs, dims=("nodeCount", "coordDim"), attrs={"units": "degrees"}
+        )
+
+        f_out["elementConn"] = xr.DataArray(
+            elem_conn_out,
+            dims=("elementCount", "maxNodePElement"),
+            attrs={"long_name": "Node indices that define the element connectivity"},
+        )
+        f_out.elementConn.encoding = {"dtype": np.int32}
+
+        f_out["numElementConn"] = xr.DataArray(
+            num_elem_conn_out,
+            dims=("elementCount"),
+            attrs={"long_name": "Number of nodes per element"},
+        )
+        f_out.numElementConn.encoding = {"dtype": np.int32}
+
+        f_out["centerCoords"] = xr.DataArray(
+            center_coords_out,
+            dims=("elementCount", "coordDim"),
+            attrs={"units": "degrees"},
+        )
+
+        # -- add mask
+        if "elementMask" in variables:
+            f_out["elementMask"] = xr.DataArray(
+                elem_mask_out, dims=("elementCount"), attrs={"units": "unitless"}
+            )
+            f_out.elementMask.encoding = {"dtype": np.int32}
+
+        if "elementArea" in variables:
+            f_out["elementArea"] = xr.DataArray(
+                elem_area_out, dims=("elementCount"), attrs={"units": "unitless"}
+            )
+
+        # -- setting fill values
+        for var in variables:
+            if "_FillValue" in f_in[var].encoding:
+                f_out[var].encoding["_FillValue"] = f_in[var].encoding["_FillValue"]
+            else:
+                f_out[var].encoding["_FillValue"] = None
+
+        # -- add global attributes
+        for attr in global_attributes:
+            if attr != "timeGenerated":
+                f_out.attrs[attr] = global_attributes[attr]
+
+        f_out.attrs = {
+            "title": "ESMF unstructured grid file for a region",
+            "created_by": "subset_data",
+            "date_created": "{}".format(datetime.now()),
+        }
+
+        f_out.to_netcdf(mesh_out)
+        logger.info("Successfully created file (mesh_out) %s", mesh_out)
+
+    def write_shell_commands(self, namelist):
+        """
+        writes out xml commands commands to a file (i.e. shell_commands) for single-point runs
+        """
+        # write_to_file surrounds text with newlines
+        with open(namelist, "w") as nl_file:
+            self.write_to_file("# Change below line if you move the subset data directory", nl_file)
+            self.write_to_file("./xmlchange {}={}".format(USRDAT_DIR, self.out_dir), nl_file)
+            self.write_to_file("./xmlchange ATM_DOMAIN_MESH={}".format(str(self.mesh)), nl_file)
+            self.write_to_file("./xmlchange LND_DOMAIN_MESH={}".format(str(self.mesh)), nl_file)
+            self.write_to_file("./xmlchange MASK_MESH={}".format(str(str(self.mesh))), nl_file)
+            self.write_to_file("./xmlchange ATM_NX={}".format(str(str(self.ni))), nl_file)
+            self.write_to_file("./xmlchange LND_NX={}".format(str(str(self.ni))), nl_file)
+            self.write_to_file("./xmlchange ATM_NY={}".format(str(str(self.nj))), nl_file)
+            self.write_to_file("./xmlchange LND_NY={}".format(str(str(self.nj))), nl_file)
diff --git a/python/ctsm/subset_data.py b/python/ctsm/subset_data.py
index 3c8f6558a3..bc9ce97962 100644
--- a/python/ctsm/subset_data.py
+++ b/python/ctsm/subset_data.py
@@ -197,7 +197,7 @@ def get_parser():
     # -- region-specific parser options
     rg_parser.add_argument(
         "--lat1",
-        help="Region start latitude. [default: %(default)s]",
+        help="Region southernmost latitude. [default: %(default)s]",
         action="store",
         dest="lat1",
         required=False,
@@ -206,7 +206,7 @@ def get_parser():
     )
     rg_parser.add_argument(
         "--lat2",
-        help="Region end latitude. [default: %(default)s]",
+        help="Region northernmost latitude. [default: %(default)s]",
         action="store",
         dest="lat2",
         required=False,
@@ -215,7 +215,7 @@ def get_parser():
     )
     rg_parser.add_argument(
         "--lon1",
-        help="Region start longitude. [default: %(default)s]",
+        help="Region westernmost longitude. [default: %(default)s]",
         action="store",
         dest="lon1",
         required=False,
@@ -224,7 +224,7 @@ def get_parser():
     )
     rg_parser.add_argument(
         "--lon2",
-        help="Region end longitude. [default: %(default)s]",
+        help="Region easternmost longitude. [default: %(default)s]",
         action="store",
         dest="lon2",
         required=False,
@@ -250,14 +250,6 @@ def get_parser():
 
     # -- common options between both subparsers
     for subparser in [pt_parser, rg_parser]:
-        subparser.add_argument(
-            "--create-domain",
-            help="Create CLM domain file at single point/region. \
-            Domain files are not needed for NUOPC cases.",
-            action="store_true",
-            dest="create_domain",
-            required=False,
-        )
         subparser.add_argument(
             "--create-surface",
             help="Create surface data file at single point/region.",
@@ -265,20 +257,39 @@ def get_parser():
             dest="create_surfdata",
             required=False,
         )
+        subparser.add_argument(
+            "--surf-year",
+            help="Year for surface data file at single point/region \
+            (and start year for land-use timeseries).",
+            action="store",
+            dest="surf_year",
+            type=int,
+            default=2000,
+            required=False,
+        )
         subparser.add_argument(
             "--create-landuse",
-            help="Create landuse data file at single point/region.",
+            help="Create landuse data file at a single point/region.",
             action="store_true",
             dest="create_landuse",
             required=False,
         )
         subparser.add_argument(
             "--create-datm",
-            help="Create DATM forcing data at single point.",
+            help="Create DATM forcing data at a single point/region.",
             action="store_true",
             dest="create_datm",
             required=False,
         )
+        subparser.add_argument(
+            "--create-domain",
+            help="Create CLM domain file for a single point/region \
+            Domain files are not needed for NUOPC cases, \
+            but are needed to create mesh files that are needed for NUOPC cases.",
+            action="store_true",
+            dest="create_domain",
+            required=False,
+        )
         subparser.add_argument(
             "--create-user-mods",
             help="Create user mods directories and files for running CTSM with the subset data.",
@@ -394,19 +405,12 @@ def check_args(args):
         )
         raise argparse.ArgumentError(None, err_msg)
 
-    if not any(
-        [
-            args.create_surfdata,
-            args.create_domain,
-            args.create_landuse,
-            args.create_datm,
-        ]
-    ):
+    if not any([args.create_surfdata, args.create_landuse, args.create_datm, args.create_domain]):
         err_msg = textwrap.dedent(
             """\
                 \n ------------------------------------
                 \n Must supply one of:
-                \n --create-surface \n --create-landuse \n --create-datm \n --create-domain \n
+                \n --create-surface \n --create-landuse \n --create-datm \n --create-domain \n \n
                 """
         )
         raise argparse.ArgumentError(None, err_msg)
@@ -429,6 +433,43 @@ def check_args(args):
         )
         raise argparse.ArgumentError(None, err_msg)
 
+    if args.create_landuse and not args.create_surfdata:
+        err_msg = textwrap.dedent(
+            """\
+                \n ------------------------------------
+                \n --create-landuse option requires the --create-surface option:
+                """
+        )
+        raise argparse.ArgumentError(None, err_msg)
+    if args.surf_year != 2000 and not args.create_surfdata:
+        err_msg = textwrap.dedent(
+            """\
+                \n ------------------------------------
+                \n --surf-year option is set to something besides the default of 2000
+                \n without the --create-surface option"
+                """
+        )
+        raise argparse.ArgumentError(None, err_msg)
+
+    if args.surf_year != 1850 and args.create_landuse:
+        err_msg = textwrap.dedent(
+            """\
+                \n ------------------------------------
+                \n --surf-year option is NOT set to 1850 and the --create-landuse option
+                \n is selected which requires it to be 1850
+                """
+        )
+        raise argparse.ArgumentError(None, err_msg)
+
+    if args.surf_year != 1850 and args.surf_year != 2000:
+        err_msg = textwrap.dedent(
+            """\
+                \n ------------------------------------
+                \n --surf-year option can only be set to 1850 or 2000
+                """
+        )
+        raise argparse.ArgumentError(None, err_msg)
+
     if args.out_surface and os.path.exists(args.out_surface) and not args.overwrite:
         err_msg = textwrap.dedent(
             """\
@@ -438,6 +479,32 @@ def check_args(args):
         )
         raise argparse.ArgumentError(None, err_msg)
 
+    if args.run_type == "region" and args.create_user_mods:
+        if not args.create_mesh:
+            err_msg = textwrap.dedent(
+                """\
+                      \n ------------------------------------
+                      \nERROR: For regional cases, you can not create user_mods
+                      \nwithout creating the mesh file.
+
+                      \nPlease rerun the script adding --create-mesh to subset the mesh file."
+                      """
+            )
+            raise argparse.ArgumentError(None, err_msg)
+
+    if args.run_type == "region" and args.create_mesh:
+        if not args.create_domain:
+            err_msg = textwrap.dedent(
+                """\
+                      \n ------------------------------------
+                      \nERROR: For regional cases, you can not create mesh files
+                      \nwithout creating the domain file.
+
+                      \nPlease rerun the script adding --create-domain to subset the domain file."
+                      """
+            )
+            raise argparse.ArgumentError(None, err_msg)
+
 
 def setup_user_mods(user_mods_dir, cesmroot):
     """
@@ -530,6 +597,10 @@ def setup_files(args, defaults, cesmroot):
             clmforcingindir,
             os.path.join(defaults.get("surfdat", "dir")),
         ),
+        "mesh_dir": os.path.join(
+            clmforcingindir,
+            os.path.join(defaults.get("surfdat", "mesh_dir")),
+        ),
         "fluse_dir": os.path.join(
             clmforcingindir,
             os.path.join(defaults.get("landuse", "dir")),
@@ -537,6 +608,7 @@ def setup_files(args, defaults, cesmroot):
         "fsurf_in": fsurf_in,
         "fsurf_out": fsurf_out,
         "fluse_in": fluse_in,
+        "mesh_surf": defaults.get("surfdat", "mesh_surf"),
         "datm_tuple": DatmFiles(
             dir_input_datm,
             dir_output_datm,
@@ -591,10 +663,6 @@ def subset_point(args, file_dict: dict):
 
     logger.debug(single_point)
 
-    # --  Create CTSM domain file
-    if single_point.create_domain:
-        single_point.create_domain_at_point(file_dict["main_dir"], file_dict["fdomain_in"])
-
     # --  Create CTSM surface data file
     if single_point.create_surfdata:
         single_point.create_surfdata_at_point(
@@ -648,6 +716,7 @@ def subset_region(args, file_dict: dict):
         create_landuse=args.create_landuse,
         create_datm=args.create_datm,
         create_user_mods=args.create_user_mods,
+        create_mesh=args.create_mesh,
         out_dir=args.out_dir,
         overwrite=args.overwrite,
     )
@@ -667,12 +736,27 @@ def subset_region(args, file_dict: dict):
             specify_fsurf_out=file_dict["fsurf_out"],
         )
 
+    # if region.create_mesh:
+    #    region.create_mesh_at_reg (file_dict["mesh_dir"], file_dict["mesh_surf"])
+
     # --  Create CTSM transient landuse data file
     if region.create_landuse:
         region.create_landuse_at_reg(
             file_dict["fluse_dir"], file_dict["fluse_in"], args.user_mods_dir
         )
 
+    # -- Write shell commands
+    if region.create_user_mods:
+        region.write_shell_commands(os.path.join(args.user_mods_dir, "shell_commands"))
+
+        print("\nFor running this regional case with the created user_mods : ")
+        print(
+            "./create_newcase --case case --res CLM_USRDAT --compset I2000Clm51BgcCrop",
+            "--run-unsupported --user-mods-dirs ",
+            args.user_mods_dir,
+            "\n\n",
+        )
+
     logger.info("Successfully ran script for a regional case.")
 
 
diff --git a/python/ctsm/test/test_advanced_sys_mesh_plotter.py b/python/ctsm/test/test_advanced_sys_mesh_plotter.py
new file mode 100755
index 0000000000..090f7806bf
--- /dev/null
+++ b/python/ctsm/test/test_advanced_sys_mesh_plotter.py
@@ -0,0 +1,61 @@
+#!/usr/bin/env python3
+"""
+Advanced System tests for mesh_plotter (requires the ctsm_pylib_wdask conda environment)
+
+"""
+
+import unittest
+import os
+import sys
+import tempfile
+import shutil
+import glob
+
+# pylint: disable=wrong-import-position
+from ctsm.path_utils import path_to_ctsm_root
+from ctsm import unit_testing
+from ctsm.mesh_plotter import main
+
+# pylint: disable=invalid-name
+
+
+class SysTestMeshMaker(unittest.TestCase):
+    """
+    Basic class for testing mesh_plotter.py.
+    """
+
+    def setUp(self):
+        """Setup for all tests"""
+        testinputs_path = os.path.join(path_to_ctsm_root(), "python/ctsm/test/testinputs")
+        self._testinputs_path = testinputs_path
+        self._infile = os.path.join(
+            testinputs_path,
+            "ESMF_mesh_5x5pt_amazon_from_domain_c230308.nc",
+        )
+        self._tempdir = tempfile.mkdtemp()
+        self.mesh_out = self._tempdir + "/mesh_out"
+
+    def tearDown(self):
+        """
+        Remove temporary directory
+        """
+        shutil.rmtree(self._tempdir, ignore_errors=True)
+
+    def test_basic(self):
+        """Do a simple basic test"""
+        sys.argv = [
+            "mesh_plotter",
+            "--input",
+            self._infile,
+            "--output",
+            self.mesh_out,
+        ]
+        main()
+        plotfiles = glob.glob(self._tempdir + "/*.png")
+        if not plotfiles:
+            self.fail("plot files were NOT created as they should have")
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_advanced_unit_mesh_plotter.py b/python/ctsm/test/test_advanced_unit_mesh_plotter.py
new file mode 100755
index 0000000000..e4930ce232
--- /dev/null
+++ b/python/ctsm/test/test_advanced_unit_mesh_plotter.py
@@ -0,0 +1,62 @@
+#!/usr/bin/env python3
+"""
+Unit tests for mesh_plotter
+
+You can run this by:
+    python -m unittest test_advanced_unit_mesh_plotter.py
+"""
+
+import unittest
+import os
+import sys
+
+# -- 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
+from ctsm import unit_testing
+from ctsm.mesh_plotter import get_parser, process_and_check_args, main
+
+# pylint: disable=invalid-name
+
+
+class TestMeshMaker(unittest.TestCase):
+    """
+    Basic class for testing mesh_plotter.py.
+    """
+
+    def test_input_file_dne(self):
+        """
+        Test that exits if input file does not exist
+        """
+        sys.argv = [
+            "mesh_plotter",
+            "--input",
+            "zztop.nc",
+            "--outdir",
+            ".",
+        ]
+        with self.assertRaisesRegex(SystemExit, "Input file not found."):
+            main()
+
+    def test_default_outfile_as_expected(self):
+        """
+        Test that the default outfile is as expected
+        """
+        infile = "ctsm/test/testinputs/default_data.cfg"
+        sys.argv = [
+            "mesh_plotter",
+            "--input",
+            infile,
+        ]
+        parser = get_parser()
+        args = parser.parse_args()
+        args = process_and_check_args(args)
+        expected_outdir = os.path.join(os.getcwd(), "meshes")
+        self.assertEqual(args.out_dir, expected_outdir, "Default out_dir is not as expected")
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_sys_mesh_maker.py b/python/ctsm/test/test_sys_mesh_maker.py
new file mode 100755
index 0000000000..fb97a31445
--- /dev/null
+++ b/python/ctsm/test/test_sys_mesh_maker.py
@@ -0,0 +1,428 @@
+#!/usr/bin/env python3
+"""
+System tests for mesh_maker
+
+"""
+
+import unittest
+import os
+import sys
+import tempfile
+import shutil
+import xarray as xr
+
+# pylint: disable=wrong-import-position
+from ctsm.path_utils import path_to_ctsm_root
+from ctsm import unit_testing
+from ctsm.mesh_maker import main
+from ctsm.mesh_maker import read_main
+
+# pylint: disable=invalid-name
+
+
+class SysTestMeshMaker(unittest.TestCase):
+    """
+    Basic class for testing mesh_maker.py.
+    """
+
+    def setUp(self):
+        """Setup for all tests"""
+        testinputs_path = os.path.join(path_to_ctsm_root(), "python/ctsm/test/testinputs")
+        self._testinputs_path = testinputs_path
+        self._infile = os.path.join(
+            testinputs_path, "surfdata_5x5_amazon_hist_16pfts_CMIP6_2000_c231031_modified.nc"
+        )
+        self._tempdir = tempfile.mkdtemp()
+        self.mesh_out = os.path.join(self._tempdir, "mesh_out.nc")
+
+    def tearDown(self):
+        """
+        Remove temporary directory
+        """
+        shutil.rmtree(self._tempdir, ignore_errors=True)
+
+    def test_basic(self):
+        """Do a simple basic test"""
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            self._infile,
+            "--lat",
+            "LATIXY",
+            "--lon",
+            "LONGXY",
+            "--output",
+            self.mesh_out,
+        ]
+        main()
+
+    def test_region(self):
+        """Do a basic test for a small regional grid"""
+        infile = os.path.join(
+            self._testinputs_path,
+            "surfdata_5x5_amazon_hist_78pfts_CMIP6_2000_c230517_modified_with_crop.nc",
+        )
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            infile,
+            "--lat",
+            "LATIXY",
+            "--lon",
+            "LONGXY",
+            "--mask",
+            "LANDFRAC_PFT",
+            "--output",
+            self.mesh_out,
+        ]
+        main()
+
+    def compare_mesh_files(self, mesh_out, expected):
+        """Compare two mesh files that you expect to be equal"""
+        self.assertEqual(
+            mesh_out.dims["coordDim"], expected.dims["coordDim"], "coordDim not the same"
+        )
+        if "origGridRank" in mesh_out.variables and "origGridRank" in expected.variables:
+            self.assertEqual(
+                mesh_out.dims["origGridRank"],
+                expected.dims["origGridRank"],
+                "origGridRank not the same",
+            )
+            equalorigGridDims = mesh_out.origGridDims == expected.origGridDims
+            self.assertTrue(equalorigGridDims.all, "origGridDims different")
+            compare_files = True
+        else:
+            # don't compare files if origGridRank isn't on one of the files
+            compare_files = False
+
+        self.assertEqual(
+            mesh_out.dims["nodeCount"], expected.dims["nodeCount"], "nodeCount not the same"
+        )
+        self.assertEqual(
+            mesh_out.dims["elementCount"],
+            expected.dims["elementCount"],
+            "elementCount not the same",
+        )
+        self.assertEqual(
+            mesh_out.dims["maxNodePElement"],
+            expected.dims["maxNodePElement"],
+            "maxNodePElement not the same",
+        )
+        equalelementConn = mesh_out.elementConn == expected.elementConn
+        equalnumElementConn = mesh_out.numElementConn == expected.numElementConn
+        equalcenterCoords = mesh_out.centerCoords == expected.centerCoords
+        equalelementMask = mesh_out.elementMask == expected.elementMask
+
+        if "elementArea" in mesh_out.variables and "elementArea" in expected.variables:
+            equalelementArea = mesh_out.elementArea == expected.elementArea
+            self.assertTrue(equalelementArea.all, "area different")
+
+        self.assertTrue(equalelementConn.all, "elementConn different")
+        self.assertTrue(equalnumElementConn.all, "numElementConn different")
+        self.assertTrue(equalcenterCoords.all, "centerCoords different")
+        self.assertTrue(equalelementMask.all, "mask different")
+        if compare_files:
+            self.assertTrue(
+                mesh_out.equals(expected),
+                "Output mesh does not compare to the expected baseline file",
+            )
+
+    def test_domainfile_region_warea(self):
+        """
+        Do a basic test for a small regional grid with a domain file
+        rather than a surfdata file including area
+        """
+        infile = os.path.join(self._testinputs_path, "domain.lnd.5x5pt-amazon_navy.090715.nc")
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            infile,
+            "--lat",
+            "yc",
+            "--lon",
+            "xc",
+            "--mask",
+            "mask",
+            "--area",
+            "area",
+            "--output",
+            self.mesh_out,
+        ]
+        main()
+        expected_mesh = os.path.join(
+            self._testinputs_path, "ESMF_mesh_5x5pt_amazon_from_domain_c230308.nc"
+        )
+        mesh_out = xr.open_dataset(self.mesh_out)
+        expected = xr.open_dataset(expected_mesh)
+        self.compare_mesh_files(mesh_out, expected)
+
+    def test_domainfile_SAmerica_region_warea(self):
+        """
+        Do a basic test for a South America regional grid with a domain file
+        rather than a surfdata file including area
+        """
+        infile = os.path.join(
+            self._testinputs_path, "domain.lnd.fv0.9x1.25_gx1v7_f09_58x45_SouthAmerica_c230522.nc"
+        )
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            infile,
+            "--lat",
+            "yc",
+            "--lon",
+            "xc",
+            "--mask",
+            "mask",
+            "--area",
+            "area",
+            "--output",
+            self.mesh_out,
+        ]
+        main()
+        expected_mesh = os.path.join(
+            self._testinputs_path,
+            "ESMF_mesh_fv0.9x1.25_gx1v7_f09_58x45_SouthAmerica_from_domain_c230522.nc",
+        )
+        mesh_out = xr.open_dataset(self.mesh_out)
+        expected = xr.open_dataset(expected_mesh)
+        self.compare_mesh_files(mesh_out, expected)
+
+    def test_domainfile_f10_warea(self):
+        """
+        Do a test of converting the f10 domain file
+        """
+        infile = os.path.join(self._testinputs_path, "domain.lnd.fv10x15_gx3v7.180321.nc")
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            infile,
+            "--lat",
+            "yc",
+            "--lon",
+            "xc",
+            "--mask",
+            "mask",
+            "--area",
+            "area",
+            "--output",
+            self.mesh_out,
+        ]
+        with self.assertRaisesRegex(SystemExit, "Expected size for element connections is wrong"):
+            main()
+
+    def test_readfile(self):
+        """
+        Test that reading a file results in the same mesh as converting one
+        """
+        infile = os.path.join(
+            self._testinputs_path, "ESMF_mesh_5x5pt_amazon_from_domain_c230308.nc"
+        )
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            infile,
+            "--lat",
+            "yc",
+            "--lon",
+            "xc",
+            "--output",
+            self.mesh_out,
+        ]
+        read_main()
+        mesh_out = xr.open_dataset(self.mesh_out)
+        expected = xr.open_dataset(infile)
+        self.compare_mesh_files(mesh_out, expected)
+
+    def test_add_mask(self):
+        """Do a simple basic test also adding mask"""
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            self._infile,
+            "--lat",
+            "LATIXY",
+            "--lon",
+            "LONGXY",
+            "--mask",
+            "LANDFRAC_PFT",
+            "--output",
+            self.mesh_out,
+        ]
+        main()
+
+    def test_noinput(self):
+        """Test with an input file that does not exist"""
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            "zztop",
+            "--lat",
+            "LATIXY",
+            "--lon",
+            "LONGXY",
+            "--output",
+            self.mesh_out,
+        ]
+        with self.assertRaisesRegex(SystemExit, "Input file not found."):
+            main()
+
+    def test_singlepoint_dies(self):
+        """Test that a single point file dies because we don't need mesh files
+        for single point cases"""
+        infile = os.path.join(
+            self._testinputs_path,
+            "surfdata_1x1_mexicocityMEX_hist_16pfts_Irrig_CMIP6_simyr2000_c221206.nc",
+        )
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            infile,
+            "--lat",
+            "LATIXY",
+            "--lon",
+            "LONGXY",
+            "--output",
+            self.mesh_out,
+        ]
+        with self.assertRaisesRegex(
+            SystemExit, r"No need to create a mesh file for a single point grid."
+        ):
+            main()
+
+    def test_nolongs(self):
+        """Bad name for longitude"""
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            self._infile,
+            "--lat",
+            "LATIXY",
+            "--lon",
+            "zztop",
+            "--output",
+            self.mesh_out,
+        ]
+        with self.assertRaisesRegex(SystemExit, "Input file does not have variable named zztop"):
+            main()
+
+    def test_nolats(self):
+        """Bad name for latitude"""
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            self._infile,
+            "--lat",
+            "zztop",
+            "--lon",
+            "LONGXY",
+            "--output",
+            self.mesh_out,
+        ]
+        with self.assertRaisesRegex(SystemExit, "Input file does not have variable named zztop"):
+            main()
+
+    def test_badareaname(self):
+        """Bad name for area"""
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            self._infile,
+            "--lat",
+            "LATIXY",
+            "--lon",
+            "LONGXY",
+            "--area",
+            "zztop",
+            "--output",
+            self.mesh_out,
+        ]
+        with self.assertRaisesRegex(
+            SystemExit, "Input file does not have area variable named zztop"
+        ):
+            main()
+
+    def test_badmaskname(self):
+        """Bad name for mask"""
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            self._infile,
+            "--lat",
+            "LATIXY",
+            "--lon",
+            "LONGXY",
+            "--mask",
+            "zztop",
+            "--output",
+            self.mesh_out,
+        ]
+        with self.assertRaisesRegex(
+            SystemExit, "Input file does not have mask variable named zztop"
+        ):
+            main()
+
+    def test_badareaunits(self):
+        """Bad area units"""
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            self._infile,
+            "--lat",
+            "LATIXY",
+            "--lon",
+            "LONGXY",
+            "--area",
+            "PCT_CROP",
+            "--output",
+            self.mesh_out,
+        ]
+        with self.assertRaisesRegex(
+            SystemExit, r"Area does NOT have the correct units of radians\^2 but has unitless"
+        ):
+            main()
+
+    def test_missingreaunits(self):
+        """Missing area units"""
+        self._infile = os.path.join(
+            self._testinputs_path,
+            "surfdata_5x5_amazon_hist_78pfts_CMIP6_2000_c230517_modified_with_crop.nc",
+        )
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            self._infile,
+            "--lat",
+            "LATIXY",
+            "--lon",
+            "LONGXY",
+            "--area",
+            "PCT_CROP",
+            "--output",
+            self.mesh_out,
+        ]
+        with self.assertRaisesRegex(SystemExit, r"Units attribute is NOT on the area variable"):
+            main()
+
+    def test_badmaskvalues(self):
+        """Bad mask values"""
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            self._infile,
+            "--lat",
+            "LATIXY",
+            "--lon",
+            "LONGXY",
+            "--mask",
+            "LAKEDEPTH",
+            "--output",
+            self.mesh_out,
+        ]
+        with self.assertRaisesRegex(SystemExit, "Mask variable is not within 0 to 1"):
+            main()
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_unit_mesh_maker.py b/python/ctsm/test/test_unit_mesh_maker.py
new file mode 100755
index 0000000000..f474e02616
--- /dev/null
+++ b/python/ctsm/test/test_unit_mesh_maker.py
@@ -0,0 +1,113 @@
+#!/usr/bin/env python3
+"""
+Unit tests for mesh_maker
+
+You can run this by:
+    python -m unittest test_unit_mesh_maker.py
+"""
+
+import unittest
+import os
+import sys
+
+# -- 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
+from ctsm import unit_testing
+from ctsm.mesh_maker import get_parser, process_and_check_args, main
+
+# pylint: disable=invalid-name
+
+
+class TestMeshMaker(unittest.TestCase):
+    """
+    Basic class for testing mesh_maker.py.
+    """
+
+    def test_input_file_dne(self):
+        """
+        Test that exits if input file does not exist
+        """
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            "zztop.nc",
+            "--lat",
+            "lsmlat",
+            "--lon",
+            "lsmlon",
+            "--outdir",
+            ".",
+        ]
+        with self.assertRaisesRegex(SystemExit, "Input file not found."):
+            main()
+
+    def test_outfile_and_outdir(self):
+        """
+        Test that exits if both outfile and outdir are provided
+        """
+        infile = "ctsm/test/testinputs/default_data.cfg"
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            infile,
+            "--lat",
+            "lsmlat",
+            "--lon",
+            "lsmlon",
+            "--outdir",
+            ".",
+            "--output",
+            "outthing.nc",
+        ]
+        with self.assertRaisesRegex(SystemExit, "You have provided both --outdir and --output."):
+            main()
+
+    def test_outfile_exists_and_no_overwrite(self):
+        """
+        Test that exits if outfile already exists and overwrite was not used
+        """
+        infile = "ctsm/test/testinputs/default_data.cfg"
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            infile,
+            "--lat",
+            "lsmlat",
+            "--lon",
+            "lsmlon",
+            "--output",
+            infile,
+        ]
+        with self.assertRaisesRegex(
+            SystemExit,
+            "output meshfile exists, please choose --overwrite to overwrite the mesh file.",
+        ):
+            main()
+
+    def test_default_outfile_as_expected(self):
+        """
+        Test that the default outfile is as expected
+        """
+        infile = "ctsm/test/testinputs/default_data.cfg"
+        sys.argv = [
+            "mesh_maker",
+            "--input",
+            infile,
+            "--lat",
+            "lsmlat",
+            "--lon",
+            "lsmlon",
+        ]
+        parser = get_parser()
+        args = parser.parse_args()
+        args = process_and_check_args(args)
+        expected_outdir = os.path.join(os.getcwd(), "meshes")
+        self.assertEqual(args.out_dir, expected_outdir, "Default out_dir is not as expected")
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_unit_mesh_type.py b/python/ctsm/test/test_unit_mesh_type.py
new file mode 100644
index 0000000000..ba59bcefe8
--- /dev/null
+++ b/python/ctsm/test/test_unit_mesh_type.py
@@ -0,0 +1,75 @@
+#!/usr/bin/env python3
+"""
+Unit tests for mesh_type
+
+You can run this by:
+    python -m unittest test_unit_mesh_type.py
+"""
+
+import os
+import sys
+import unittest
+import numpy as np
+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
+from ctsm import unit_testing
+
+# from ctsm.site_and_regional.mesh_type import MeshType
+from ctsm.site_and_regional.mesh_type import MeshType
+
+# Allow test names that pylint doesn't like; otherwise hard to make them
+# readable
+# pylint: disable=invalid-name
+
+
+class TestMeshType(unittest.TestCase):
+    """
+    Basic class for testing mesh_type.py.
+    """
+
+    def setUp(self):
+        """Setup for all tests"""
+        lon0 = np.array([120.0])
+        lat0 = np.array([45.0])
+        x_dim = "lon"
+        y_dim = "lat"
+        self.lons = xr.DataArray(lon0, name="lon", dims=x_dim, coords={x_dim: lon0})
+        self.lats = xr.DataArray(lat0, name="lat", dims=y_dim, coords={y_dim: lat0})
+
+        self.mesh = MeshType(self.lats, self.lons)
+
+    def test_read_file_fails_notXarrayDataset(self):
+        """Test that read_file properly fails if not given an X-array Dataset"""
+        with self.assertRaisesRegex(SystemExit, "Input file is not a X-Array DataSet type"):
+            self.mesh.read_file(1.0)
+
+    def test_read_file_fails_badvarlist(self):
+        """Test that read_file properly fails if input dataset does not have required variables"""
+        ds = xr.Dataset(
+            {
+                "PCT_NATVEG": xr.DataArray(
+                    data=np.random.rand(1, 1),
+                    dims=["lat", "lon"],
+                    coords={"lat": self.lats, "lon": self.lons},
+                    attrs={
+                        "long_name": "total percent natural vegetation landunit",
+                        "units": "unitless",
+                    },
+                ),
+            },
+            attrs={"Conventions": "test data only"},
+        )
+        with self.assertRaisesRegex(
+            SystemExit, "Variable expected to be on an ESMF mesh file is NOT found: nodeCoords"
+        ):
+            self.mesh.read_file(ds)
+
+
+if __name__ == "__main__":
+    unit_testing.setup_for_tests()
+    unittest.main()
diff --git a/python/ctsm/test/test_unit_subset_data.py b/python/ctsm/test/test_unit_subset_data.py
index 1994a29102..143b632e8c 100755
--- a/python/ctsm/test/test_unit_subset_data.py
+++ b/python/ctsm/test/test_unit_subset_data.py
@@ -120,7 +120,7 @@ def test_check_args_outsurfdat_fails_without_create_surface(self):
         Test that check args does not allow an output surface dataset to be specified
         when create-surface is not on
         """
-        sys.argv = ["subset_data", "point", "--create-landuse", "--out-surface", "outputsurface.nc"]
+        sys.argv = ["subset_data", "point", "--create-datm", "--out-surface", "outputsurface.nc"]
         self.args = self.parser.parse_args()
         with self.assertRaisesRegex(
             argparse.ArgumentError,
@@ -128,6 +128,56 @@ def test_check_args_outsurfdat_fails_without_create_surface(self):
         ):
             check_args(self.args)
 
+    def test_check_args_fails_for_timeseries_without_correct_surface_year(self):
+        """
+        Test that check args does not allow landuse-timeseries to be used
+        without providing the correct start surface year
+        """
+        sys.argv = ["subset_data", "point", "--create-landuse", "--create-surface"]
+        self.args = self.parser.parse_args()
+        with self.assertRaisesRegex(
+            argparse.ArgumentError,
+            "--surf-year option is NOT set to 1850 and the --create-landuse option",
+        ):
+            check_args(self.args)
+
+    def test_check_args_fails_for_surf_year_without_surface(self):
+        """
+        Test that check args does not allow surf_year to be set
+        without the create-surface option
+        """
+        sys.argv = ["subset_data", "point", "--create-datm", "--surf-year", "1850"]
+        self.args = self.parser.parse_args()
+        with self.assertRaisesRegex(
+            argparse.ArgumentError,
+            "--surf-year option is set to something besides the default of 2000",
+        ):
+            check_args(self.args)
+
+    def test_check_args_fails_for_landuse_without_surface(self):
+        """
+        Test that check args does not allow landuse to be set
+        without the create-surface option
+        """
+        sys.argv = ["subset_data", "point", "--create-landuse"]
+        self.args = self.parser.parse_args()
+        with self.assertRaisesRegex(
+            argparse.ArgumentError,
+            "--create-landuse option requires the --create-surface option",
+        ):
+            check_args(self.args)
+
+    def test_check_args_fails_bad_surface_year(self):
+        """
+        Test that check args does not allow --surf-year to be bad
+        """
+        sys.argv = ["subset_data", "point", "--create-surface", "--surf-year", "2305"]
+        self.args = self.parser.parse_args()
+        with self.assertRaisesRegex(
+            argparse.ArgumentError, "--surf-year option can only be set to 1850 or 2000"
+        ):
+            check_args(self.args)
+
     def test_check_args_outsurfdat_fails_without_overwrite(self):
         """
         Test that check args does not allow an output surface dataset to be specified
@@ -157,6 +207,58 @@ def test_inputdata_setup_files_bad_inputdata_arg(self):
         with self.assertRaisesRegex(SystemExit, "inputdata directory does not exist"):
             setup_files(self.args, self.defaults, self.cesmroot)
 
+    def test_create_user_mods_without_create_mesh(self):
+        """
+        Test that you can't run create user mods without also doing create_mesh
+        """
+        sys.argv = ["subset_data", "region", "--create-user-mods", "--create-surface"]
+        self.args = self.parser.parse_args()
+        with self.assertRaisesRegex(
+            argparse.ArgumentError, "For regional cases, you can not create user_mods"
+        ):
+            check_args(self.args)
+
+    def test_create_mesh_without_domain(self):
+        """
+        Test that you can't run create mesh without domain
+        """
+        sys.argv = [
+            "subset_data",
+            "region",
+            "--create-user-mods",
+            "--create-surface",
+            "--create-mesh",
+        ]
+        self.args = self.parser.parse_args()
+        with self.assertRaisesRegex(
+            argparse.ArgumentError, "For regional cases, you can not create mesh files"
+        ):
+            check_args(self.args)
+
+    def test_complex_option_works(self):
+        """
+        Test that check_args won't flag a set of complex options that is valid
+        Do user-mods, surface and landuse-timeseries, as well as DATM, for verbose with crop
+        """
+        sys.argv = [
+            "subset_data",
+            "region",
+            "--reg",
+            "testname",
+            "--create-user-mods",
+            "--create-surface",
+            "--create-landuse",
+            "--surf-year",
+            "1850",
+            "--create-mesh",
+            "--create-domain",
+            "--create-datm",
+            "--verbose",
+            "--crop",
+        ]
+        self.args = self.parser.parse_args()
+        check_args(self.args)
+
 
 if __name__ == "__main__":
     unit_testing.setup_for_tests()
diff --git a/python/ctsm/test/testinputs/ESMF_mesh_5x5pt_amazon_from_domain_c230308.nc b/python/ctsm/test/testinputs/ESMF_mesh_5x5pt_amazon_from_domain_c230308.nc
new file mode 100644
index 0000000000..863d96c2c4
--- /dev/null
+++ b/python/ctsm/test/testinputs/ESMF_mesh_5x5pt_amazon_from_domain_c230308.nc
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9f791e4a1ad27dc22f4779fbf74c6023a1624beaf3102ff750e023e4341761b2
+size 14920
diff --git a/python/ctsm/test/testinputs/ESMF_mesh_fv0.9x1.25_gx1v7_f09_58x45_SouthAmerica_from_domain_c230522.nc b/python/ctsm/test/testinputs/ESMF_mesh_fv0.9x1.25_gx1v7_f09_58x45_SouthAmerica_from_domain_c230522.nc
new file mode 100644
index 0000000000..5dff70e72e
--- /dev/null
+++ b/python/ctsm/test/testinputs/ESMF_mesh_fv0.9x1.25_gx1v7_f09_58x45_SouthAmerica_from_domain_c230522.nc
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:275b0915db6328d872c6a5ea452c1ac920bdbad19da37ce12f3ce4b5b1f4bab7
+size 173511
diff --git a/python/ctsm/test/testinputs/default_data.cfg b/python/ctsm/test/testinputs/default_data.cfg
new file mode 100644
index 0000000000..d736f2e933
--- /dev/null
+++ b/python/ctsm/test/testinputs/default_data.cfg
@@ -0,0 +1,30 @@
+[main]
+clmforcingindir = /glade/campaign/cesm/cesmdata/cseg/inputdata
+
+[datm_gswp3]
+dir = atm/datm7/atm_forcing.datm7.GSWP3.0.5d.v1.c170516
+domain = domain.lnd.360x720_gswp3.0v1.c170606.nc
+solardir = Solar
+precdir = Precip
+tpqwdir = TPHWL
+solartag = clmforc.GSWP3.c2011.0.5x0.5.Solr.
+prectag = clmforc.GSWP3.c2011.0.5x0.5.Prec.
+tpqwtag = clmforc.GSWP3.c2011.0.5x0.5.TPQWL.
+solarname = CLMGSWP3v1.Solar
+precname = CLMGSWP3v1.Precip
+tpqwname = CLMGSWP3v1.TPQW
+
+[surfdat]
+dir = lnd/clm2/surfdata_esmf/ctsm5.2.0
+surfdat_16pft = surfdata_0.9x1.25_hist_2000_16pfts_c240216.nc
+surfdat_78pft = surfdata_0.9x1.25_hist_2000_78pfts_c240216.nc
+mesh_dir = share/meshes/
+mesh_surf = fv0.9x1.25_141008_ESMFmesh.nc
+
+[landuse]
+dir = lnd/clm2/surfdata_esmf/ctsm5.2.0
+landuse_16pft = landuse.timeseries_0.9x1.25_SSP2-4.5_1850-2100_78pfts_c240216.nc
+landuse_78pft = landuse.timeseries_0.9x1.25_SSP2-4.5_1850-2100_78pfts_c240216.nc
+
+[domain]
+file = share/domains/domain.lnd.fv0.9x1.25_gx1v7.151020.nc
diff --git a/python/ctsm/test/testinputs/domain.lnd.5x5pt-amazon_navy.090715.nc b/python/ctsm/test/testinputs/domain.lnd.5x5pt-amazon_navy.090715.nc
new file mode 100644
index 0000000000..450536c2fc
--- /dev/null
+++ b/python/ctsm/test/testinputs/domain.lnd.5x5pt-amazon_navy.090715.nc
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:714824c2ded32b98a7db6b3c62b66267b2ecea8a9a0af61103bd6a816eab3e9d
+size 5136
diff --git a/python/ctsm/test/testinputs/domain.lnd.fv0.9x1.25_gx1v7_f09_58x45_SouthAmerica_c230522.nc b/python/ctsm/test/testinputs/domain.lnd.fv0.9x1.25_gx1v7_f09_58x45_SouthAmerica_c230522.nc
new file mode 100644
index 0000000000..c3196bf7c4
--- /dev/null
+++ b/python/ctsm/test/testinputs/domain.lnd.fv0.9x1.25_gx1v7_f09_58x45_SouthAmerica_c230522.nc
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:177f512ee584cf4a33109f0163f36869a0ebce8c4b09212ae21fa4e3d5f2bfaa
+size 265248
diff --git a/python/ctsm/test/testinputs/domain.lnd.fv10x15_gx3v7.180321.nc b/python/ctsm/test/testinputs/domain.lnd.fv10x15_gx3v7.180321.nc
new file mode 100644
index 0000000000..38e10d5f32
--- /dev/null
+++ b/python/ctsm/test/testinputs/domain.lnd.fv10x15_gx3v7.180321.nc
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:14999fd73a1314bf98c50952ca554c47e95b91c06482d70af3211d0ab3c161c5
+size 47952
diff --git a/test/tools/README.testnames b/test/tools/README.testnames
index 983f883490..74dbe8e5f3 100644
--- a/test/tools/README.testnames
+++ b/test/tools/README.testnames
@@ -20,7 +20,7 @@ n is the configuration type:
 6 -- unused
 7 -- unused
 8 -- unused
-9 -- unused
+9 -- mesh_maker
 0 -- run_neon
 a -- modify_data
 b -- subset_data
@@ -44,9 +44,11 @@ m is the resolution
 a -- NEON YELL
 b -- NEON KONA
 c -- NEON OSBS
-d -- region1
-c -- single point from the 0.9x1.25 grid
-g -- unused
+d -- SouthAmerica
+e -- 1850PanTropics
+f -- PanBoreal
+g -- AlaskaTananaValley
+h -- single point from the 0.9x1.25 grid (Township SD)
 y -- 1.9x2.5 with transient 1850-2100 for rcp=2.6 and glacier-MEC on
 T -- 1x1_numaIA
 Z -- 10x15 with crop on
diff --git a/test/tools/input_tests_master b/test/tools/input_tests_master
index 3345aebbba..7da8c19803 100644
--- a/test/tools/input_tests_master
+++ b/test/tools/input_tests_master
@@ -14,10 +14,19 @@ smbb1 TSMscript_tools.sh site_and_regional subset_data subset_data_KONA
 blbb1 TBLscript_tools.sh site_and_regional subset_data subset_data_KONA
 smb81 TSMscript_tools.sh site_and_regional subset_data subset_data_US-UMB
 blb81 TBLscript_tools.sh site_and_regional subset_data subset_data_US-UMB
-smbc1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_US_pt
-blbc1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_US_pt
-smbd1 TSMscript_tools.sh site_and_regional subset_data subset_data_region1
-blbd1 TBLscript_tools.sh site_and_regional subset_data subset_data_region1
+smbh1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_1x1pt_townshipSD
+blbh1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_1x1pt_townshipSD
+smbd1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_58x45pt_SouthAmerica
+blbd1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_58x45pt_SouthAmerica
+smbe1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_90x288pt_1850PanTropics
+blbe1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_90x288pt_1850PanTropics
+smbf1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_37x288pt_PanBoreal
+blbf1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_37x288pt_PanBoreal
+smbg1 TSMscript_tools.sh site_and_regional subset_data subset_data_f09_4x9pt_AlaskaTananaValley
+blbg1 TBLscript_tools.sh site_and_regional subset_data subset_data_f09_4x9pt_AlaskaTananaValley
+
+sm901 TSMscript_tools.sh site_and_regional mesh_maker mesh_maker_fv09
+bl901 TBLscript_tools.sh site_and_regional mesh_maker mesh_maker_fv09
 
 smaa2 TSMscript_tools.sh site_and_regional modify_singlept_site_neon.py modify_data_YELL
 blaa2 TBLscript_tools.sh site_and_regional modify_singlept_site_neon.py modify_data_YELL
diff --git a/test/tools/nl_files/mesh_maker_fv09 b/test/tools/nl_files/mesh_maker_fv09
new file mode 100644
index 0000000000..7de951fee1
--- /dev/null
+++ b/test/tools/nl_files/mesh_maker_fv09
@@ -0,0 +1 @@
+ --input CSMDATA/atm/datm7/domain.lnd.fv0.9x1.25_gx1v6.090309.nc --lat yc --lon xc --overwrite --mask mask --area area --verbose
diff --git a/test/tools/nl_files/subset_data_KONA b/test/tools/nl_files/subset_data_KONA
index c3be007869..0df59b1b17 100644
--- a/test/tools/nl_files/subset_data_KONA
+++ b/test/tools/nl_files/subset_data_KONA
@@ -1 +1 @@
-point --lon 263.38956 --lat 39.1082 --site KONA --dompft 17 19 23 45 --pctpft 28 12 32 28 --crop --create-domain --create-surface --outdir EXEDIR/KONA_user-mod_and_data --user-mods-dir EXEDIR/KONA_user-mod_and_data --verbose  --inputdata-dir CSMDATA
+point --lon 263.38956 --lat 39.1082 --site KONA --dompft 17 19 23 45 --pctpft 28 12 32 28 --crop --create-surface --outdir EXEDIR/KONA_user-mod_and_data --user-mods-dir EXEDIR/KONA_user-mod_and_data --verbose  --inputdata-dir CSMDATA
diff --git a/test/tools/nl_files/subset_data_YELL b/test/tools/nl_files/subset_data_YELL
index 8295830c25..0d6960e7f5 100644
--- a/test/tools/nl_files/subset_data_YELL
+++ b/test/tools/nl_files/subset_data_YELL
@@ -1 +1 @@
-point --lon 250.45804 --lat 44.95597 --site YELL --dompft 1 --crop --create-domain --create-surface --outdir EXEDIR/YELL_user-mod_and_data --user-mods-dir EXEDIR/YELL_user-mod_and_data --verbose --inputdata-dir CSMDATA
+point --lon 250.45804 --lat 44.95597 --site YELL --dompft 1 --crop --uniform-snowpack --cap-saturation --create-surface --outdir EXEDIR/YELL_user-mod_and_data --user-mods-dir EXEDIR/YELL_user-mod_and_data --silent --inputdata-dir CSMDATA
diff --git a/test/tools/nl_files/subset_data_f09_1x1pt_townshipSD b/test/tools/nl_files/subset_data_f09_1x1pt_townshipSD
new file mode 100644
index 0000000000..aa25c07d1e
--- /dev/null
+++ b/test/tools/nl_files/subset_data_f09_1x1pt_townshipSD
@@ -0,0 +1 @@
+point --lon 257.5 --lat 43.822 --site f09_1x1pt_townshipSD  --include-nonveg --crop --create-datm --create-user-mods --datm-syr 2000 --datm-eyr 2000 --create-surface --outdir EXEDIR/f09_US_pt_user-mod_and_data --user-mods-dir EXEDIR/f09_US_pt_user-mod_and_data --verbose --inputdata-dir CSMDATA
diff --git a/test/tools/nl_files/subset_data_f09_37x288pt_PanBoreal b/test/tools/nl_files/subset_data_f09_37x288pt_PanBoreal
new file mode 100644
index 0000000000..448b5052d6
--- /dev/null
+++ b/test/tools/nl_files/subset_data_f09_37x288pt_PanBoreal
@@ -0,0 +1 @@
+region --lat1 55 --lat2 89.1 --lon1 0 --lon2 360 --create-mesh --create-surface --create-domain --create-user-mods --verbose --overwrite --reg f09_37x288pt_PanBoreal --inputdata-dir CSMDATA
diff --git a/test/tools/nl_files/subset_data_f09_4x9pt_AlaskaTananaValley b/test/tools/nl_files/subset_data_f09_4x9pt_AlaskaTananaValley
new file mode 100644
index 0000000000..9928d78429
--- /dev/null
+++ b/test/tools/nl_files/subset_data_f09_4x9pt_AlaskaTananaValley
@@ -0,0 +1 @@
+region --lat1 62 --lat2 66 --lon1 -152 --lon2 -141 --create-mesh --create-domain --create-surface --create-user-mods --verbose --overwrite --reg f09_4x9pt_AlaskaTananaValley --inputdata-dir CSMDATA
diff --git a/test/tools/nl_files/subset_data_f09_58x45pt_SouthAmerica b/test/tools/nl_files/subset_data_f09_58x45pt_SouthAmerica
new file mode 100644
index 0000000000..201dd2c76c
--- /dev/null
+++ b/test/tools/nl_files/subset_data_f09_58x45pt_SouthAmerica
@@ -0,0 +1 @@
+region --lat1 -40 --lat2 15 --lon1 275 --lon2 330 --create-mesh --create-surface --create-user-mods --create-domain --create-landuse --verbose --overwrite --reg f09_58x45_SouthAmerica --inputdata-dir CSMDATA
diff --git a/test/tools/nl_files/subset_data_f09_90x288pt_1850PanTropics b/test/tools/nl_files/subset_data_f09_90x288pt_1850PanTropics
new file mode 100644
index 0000000000..1c9d5eace9
--- /dev/null
+++ b/test/tools/nl_files/subset_data_f09_90x288pt_1850PanTropics
@@ -0,0 +1 @@
+region --lat1 -55 --lat2 30 --lon1 0 --lon2 360 --crop --create-surface --create-domain --create-mesh --overwrite --reg f09_90x288pt_1850PanTropics --inputdata-dir CSMDATA --cfg-file CTSM_ROOT/tools/mksurfdata_map/default_data_1850.cfg --verbose
diff --git a/test/tools/nl_files/subset_data_f09_US_pt b/test/tools/nl_files/subset_data_f09_US_pt
deleted file mode 100644
index bf6d5e2861..0000000000
--- a/test/tools/nl_files/subset_data_f09_US_pt
+++ /dev/null
@@ -1 +0,0 @@
-point --lon 257.5 --lat 43.822 --site 1x1_  --include-nonveg --crop --create-landuse --create-datm --create-user-mods --datm-syr 2000 --datm-eyr 2000 --create-surface --outdir EXEDIR/f09_US_pt_user-mod_and_data --user-mods-dir EXEDIR/f09_US_pt_user-mod_and_data --verbose --inputdata-dir CSMDATA
diff --git a/test/tools/nl_files/subset_data_region1 b/test/tools/nl_files/subset_data_region1
deleted file mode 100644
index fce83f0e2e..0000000000
--- a/test/tools/nl_files/subset_data_region1
+++ /dev/null
@@ -1 +0,0 @@
-region --lat1 -40 --lat2 15 --lon1 275 --lon2 330 --create-domain --create-surface --create-landuse --verbose --overwrite --reg test1 --inputdata-dir CSMDATA
diff --git a/test/tools/test_driver.sh b/test/tools/test_driver.sh
index 6b1811c8ab..f93301a530 100755
--- a/test/tools/test_driver.sh
+++ b/test/tools/test_driver.sh
@@ -431,7 +431,7 @@ else
     fi
 fi
 
-# Setup conda environement
+# Setup conda environment
 conda activate ctsm_pylib
 if [ \$? -ne 0 ]; then
    echo "ERROR: Trouble activating the ctsm_pylib conda environment, be sure it's setup with \$CLM_ROOT/py_env_create, then rerun"
diff --git a/tools/site_and_regional/README b/tools/site_and_regional/README
index 7b36b7d2b9..1fa730c96c 100644
--- a/tools/site_and_regional/README
+++ b/tools/site_and_regional/README
@@ -20,6 +20,12 @@ subset_data
     For extracting domain files, surface dataset, and DATM files at a region, use:
         ./subset_data region
 
+mesh_maker
+    This script creates a mesh file from a netcdf file with valid lats and lons.
+
+mesh_plotter
+    This script plots a mesh file
+
 modify_singlept_site_neon.py
     After running subset_data.py overwrite some fields with site-specific
     data for neon sites.
diff --git a/tools/site_and_regional/default_data_2000.cfg b/tools/site_and_regional/default_data_2000.cfg
index 0ba7f8758b..d736f2e933 100644
--- a/tools/site_and_regional/default_data_2000.cfg
+++ b/tools/site_and_regional/default_data_2000.cfg
@@ -18,6 +18,8 @@ tpqwname = CLMGSWP3v1.TPQW
 dir = lnd/clm2/surfdata_esmf/ctsm5.2.0
 surfdat_16pft = surfdata_0.9x1.25_hist_2000_16pfts_c240216.nc
 surfdat_78pft = surfdata_0.9x1.25_hist_2000_78pfts_c240216.nc
+mesh_dir = share/meshes/
+mesh_surf = fv0.9x1.25_141008_ESMFmesh.nc
 
 [landuse]
 dir = lnd/clm2/surfdata_esmf/ctsm5.2.0
diff --git a/tools/site_and_regional/mesh_maker b/tools/site_and_regional/mesh_maker
new file mode 100755
index 0000000000..24b36de995
--- /dev/null
+++ b/tools/site_and_regional/mesh_maker
@@ -0,0 +1,40 @@
+#!/usr/bin/env python3
+"""
+This is a just top-level skeleton script that calls
+mesh_maker.py.
+The original code (mesh_maker.py) is located under the
+python/ctsm folder.
+
+For full instructions on how to run the code and different options,
+please check the python/ctsm/mesh_maker.py file.
+
+----------------------------------------------------------------
+Instructions for running using conda python environments:
+
+../../py_env_create
+conda activate ctsm_py
+|------------------------------------------------------------------|
+|---------------------  Instructions  -----------------------------|
+|------------------------------------------------------------------|
+This script creates ESMF unstructured GRID (mesh file) from a netcdf
+file with valid lats and lons. Provided lats and lons can be 1D or 2D.
+
+For example for running WRF-CTSM cases, the user can create a mesh
+file for their domain :
+    ./mesh_maker --input wrfinput_d01 --output my_region
+        --lat XLAT --lon XLONG --verbose
+
+"""
+import os
+import sys
+
+# -- add python/ctsm  to path
+_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.mesh_maker import main
+
+if __name__ == "__main__":
+    main()
diff --git a/tools/site_and_regional/mesh_plotter b/tools/site_and_regional/mesh_plotter
new file mode 100755
index 0000000000..2f195a2b31
--- /dev/null
+++ b/tools/site_and_regional/mesh_plotter
@@ -0,0 +1,33 @@
+#!/usr/bin/env python3
+"""
+This is a just top-level skeleton script that calls
+mesh_plotter.py.
+The original code (mesh_plotter.py) is located under the
+python/ctsm folder.
+
+For full instructions on how to run the code and different options,
+please check the python/ctsm/mesh_plotter.py file.
+
+----------------------------------------------------------------
+Instructions for running using conda python environments:
+
+../../py_env_create --dask
+conda activate ctsm_py_wdask
+|------------------------------------------------------------------|
+|---------------------  Instructions  -----------------------------|
+|------------------------------------------------------------------|
+This script plots a  ESMF unstructured GRID (mesh file).
+"""
+import os
+import sys
+
+# -- add python/ctsm  to path
+_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.mesh_plotter import main
+
+if __name__ == "__main__":
+    main()