diff --git a/Makefile b/Makefile
index 47a7c0b..735f5ac 100644
--- a/Makefile
+++ b/Makefile
@@ -16,11 +16,11 @@ gnu:   # BUILDTARGET GNU Fortran, C, and C++ compilers
 	"CC_SERIAL = gcc" \
 	"CXX_SERIAL = g++" \
 	"FFLAGS_PROMOTION = -fdefault-real-8 -fdefault-double-8" \
-	"FFLAGS_OPT = -std=f2008 -O3 -ffree-line-length-none -fconvert=big-endian -ffree-form" \
+	"FFLAGS_OPT = -std=f2008 -O3 -fallow-argument-mismatch -ffree-line-length-none -fconvert=big-endian -ffree-form" \
 	"CFLAGS_OPT = -O3" \
 	"CXXFLAGS_OPT = -O3" \
 	"LDFLAGS_OPT = -O3" \
-	"FFLAGS_DEBUG = -std=f2008 -g -ffree-line-length-none -fconvert=big-endian -ffree-form -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow" \
+	"FFLAGS_DEBUG = -g -ffree-line-length-none -fallow-argument-mismatch -fconvert=big-endian -ffree-form -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow" \
 	"CFLAGS_DEBUG = -g" \
 	"CXXFLAGS_DEBUG = -g" \
 	"LDFLAGS_DEBUG = -g" \
@@ -154,7 +154,7 @@ nvhpc:   # BUILDTARGET NVIDIA HPC SDK
 	"FFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -byteswapio -Mfree -Ktrap=divz,fp,inv,ovf -traceback" \
 	"CFLAGS_DEBUG = -O0 -g -traceback" \
 	"CXXFLAGS_DEBUG = -O0 -g -traceback" \
-	"LDFLAGS_DEBUG = -O0 -g -Mbounds -Ktrap=divz,fp,inv,ovf -traceback" \
+	"LDFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -Ktrap=divz,fp,inv,ovf -traceback" \
 	"FFLAGS_OMP = -mp" \
 	"CFLAGS_OMP = -mp" \
 	"FFLAGS_ACC = -Mnofma -acc -gpu=cc70,cc80 -Minfo=accel" \
@@ -184,7 +184,7 @@ pgi:   # BUILDTARGET PGI compiler suite
 	"FFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -byteswapio -Mfree -Ktrap=divz,fp,inv,ovf -traceback" \
 	"CFLAGS_DEBUG = -O0 -g -traceback" \
 	"CXXFLAGS_DEBUG = -O0 -g -traceback" \
-	"LDFLAGS_DEBUG = -O0 -g -Mbounds -Ktrap=divz,fp,inv,ovf -traceback" \
+	"LDFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -Ktrap=divz,fp,inv,ovf -traceback" \
 	"FFLAGS_OMP = -mp" \
 	"CFLAGS_OMP = -mp" \
 	"FFLAGS_ACC = -Mnofma -acc -Minfo=accel" \
@@ -216,7 +216,7 @@ pgi-summit:   # BUILDTARGET PGI compiler suite w/OpenACC options for ORNL Summit
 	"FFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -byteswapio -Mfree -Ktrap=divz,fp,inv,ovf -traceback" \
 	"CFLAGS_DEBUG = -O0 -g -traceback" \
 	"CXXFLAGS_DEBUG = -O0 -g -traceback" \
-	"LDFLAGS_DEBUG = -O0 -g -Mbounds -Ktrap=divz,fp,inv,ovf -traceback" \
+	"LDFLAGS_DEBUG = -O0 -g -Mbounds -Mchkptr -Ktrap=divz,fp,inv,ovf -traceback" \
 	"FFLAGS_OMP = -mp" \
 	"CFLAGS_OMP = -mp" \
 	"PICFLAG = -fpic" \
@@ -667,10 +667,10 @@ intel:   # BUILDTARGET Intel oneAPI Fortran, C, and C++ compiler suite
 	"CFLAGS_OPT = -O3" \
 	"CXXFLAGS_OPT = -O3" \
 	"LDFLAGS_OPT = -O3" \
-	"FFLAGS_DEBUG = -g -convert big_endian -free -check bounds,pointers,arg_temp_created,format,shape,contiguous -fpe0 -traceback" \
+	"FFLAGS_DEBUG = -g -convert big_endian -free -check all -fpe0 -traceback" \
 	"CFLAGS_DEBUG = -g -traceback" \
 	"CXXFLAGS_DEBUG = -g -traceback" \
-	"LDFLAGS_DEBUG = -g -traceback" \
+	"LDFLAGS_DEBUG = -g -fpe0 -traceback" \
 	"FFLAGS_OMP = -qopenmp" \
 	"CFLAGS_OMP = -qopenmp" \
 	"PICFLAG = -fpic" \
@@ -681,6 +681,86 @@ intel:   # BUILDTARGET Intel oneAPI Fortran, C, and C++ compiler suite
 	"OPENMP = $(OPENMP)" \
 	"CPPFLAGS = $(MODEL_FORMULATION) -D_MPI" )
 
+intel-xd2000:
+	( $(MAKE) all \
+        "FC_PARALLEL = ftn" \
+        "CC_PARALLEL = cc" \
+        "CXX_PARALLEL = CC" \
+        "FC_SERIAL = ftn" \
+        "CC_SERIAL = cc" \
+        "CXX_SERIAL = CC" \
+        "FFLAGS_PROMOTION = -real-size 64" \
+        "FFLAGS_OPT = -O3 -traceback -convert big_endian -FR -march=core-avx2 -mtune=core-avx2" \
+        "CFLAGS_OPT = -O3 -traceback -std=gnu90" \
+        "CXXFLAGS_OPT = -O3 -traceback" \
+        "LDFLAGS_OPT = -O3 -traceback" \
+        "FFLAGS_DEBUG = -g -convert big_endian -FR -CU -CB -check all -fpe0 -traceback" \
+        "CFLAGS_DEBUG = -g -traceback" \
+        "CXXFLAGS_DEBUG = -g -traceback" \
+        "LDFLAGS_DEBUG = -g -fpe0 -traceback" \
+        "FFLAGS_OMP = -qopenmp" \
+        "CFLAGS_OMP = -qopenmp" \
+        "CORE = $(CORE)" \
+        "DEBUG = $(DEBUG)" \
+        "USE_PAPI = $(USE_PAPI)" \
+        "OPENMP = $(OPENMP)" \
+        "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE" )
+
+intel2-xd2000:
+	( $(MAKE) all \
+        "FC_PARALLEL = ftn" \
+        "CC_PARALLEL = cc" \
+        "CXX_PARALLEL = CC" \
+        "FC_SERIAL = ftn" \
+        "CC_SERIAL = cc" \
+        "CXX_SERIAL = CC" \
+        "FFLAGS_PROMOTION = -real-size 64" \
+        "FFLAGS_OPT = -O2 -traceback -convert big_endian -FR -march=core-avx2 -mtune=core-avx2" \
+        "CFLAGS_OPT = -O2 -traceback -std=gnu90" \
+        "CXXFLAGS_OPT = -O2 -traceback" \
+        "LDFLAGS_OPT = -O2 -traceback" \
+        "FFLAGS_DEBUG = -g -convert big_endian -FR -CU -CB -check all -fpe0 -traceback" \
+        "CFLAGS_DEBUG = -g -traceback" \
+        "CXXFLAGS_DEBUG = -g -traceback" \
+        "LDFLAGS_DEBUG = -g -fpe0 -traceback" \
+        "FFLAGS_OMP = -qopenmp" \
+        "CFLAGS_OMP = -qopenmp" \
+        "CORE = $(CORE)" \
+        "DEBUG = $(DEBUG)" \
+        "USE_PAPI = $(USE_PAPI)" \
+        "OPENMP = $(OPENMP)" \
+        "CPPFLAGS = $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE" )
+
+gfortran-xd2000:   # BUILDTARGET GNU Fortran, C, and C++ compilers
+	( $(MAKE) all \
+	"FC_PARALLEL = ftn" \
+	"CC_PARALLEL = cc" \
+	"CXX_PARALLEL = CC" \
+	"FC_SERIAL = ftn" \
+	"CC_SERIAL = cc" \
+	"CXX_SERIAL = CC" \
+	"FFLAGS_PROMOTION = -fdefault-real-8 -fdefault-double-8" \
+	"FFLAGS_OPT = -O3 -ffree-line-length-none -fconvert=big-endian -ffree-form -fallow-argument-mismatch" \
+	"CFLAGS_OPT = -O3" \
+	"CXXFLAGS_OPT = -O3" \
+	"LDFLAGS_OPT = -O3" \
+	"FFLAGS_DEBUG = -g -ffree-line-length-none -fconvert=big-endian -ffree-form -fallow-argument-mismatch -fcheck=all -fbacktrace -ffpe-trap=invalid,zero,overflow" \
+	"CFLAGS_DEBUG = -g" \
+	"CXXFLAGS_DEBUG = -g" \
+	"LDFLAGS_DEBUG = -g" \
+	"FFLAGS_OMP = -fopenmp" \
+	"CFLAGS_OMP = -fopenmp" \
+	"FFLAGS_ACC =" \
+	"CFLAGS_ACC =" \
+	"PICFLAG = -fPIC" \
+	"BUILD_TARGET = $(@)" \
+	"CORE = $(CORE)" \
+	"DEBUG = $(DEBUG)" \
+	"USE_PAPI = $(USE_PAPI)" \
+	"OPENMP = $(OPENMP)" \
+	"OPENACC = $(OPENACC)" \
+	"CPPFLAGS = $(MODEL_FORMULATION) -D_MPI" )
+
 CPPINCLUDES =
 FCINCLUDES =
 LIBS =
@@ -690,6 +770,9 @@ ifneq "$(PIO)" ""
 # Regardless of PIO library version, look for a lib subdirectory of PIO path
 # NB: PIO_LIB is used later, so we don't just set LIBS directly
 #
+#ifneq ($(wildcard $(PIO)/lib64), )
+#	PIO_LIB = $(PIO)/lib64
+#else
 ifneq ($(wildcard $(PIO)/lib), )
 	PIO_LIB = $(PIO)/lib
 else
@@ -737,10 +820,10 @@ else # Not using PIO, using SMIOL
 endif
 
 ifneq "$(NETCDF)" ""
-ifneq ($(wildcard $(NETCDF)/lib/libnetcdf.*), )
+ifneq ($(wildcard $(NETCDF)/lib), )
 	NETCDFLIBLOC = lib
 endif
-ifneq ($(wildcard $(NETCDF)/lib64/libnetcdf.*), )
+ifneq ($(wildcard $(NETCDF)/lib64), )
 	NETCDFLIBLOC = lib64
 endif
 	CPPINCLUDES += -I$(NETCDF)/include
@@ -761,10 +844,10 @@ endif
 
 
 ifneq "$(PNETCDF)" ""
-ifneq ($(wildcard $(PNETCDF)/lib/libpnetcdf.*), )
+ifneq ($(wildcard $(PNETCDF)/lib), )
 	PNETCDFLIBLOC = lib
 endif
-ifneq ($(wildcard $(PNETCDF)/lib64/libpnetcdf.*), )
+ifneq ($(wildcard $(PNETCDF)/lib64), )
 	PNETCDFLIBLOC = lib64
 endif
 	CPPINCLUDES += -I$(PNETCDF)/include
@@ -958,6 +1041,28 @@ else
 	OPENACC_MESSAGE="MPAS was built without OpenACC accelerator support."
 endif
 
+ifneq ($(wildcard .mpas_core_*), ) # CHECK FOR BUILT CORE
+
+ifneq ($(wildcard .mpas_core_$(CORE)), ) # CHECK FOR SAME CORE AS ATTEMPTED BUILD.
+	override AUTOCLEAN=false
+	CONTINUE=true
+else
+	LAST_CORE=`cat .mpas_core_*`
+
+ifeq "$(AUTOCLEAN)" "true" # CHECK FOR CLEAN PRIOR TO BUILD OF A NEW CORE.
+	CONTINUE=true
+	AUTOCLEAN_MESSAGE="Infrastructure was cleaned prior to building ."
+else
+	CONTINUE=false
+endif # END OF AUTOCLEAN CHECK
+
+endif # END OF CORE=LAST_CORE CHECK
+
+else
+
+	override AUTOCLEAN=false
+	CONTINUE=true
+endif # END IF BUILT CORE CHECK
 
 ifneq ($(wildcard namelist.$(NAMELIST_SUFFIX)), ) # Check for generated namelist file.
 	NAMELIST_MESSAGE="A default namelist file (namelist.$(NAMELIST_SUFFIX).defaults) has been generated, but namelist.$(NAMELIST_SUFFIX) has not been modified."
@@ -1014,119 +1119,12 @@ report_builds:
 	@echo "CORE=$(CORE)"
 endif
 
+ifeq "$(CONTINUE)" "true"
 all: mpas_main
-
+else
+all: clean_core
 endif
 
-#
-# The rebuild_check target determines whether the shared framework or $(CORE) were
-# previously compiled with incompatible options, and stops the build with an error
-# message if so.
-#
-rebuild_check:
-	@#
-	@# Write current build options to a file .build_opts.tmp, to later be
-	@# compared with build options use for the shared framework or core.
-	@# Only build options that affect compatibility are written, while options
-	@# like $(RM), $(BUILD_TARGET), and $(CORE) are not.
-	@#
-	$(shell printf "FC=$(FC)\n$\
-	CC=$(CC)\n$\
-	CXX=$(CXX)\n$\
-	SFC=$(SFC)\n$\
-	SCC=$(SCC)\n$\
-	CFLAGS=$(CFLAGS)\n$\
-	CXXFLAGS=$(CXXFLAGS)\n$\
-	FFLAGS=$(FFLAGS)\n$\
-	LDFLAGS=$(LDFLAGS)\n$\
-	CPPFLAGS=$(CPPFLAGS)\n$\
-	LIBS=$(LIBS)\n$\
-	CPPINCLUDES=$(CPPINCLUDES)\n$\
-	OPENMP=$(OPENMP)\n$\
-	OPENMP_OFFLOAD=$(OPENMP_OFFLOAD)\n$\
-	OPENACC=$(OPENACC)\n$\
-	TAU=$(TAU)\n$\
-	PICFLAG=$(PICFLAG)\n$\
-	TIMER_LIB=$(TIMER_LIB)\n$\
-	GEN_F90=$(GEN_F90)\n" | sed 's/-DMPAS_EXE_NAME=[^[:space:]]*//' | sed 's/-DMPAS_NAMELIST_SUFFIX=[^[:space:]]*//' | sed 's/-DCORE_[^[:space:]]*//' | sed 's/-DMPAS_GIT_VERSION=[^[:space:]]*//' > .build_opts.tmp )
-
-	@#
-	@# PREV_BUILD is set to "OK" if the shared framework and core are either
-	@# clean or were previously compiled with compatible options. Otherwise,
-	@# PREV_BUILD is set to "shared framework" if the shared framework was
-	@# built with incompatible options, or "$(CORE) core" if the core was
-	@# built with incompatible options.
-	@#
-	$(eval PREV_BUILD := $(shell $\
-		if [ -f ".build_opts.framework" ]; then $\
-			cmp -s .build_opts.tmp .build_opts.framework; $\
-			if [ $$? -eq 0 ]; then $\
-				stat=0; $\
-			else $\
-				stat=1; $\
-				x="shared framework"; $\
-				if [ "$(AUTOCLEAN)" = "true" ]; then $\
-					cp .build_opts.tmp .build_opts.framework; $\
-				fi; $\
-			fi $\
-		else $\
-			stat=0; $\
-			cp .build_opts.tmp .build_opts.framework; $\
-		fi; $\
-                : ; $\
-                : At this this point, stat is already set, and we should only ; $\
-                : set it to 1 but never to 0, as that might mask an incompatibility ; $\
-                : in the framework build. ; $\
-                : ; $\
-		if [ -f ".build_opts.$(CORE)" ]; then $\
-			cmp -s .build_opts.tmp .build_opts.$(CORE); $\
-			if [ $$? -ne 0 ]; then $\
-				stat=1; $\
-				if [ "$$x" = "" ]; then $\
-					x="$(CORE) core"; $\
-				else $\
-					x="$$x and $(CORE) core"; $\
-				fi; $\
-				if [ "$(AUTOCLEAN)" = "true" ]; then $\
-					cp .build_opts.tmp .build_opts.$(CORE); $\
-				fi; $\
-			fi; $\
-		else $\
-			if [ $$stat -eq 0 ]; then $\
-				cp .build_opts.tmp .build_opts.$(CORE); $\
-			fi; $\
-		fi; $\
-		rm -f .build_opts.tmp; $\
-		if [ $$stat -eq 1 ]; then $\
-			printf "$$x"; $\
-		else $\
-			printf "OK"; $\
-		fi; $\
-	))
-
-	$(if $(findstring and,$(PREV_BUILD)),$(eval VERB=were),$(eval VERB=was))
-ifeq "$(AUTOCLEAN)" "true"
-	$(if $(findstring framework,$(PREV_BUILD)),$(eval AUTOCLEAN_DEPS+=clean_shared))
-	$(if $(findstring core,$(PREV_BUILD)),$(eval AUTOCLEAN_DEPS+=clean_core))
-	$(if $(findstring OK,$(PREV_BUILD)), $(eval override AUTOCLEAN=false), )
-	$(eval AUTOCLEAN_MESSAGE=The $(PREV_BUILD) $(VERB) cleaned and re-compiled.)
-else
-	$(if $(findstring OK,$(PREV_BUILD)), \
-	, \
-	$(info ************************************************************************) \
-	$(info The $(PREV_BUILD) $(VERB) previously compiled with ) \
-	$(info incompatible options. Please do one of the following:) \
-	$(info ) \
-	$(info   - Clean the $(CORE) core, which will also cause the shared) \
-	$(info     framework to be cleaned; then compile the $(CORE) core.) \
-	$(info ) \
-	$(info   or)\
-	$(info ) \
-	$(info   - Add AUTOCLEAN=true to the build command to automatically clean) \
-	$(info     and re-compile the $(PREV_BUILD).) \
-	$(info ) \
-	$(info ************************************************************************) \
-	$(error ))
 endif
 
 
@@ -1266,7 +1264,7 @@ ifeq "$(OPENACC)" "true"
 endif # OPENACC eq true
 
 
-pio_test: openmp_test openacc_test pnetcdf_test
+pio_test: openmp_test openacc_test
 	@#
 	@# PIO_VERS will be set to:
 	@#  0 if no working PIO library was detected (and .piotest.log will contain error messages)
@@ -1352,10 +1350,9 @@ mpi_f08_test:
 	$(info Checking for mpi_f08 support...)
 	$(eval MPAS_MPI_F08 := $(shell $\
 		printf "program main\n$\
-		        &   use mpi_f08, only : MPI_Init, MPI_Comm, MPI_INTEGER, MPI_Datatype\n$\
+		        &   use mpi_f08, only : MPI_Init, MPI_Comm\n$\
 		        &   integer :: ierr\n$\
 		        &   type (MPI_Comm) :: comm\n$\
-		        &   type (MPI_Datatype), parameter :: MPI_INTEGERKIND = MPI_INTEGER\n$\
 		        &   call MPI_Init(ierr)\n$\
 		        end program main\n" | sed 's/&/ /' > mpi_f08.f90; $\
 		$\
@@ -1374,57 +1371,20 @@ mpi_f08_test:
 	$(if $(findstring 1,$(MPAS_MPI_F08)), $(eval MPI_F08_MESSAGE = "Using the mpi_f08 module."), )
 	$(if $(findstring 1,$(MPAS_MPI_F08)), $(info mpi_f08 module detected.))
 
-
-pnetcdf_test:
-	@#
-	@# Create test C programs that look for PNetCDF header file and some symbols in it
-	@#
-ifneq "$(PNETCDF)" ""
-	@echo "Checking for a working PnetCDF library..."
-	@printf "#include \"pnetcdf.h\"\n\
-			&#include \"mpi.h\"\n\
-			&int main(){\n\
-			&    int err, ncid;\n\
-			&    err = ncmpi_create(MPI_COMM_WORLD, \"foo.nc\", NC_NOCLOBBER, MPI_INFO_NULL, &ncid);\n\
-			&    return 0;\n\
-			&}\n" | sed 's/&/ /' > pnetcdf.c
-	@( $(CC) pnetcdf.c $(CPPINCLUDES) $(CFLAGS) $(LDFLAGS) -L$(PNETCDF)/$(PNETCDFLIBLOC) -lpnetcdf  -o pnetcdf.out > pnetcdf.log 2>&1; \
-	   if [ $$? -eq 0 ] ; then \
-	       echo "$(CC) can compile test PnetCDF C program."; \
-	   else \
-	       echo "*********************************************************"; \
-	       echo "ERROR: Test PnetCDF C program could not be compiled by $(CC)."; \
-	       echo "Please ensure you have a working PnetCDF library installed."; \
-	       echo ""; \
-	       echo "The following compilation command failed with errors:" ; \
-	       echo "$(CC) pnetcdf.c $(CPPINCLUDES) $(CFLAGS) $(LDFLAGS) -L$(PNETCDF)/$(PNETCDFLIBLOC) -lpnetcdf -o pnetcdf.out"; \
-	       echo ""; \
-	       echo "Test program pnetcdf.c and output pnetcdf.log have been left"; \
-	       echo "in the top-level MPAS directory for further debugging"; \
-	       echo "*********************************************************"; \
-	       rm -f pnetcdf.out; exit 1; \
-	   fi )
-
-	@rm -f pnetcdf.c pnetcdf.out pnetcdf.log
-else
-	@echo "*********************************************************"; \
-	 echo "ERROR: The PNETCDF environment variable isn't set."; \
-	 echo "Please set this variable to where PnetCDF is installed."; \
-	 echo "*********************************************************"; \
-	 exit 1
-endif
-
-
 ifneq "$(PIO)" ""
-MAIN_DEPS = rebuild_check openmp_test openacc_test pnetcdf_test pio_test mpi_f08_test
+MAIN_DEPS = openmp_test openacc_test pio_test mpi_f08_test
 override CPPFLAGS += "-DMPAS_PIO_SUPPORT"
 else
-MAIN_DEPS = rebuild_check openmp_test openacc_test pnetcdf_test mpi_f08_test
+MAIN_DEPS = openmp_test openacc_test mpi_f08_test
 IO_MESSAGE = "Using the SMIOL library."
 override CPPFLAGS += "-DMPAS_SMIOL_SUPPORT"
 endif
 
+
 mpas_main: $(MAIN_DEPS)
+ifeq "$(AUTOCLEAN)" "true"
+	$(RM) .mpas_core_*
+endif
 	cd src; $(MAKE) FC="$(FC)" \
                  CC="$(CC)" \
                  CXX="$(CXX)" \
@@ -1443,11 +1403,11 @@ mpas_main: $(MAIN_DEPS)
                  FCINCLUDES="$(FCINCLUDES)" \
                  CORE="$(CORE)"\
                  AUTOCLEAN="$(AUTOCLEAN)" \
-                 AUTOCLEAN_DEPS="$(AUTOCLEAN_DEPS)" \
                  GEN_F90="$(GEN_F90)" \
                  NAMELIST_SUFFIX="$(NAMELIST_SUFFIX)" \
                  EXE_NAME="$(EXE_NAME)"
 
+	@echo "$(EXE_NAME)" > .mpas_core_$(CORE)
 	if [ -e src/$(EXE_NAME) ]; then mv src/$(EXE_NAME) .; fi
 	( cd src/core_$(CORE); $(MAKE) ROOT_DIR="$(PWD)" post_build )
 	@echo "*******************************************************************************"
@@ -1469,13 +1429,11 @@ endif
 	@echo $(IO_MESSAGE)
 	@echo "*******************************************************************************"
 clean:
-	cd src; $(MAKE) clean RM="$(RM)" CORE="$(CORE)" AUTOCLEAN="$(AUTOCLEAN)"
+	cd src; $(MAKE) clean RM="$(RM)" CORE="$(CORE)"
+	$(RM) .mpas_core_*
 	$(RM) $(EXE_NAME)
 	$(RM) namelist.$(NAMELIST_SUFFIX).defaults
 	$(RM) streams.$(NAMELIST_SUFFIX).defaults
-	if [ -f .build_opts.framework ]; then $(RM) .build_opts.framework; fi
-	if [ -f .build_opts.$(CORE) ]; then $(RM) .build_opts.$(CORE); fi
-
 core_error:
 	@echo ""
 	@echo "*******************************************************************************"
@@ -1486,6 +1444,26 @@ core_error:
 	exit 1
 error: errmsg
 
+clean_core:
+	@echo ""
+	@echo "*******************************************************************************"
+	@echo " The MPAS infrastructure is currently built for the $(LAST_CORE) core."
+	@echo " Before building the $(CORE) core, please do one of the following."
+	@echo ""
+	@echo ""
+	@echo " To remove the $(LAST_CORE)_model executable and clean the MPAS infrastructure, run:"
+	@echo "      make clean CORE=$(LAST_CORE)"
+	@echo ""
+	@echo " To preserve all executables except $(CORE)_model and clean the MPAS infrastructure, run:"
+	@echo "      make clean CORE=$(CORE)"
+	@echo ""
+	@echo " Alternatively, AUTOCLEAN=true can be appended to the make command to force a clean,"
+	@echo " build a new $(CORE)_model executable, and preserve all other executables."
+	@echo ""
+	@echo "*******************************************************************************"
+	@echo ""
+	exit 1
+
 else # CORE IF
 
 all: error
@@ -1513,7 +1491,7 @@ errmsg:
 	@echo "    DEBUG=true    - builds debug version. Default is optimized version."
 	@echo "    USE_PAPI=true - builds version using PAPI for timers. Default is off."
 	@echo "    TAU=true      - builds version using TAU hooks for profiling. Default is off."
-	@echo "    AUTOCLEAN=true - Enables automatic cleaning and re-compilation of code as needed."
+	@echo "    AUTOCLEAN=true    - forces a clean of infrastructure prior to build new core."
 	@echo "    GEN_F90=true  - Generates intermediate .f90 files through CPP, and builds with them."
 	@echo "    TIMER_LIB=opt - Selects the timer library interface to be used for profiling the model. Options are:"
 	@echo "                    TIMER_LIB=native - Uses native built-in timers in MPAS"
@@ -1530,4 +1508,3 @@ errmsg:
 ifdef CORE
 	exit 1
 endif
-
diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 1dae832..048fc58 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -508,6 +508,7 @@
 #endif
 #ifdef DO_PHYSICS
 			<var name="isltyp"/>
+			<var name="isctyp"/>
 			<var name="ivgtyp"/>
 			<var name="mminlu"/>
 			<var name="isice_lu"/>
@@ -4462,7 +4463,10 @@
                 <!-- ================================================================================================== -->
 
                 <var name="isltyp" type="integer" dimensions="nCells" units="unitless"
-                     description="dominant soil category"/>
+                     description="dominant soil texture category"/>
+
+                <var name="isctyp" type="integer" dimensions="nCells" units="unitless"
+                     description="dominant soil colour category"/>
 
                 <var name="ivgtyp" type="integer" dimensions="nCells" units="unitless"
                      description="dominant vegetation category"/>
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F
index 0116dcf..ba5faf2 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F
@@ -139,6 +139,7 @@ subroutine allocate_lsm
  if(.not.allocated(cpm_p)       ) allocate(cpm_p(ims:ime,jms:jme)       )
  if(.not.allocated(cqs2_p)      ) allocate(cqs2_p(ims:ime,jms:jme)      )
  if(.not.allocated(isltyp_p)    ) allocate(isltyp_p(ims:ime,jms:jme)    )
+ if(.not.allocated(isctyp_p)    ) allocate(isctyp_p(ims:ime,jms:jme)    )
  if(.not.allocated(ivgtyp_p)    ) allocate(ivgtyp_p(ims:ime,jms:jme)    )
  if(.not.allocated(glw_p)       ) allocate(glw_p(ims:ime,jms:jme)       )
  if(.not.allocated(grdflx_p)    ) allocate(grdflx_p(ims:ime,jms:jme)    )
@@ -216,6 +217,7 @@ subroutine deallocate_lsm
  if(allocated(gsw_p)       ) deallocate(gsw_p       )
  if(allocated(hfx_p)       ) deallocate(hfx_p       )
  if(allocated(isltyp_p)    ) deallocate(isltyp_p    )
+ if(allocated(isctyp_p)    ) deallocate(isctyp_p    )
  if(allocated(ivgtyp_p)    ) deallocate(ivgtyp_p    )
  if(allocated(lai_p)       ) deallocate(lai_p       )
  if(allocated(lh_p)        ) deallocate(lh_p        )
@@ -279,7 +281,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
  character(len=StrKIND),pointer:: config_microp_scheme,    &
                                   config_convection_scheme
 
- integer,dimension(:),pointer:: isltyp,ivgtyp
+ integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp
 
  real(kind=RKIND),dimension(:),pointer  :: acsnom,acsnow,canwat,chs,chs2,chklowq,cpm,cqs2,glw,   &
                                            grdflx,gsw,hfx,lai,lh,noahres,potevp,qfx,qgh,qsfc,    &
@@ -336,6 +338,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
  call mpas_pool_get_array(diag_physics,'znt'       ,znt       )
 
  call mpas_pool_get_array(sfc_input,'isltyp'    ,isltyp    )
+ call mpas_pool_get_array(sfc_input,'isctyp'    ,isctyp    )
  call mpas_pool_get_array(sfc_input,'ivgtyp'    ,ivgtyp    )
  call mpas_pool_get_array(sfc_input,'shdmin'    ,shdmin    )
  call mpas_pool_get_array(sfc_input,'shdmax'    ,shdmax    )
@@ -409,6 +412,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
     znt_p(i,j)        = znt(i)
 
     isltyp_p(i,j)     = isltyp(i)
+    isctyp_p(i,j)     = isctyp(i)
     ivgtyp_p(i,j)     = ivgtyp(i)
     shdmin_p(i,j)     = shdmin(i)
     shdmax_p(i,j)     = shdmax(i)
@@ -490,7 +494,7 @@ subroutine lsm_to_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
 !local pointers:
  character(len=StrKIND),pointer:: config_microp_scheme
 
- integer,dimension(:),pointer:: isltyp,ivgtyp
+ integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp
 
  real(kind=RKIND),dimension(:),pointer  :: acsnom,acsnow,canwat,chs,chs2,chklowq,cpm,cqs2,glw, &
                                            grdflx,gsw,hfx,lai,lh,noahres,potevp,qfx,qgh,qsfc,  &
@@ -544,6 +548,7 @@ subroutine lsm_to_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
  call mpas_pool_get_array(diag_physics,'znt'       ,znt       )
 
  call mpas_pool_get_array(sfc_input,'isltyp'    ,isltyp    )
+ call mpas_pool_get_array(sfc_input,'isctyp'    ,isctyp    )
  call mpas_pool_get_array(sfc_input,'ivgtyp'    ,ivgtyp    )
  call mpas_pool_get_array(sfc_input,'shdmin'    ,shdmin    )
  call mpas_pool_get_array(sfc_input,'shdmax'    ,shdmax    )
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F
index 7b93e7c..b2ba23c 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F
@@ -58,7 +58,7 @@ subroutine lsm_noahmp_fromMPAS(configs,mesh,diag,diag_physics,diag_physics_noahm
 
  integer:: i,its,ite
  integer:: n,ns,nsoil,nsnow,nzsnow
- integer,dimension(:),pointer:: isltyp,ivgtyp
+ integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp
 
  real(kind=RKIND),dimension(:),pointer:: latCell,lonCell
  real(kind=RKIND),dimension(:),pointer:: shdmax,shdmin,vegfra,tmn,xice,xland
diff --git a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F
index 07481fd..bce73fe 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F
@@ -91,7 +91,7 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input)
  character(len=StrKIND),pointer:: mminlu,mminsl
 
  integer,pointer:: nCells,nSoilLevels
- integer,dimension(:),pointer:: ivgtyp,isltyp
+ integer,dimension(:),pointer:: ivgtyp,isltyp,isctyp
 
  real(kind=RKIND),dimension(:),pointer:: snoalb,snow,snowh
  real(kind=RKIND),dimension(:,:),pointer:: tslb,smois,sh2o
@@ -118,6 +118,7 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input)
  call mpas_pool_get_dimension(mesh,'nSoilLevels',nSoilLevels)
 
  call mpas_pool_get_array(sfc_input,'isltyp', isltyp)
+ call mpas_pool_get_array(sfc_input,'isctyp', isctyp)
  call mpas_pool_get_array(sfc_input,'ivgtyp', ivgtyp)
  call mpas_pool_get_array(sfc_input,'sh2o'  , sh2o  )
  call mpas_pool_get_array(sfc_input,'smois' , smois )
@@ -149,6 +150,20 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input)
        call physics_error_fatal("module_sf_noahlsm.F: lsminit: out of range value "// &
                             "of ISLTYP. Is this field in the input?" )
 
+    ! Make sure all cells have reasonable soil colour indices.
+    errflag = 0
+    do iCell = 1, nCells
+       if(isctyp(iCell) < 1) then
+          errflag = 1
+          write(err_message,*) "module_sf_noahlsm.F: lsminit: out of range ISCTYP ", &
+                               iCell,isctyp(iCell),isltyp(iCell),ivgtyp(iCell)
+          call physics_message(err_message)
+       endif
+    end do
+    if (errflag == 1) &
+       call physics_error_fatal("module_sf_noahlsm.F: lsminit: out of range value "// &
+                                "of ISCTYP. Is this field in the input?" )
+
 !initializes soil liquid water content SH2O:
     do iCell = 1, nCells
        bx = bb(isltyp(iCell))
diff --git a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F
index 19ad287..63d7cc5 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F
@@ -234,7 +234,7 @@ subroutine noahmp_init(configs,mesh,clock,diag_physics,diag_physics_noahmp,outpu
  logical,pointer:: urban_physics
 
  integer,pointer:: nsoilcomps
- integer,dimension(:),pointer:: isltyp,ivgtyp
+ integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp
  integer,dimension(:),pointer:: isnowxy
  integer,dimension(:),pointer:: irnumsi,irnummi,irnumfi
  
@@ -322,11 +322,14 @@ subroutine noahmp_init(configs,mesh,clock,diag_physics,diag_physics_noahmp,outpu
 !--- initialization of time-invariant surface variables needed in subroutine NoahmpInitMain:
  call mpas_pool_get_array(sfc_input,'dzs'   ,dzs   )
  call mpas_pool_get_array(sfc_input,'isltyp',isltyp)
+ call mpas_pool_get_array(sfc_input,'isctyp',isctyp)
  call mpas_pool_get_array(sfc_input,'ivgtyp',ivgtyp)
 
  do i = its, ite
     mpas_noahmp%isltyp(i) = isltyp(i)
+    mpas_noahmp%isctyp(i) = isctyp(i)
     mpas_noahmp%ivgtyp(i) = ivgtyp(i)
+
  enddo
  do ns = 1, nsoil
     mpas_noahmp%dzs(ns) = dzs(ns,its)
diff --git a/src/core_atmosphere/physics/mpas_atmphys_update_surface.F b/src/core_atmosphere/physics/mpas_atmphys_update_surface.F
index 6e2057d..4d6a59e 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_update_surface.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_update_surface.F
@@ -151,7 +151,7 @@ subroutine physics_update_sst(dminfo,config_frac_seaice,mesh,sfc_input,diag_phys
 
  call mpas_pool_get_array(sfc_input,'isice'     ,isice     )
  call mpas_pool_get_array(sfc_input,'iswater'   ,iswater   )
- call mpas_pool_get_array(sfc_input,'isltyp'    ,isltyp    )    
+ call mpas_pool_get_array(sfc_input,'isltyp'    ,isltyp    )
  call mpas_pool_get_array(sfc_input,'ivgtyp'    ,ivgtyp    )
  call mpas_pool_get_array(sfc_input,'landmask'  ,landmask  )
  call mpas_pool_get_array(sfc_input,'vegfra'    ,vegfra    )
diff --git a/src/core_atmosphere/physics/mpas_atmphys_vars.F b/src/core_atmosphere/physics/mpas_atmphys_vars.F
index c6ebb80..2e8a26a 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_vars.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_vars.F
@@ -882,6 +882,7 @@ module mpas_atmphys_vars
     
  integer,dimension(:,:),allocatable:: &
     isltyp_p,         &!dominant soil type category                                                            [-]
+    isctyp_p,         &!dominant soil colour category                                                            [-]
     ivgtyp_p           !dominant vegetation category                                                           [-]
 
  real(kind=RKIND),dimension(:),allocatable:: &
diff --git a/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md b/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md
index ae45c39..6427e6b 100644
--- a/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md
+++ b/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md
@@ -400,7 +400,7 @@
 
   - but in the file listed in the namelist as: HRLDAS_SETUP_FILE = "
   - The initialization fields are: SNOW,CANWAT,TSK,TSLB,SMOIS
-  - This file also contains the static grid/domain information: XLAT,XLONG,TMN,HGT,SEAICE,MAPFAC_MX,MAPFAC_MY,SHDMAX,SHDMIN,XLAND,IVGTYP,ISLTYP,DZS,ZS
+  - This file also contains the static grid/domain information: XLAT,XLONG,TMN,HGT,SEAICE,MAPFAC_MX,MAPFAC_MY,SHDMAX,SHDMIN,XLAND,IVGTYP,ISLTYP,ISCTYP,DZS,ZS
   - This file can also contains some optional fields: LAI  
   - NOTE: a WRF input file can be used as a HRLDAS_SETUP_FILE
 
diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90
index 2de35ed..1cba671 100644
--- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90
@@ -62,7 +62,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO)
     ! config domain variable
     noahmp%config%domain%SurfaceType                 = 1
     noahmp%config%domain%NumSwRadBand                = 2
-    noahmp%config%domain%SoilColor                   = 4
+    noahmp%config%domain%SoilColor                   = NoahmpIO%ISCTYP(I)
     noahmp%config%domain%NumCropGrowStage            = 8
     noahmp%config%domain%FlagSoilProcess             = NoahmpIO%calculate_soil
     noahmp%config%domain%NumSoilTimeStep             = NoahmpIO%soil_update_steps
diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90
index 66458d6..49a709f 100644
--- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90
@@ -37,6 +37,7 @@ subroutine NoahmpIOVarFinalizeDefault(NoahmpIO)
     if ( allocated (NoahmpIO%zsoil)     ) deallocate ( NoahmpIO%zsoil               ) ! depth to soil interfaces [m]
     if ( allocated (NoahmpIO%ivgtyp)    ) deallocate ( NoahmpIO%ivgtyp              ) ! vegetation type
     if ( allocated (NoahmpIO%isltyp)    ) deallocate ( NoahmpIO%isltyp              ) ! soil type
+    if ( allocated (NoahmpIO%isctyp)    ) deallocate ( NoahmpIO%isctyp              ) ! soil colour class
     if ( allocated (NoahmpIO%vegfra)    ) deallocate ( NoahmpIO%vegfra              ) ! vegetation fraction []
     if ( allocated (NoahmpIO%tmn)       ) deallocate ( NoahmpIO%tmn                 ) ! deep soil temperature [K]
     if ( allocated (NoahmpIO%xland)     ) deallocate ( NoahmpIO%xland               ) ! =2 ocean; =1 land/seaice
diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90
index 4f3c3f4..cfd91bb 100644
--- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90
@@ -41,6 +41,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO)
     if ( .not. allocated (NoahmpIO%zsoil)     ) allocate ( NoahmpIO%zsoil      (1:nsoil            )      ) ! depth to soil interfaces [m]
     if ( .not. allocated (NoahmpIO%ivgtyp)    ) allocate ( NoahmpIO%ivgtyp     (its:ite            )      ) ! vegetation type
     if ( .not. allocated (NoahmpIO%isltyp)    ) allocate ( NoahmpIO%isltyp     (its:ite            )      ) ! soil type
+    if ( .not. allocated (NoahmpIO%isctyp)    ) allocate ( NoahmpIO%isctyp     (its:ite            )      ) ! soil colour class
     if ( .not. allocated (NoahmpIO%vegfra)    ) allocate ( NoahmpIO%vegfra     (its:ite            )      ) ! vegetation fraction []
     if ( .not. allocated (NoahmpIO%tmn)       ) allocate ( NoahmpIO%tmn        (its:ite            )      ) ! deep soil temperature [K]
     if ( .not. allocated (NoahmpIO%xland)     ) allocate ( NoahmpIO%xland      (its:ite            )      ) ! =2 ocean; =1 land/seaice
@@ -469,6 +470,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO)
     NoahmpIO%ice             = undefined_int
     NoahmpIO%ivgtyp          = undefined_int
     NoahmpIO%isltyp          = undefined_int
+    NoahmpIO%isctyp          = undefined_int
     NoahmpIO%isnowxy         = undefined_int
     NoahmpIO%coszen          = undefined_real
     NoahmpIO%xlat            = undefined_real
diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90
index 05a29d7..1852d83 100644
--- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90
@@ -68,6 +68,7 @@ module NoahmpIOVarType
     integer                                                ::  soil_update_steps   ! number of model time steps to update soil process
     integer,                allocatable, dimension(:)      ::  ivgtyp              ! vegetation type
     integer,                allocatable, dimension(:)      ::  isltyp              ! soil type
+    integer,                allocatable, dimension(:)      ::  isctyp              ! soil colour class
     real(kind=kind_noahmp), allocatable, dimension(:)      ::  coszen              ! cosine zenith angle
     real(kind=kind_noahmp), allocatable, dimension(:)      ::  xlat                ! latitude [rad]
     real(kind=kind_noahmp), allocatable, dimension(:,:)    ::  dz8w                ! thickness of atmo layers [m]
diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90
index eb01ceb..7339ba2 100644
--- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90
@@ -28,7 +28,7 @@ subroutine NoahmpReadTable(NoahmpIO)
     !-------------------------------------------------------
     integer, parameter :: MVT         = 27   ! number of vegetation types
     integer, parameter :: MBAND       = 2    ! number of radiation bands
-    integer, parameter :: MSC         = 8    ! number of soil texture
+    integer, parameter :: MSC         = 21   ! number of soil colour types
     integer, parameter :: MAX_SOILTYP = 30   ! max number of soil types
     integer, parameter :: NCROP       = 5    ! number of crop types
     integer, parameter :: NSTAGE      = 8    ! number of crop growth stages
diff --git a/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL b/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL
index c9d37c5..4ab3b66 100644
--- a/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL
+++ b/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL
@@ -405,13 +405,19 @@
 /
 
 &noahmp_rad_parameters
- !------------------------------------------------------------------------------
- ! soil color:    1       2       3       4       5       6       7       8     soil color index for soil albedo
- !------------------------------------------------------------------------------
- ALBSAT_VIS =   0.15,   0.11,   0.10,   0.09,   0.08,   0.07,   0.06,   0.05    ! saturated soil albedo at visible band
- ALBSAT_NIR =   0.30,   0.22,   0.20,   0.18,   0.16,   0.14,   0.12,   0.10    ! saturated soil albedo at NIR band
- ALBDRY_VIS =   0.27,   0.22,   0.20,   0.18,   0.16,   0.14,   0.12,   0.10    ! dry soil albedo at visible band
- ALBDRY_NIR =   0.54,   0.44,   0.40,   0.36,   0.32,   0.28,   0.24,   0.20    ! dry soil albedo at NIR band
+ !------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+ ! soil color index for soil albedo
+ ! soil color:    1       2       3       4       5       6       7       8       9      10      11      12      13      14      15      16      17      18      19      20      21   
+ !------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+ ! saturated soil albedo at visible band
+ ALBSAT_VIS =   0.25,   0.23,   0.21,   0.20,   0.19,   0.18,   0.17,   0.16,   0.15,   0.14,   0.13,   0.12,   0.11,   0.10,   0.09,   0.08,   0.07,   0.06,   0.05,   0.04,   0.00
+ ! saturated soil albedo at NIR band
+ ALBSAT_NIR =   0.50,   0.46,   0.42,   0.40,   0.38,   0.36,   0.34,   0.32,   0.30,   0.28,   0.26,   0.24,   0.22,   0.20,   0.18,   0.16,   0.14,   0.12,   0.10,   0.08,   0.00
+ ! dry soil albedo at visible band
+ ALBDRY_VIS =   0.36,   0.34,   0.32,   0.31,   0.30,   0.29,   0.28,   0.27,   0.26,   0.25,   0.24,   0.23,   0.22,   0.20,   0.18,   0.16,   0.14,   0.12,   0.10,   0.08,   0.00
+ !  dry soil albedo at NIR band
+ ALBDRY_NIR =   0.61,   0.57,   0.53,   0.51,   0.49,   0.48,   0.45,   0.43,   0.41,   0.39,   0.37,   0.35,   0.33,   0.31,   0.29,   0.27,   0.25,   0.23,   0.21,   0.16,   0.00
+ !------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  ALBICE     =   0.80,   0.55                                                    ! albedo land ice: 1=vis, 2=nir
  ALBLAK     =   0.60,   0.40                                                    ! albedo frozen lakes: 1=vis, 2=nir
  OMEGAS     =   0.8 ,   0.4                                                     ! two-stream parameter omega for snow
diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml
index ceec72e..14bb475 100644
--- a/src/core_init_atmosphere/Registry.xml
+++ b/src/core_init_atmosphere/Registry.xml
@@ -442,6 +442,7 @@
                         <var name="iswater_lu" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="ivgtyp" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="isltyp" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
+                        <var name="isctyp" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="snoalb" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="soiltemp" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="greenfrac" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
@@ -547,6 +548,7 @@
                         <var name="iswater_lu"/>
                         <var name="ivgtyp"/>
                         <var name="isltyp"/>
+                        <var name="isctyp"/>
                         <var name="snoalb"/>
                         <var name="soiltemp"/>
                         <var name="greenfrac"/>
@@ -847,6 +849,9 @@
                 <var name="isltyp" name_in_code="soilcat_top"  type="integer"  dimensions="nCells" units="unitless"
                      description="dominate soil category"/>
 
+                <var name="isctyp" name_in_code="soilcol_idx"  type="integer"  dimensions="nCells" units="unitless"
+                     description="dominant soil colour category"/>
+
                 <var name="snoalb" type="real" dimensions="nCells" units="unitless"
                      description="annual maximum snow albedo"/>
 
diff --git a/src/core_init_atmosphere/mpas_geotile_manager.F b/src/core_init_atmosphere/mpas_geotile_manager.F
index 04f9c60..b6b1663 100644
--- a/src/core_init_atmosphere/mpas_geotile_manager.F
+++ b/src/core_init_atmosphere/mpas_geotile_manager.F
@@ -98,6 +98,7 @@ module mpas_geotile_manager
     !>     * isice = 24
     !>     * isurban = 1
     !>     * isoilwater = 14
+    !>     * islcolwater = 21
     !
     !-----------------------------------------------------------------------
     function mpas_geotile_mgr_init(mgr, path) result(ierr)
@@ -119,7 +120,7 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr)
         integer, pointer :: tile_z_start, tile_z_end
         integer, pointer :: signed
         integer, pointer :: tile_bdr
-        integer, pointer :: iswater, islake, isice, isurban, isoilwater
+        integer, pointer :: iswater, islake, isice, isurban, isoilwater, islcolwater
         integer, pointer :: category_min, category_max
         integer :: err_level
         real (kind=RKIND), pointer :: dx ! Grid spacing in the x-direction
@@ -188,7 +189,7 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr)
 
         !
         ! If this is a categorical field, then check to see if it has category_max and category_min,
-        ! and then set the defaults of iswater, islake, isice, isurban and isoilwater
+        ! and then set the defaults of iswater, islake, isice, isurban, isoilwater and islcolwater
         !
         call mpas_pool_get_config(mgr % pool, 'type', fieldType)
         if (fieldType == 'categorical') then
@@ -219,12 +220,14 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr)
             isice => null()
             isurban => null()
             isoilwater => null()
+            islcolwater => null()
 
             call mpas_pool_get_config(mgr % pool, 'iswater', iswater)
             call mpas_pool_get_config(mgr % pool, 'islake', islake)
             call mpas_pool_get_config(mgr % pool, 'isice', isice)
             call mpas_pool_get_config(mgr % pool, 'isurban', isurban)
             call mpas_pool_get_config(mgr % pool, 'isoilwater', isoilwater)
+            call mpas_pool_get_config(mgr % pool, 'islcolwater', islcolwater)
 
             if (.not. associated(iswater)) then
                 call mpas_pool_add_config(mgr % pool, 'iswater', 16)
@@ -245,6 +248,10 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr)
             if (.not. associated(isoilwater)) then
                 call mpas_pool_add_config(mgr % pool, 'isoilwater', 14)
             endif
+
+            if (.not. associated(isoilwater)) then
+                call mpas_pool_add_config(mgr % pool, 'islcolwater', 21)
+            end if
         endif
 
         !
diff --git a/src/core_init_atmosphere/mpas_init_atm_static.F b/src/core_init_atmosphere/mpas_init_atm_static.F
index 0c54009..291f068 100644
--- a/src/core_init_atmosphere/mpas_init_atm_static.F
+++ b/src/core_init_atmosphere/mpas_init_atm_static.F
@@ -93,6 +93,7 @@ end subroutine interp_accumulation_function
  real (kind=RKIND) :: soilcomp_msgval = 255.0_RKIND   ! Modified later based on index file for soilcomp
  integer, dimension(:), pointer :: lu_index
  integer, dimension(:), pointer :: soilcat_top
+ integer, dimension(:), pointer :: soilcol_idx
  integer, dimension(:), pointer :: nhs
  integer, dimension(:,:), allocatable:: ncat
  ! Landmask is used by the accumulation function for maxsnoalb and soilcomp,
@@ -155,6 +156,7 @@ subroutine init_atm_static(mesh, dims, configs)
 
  integer, pointer :: isice_lu, iswater_lu
  integer :: iswater_soil
+ integer :: iswater_slcol
  integer, pointer :: nCells, nCellsSolve, nEdges, nVertices, maxEdges
  logical, pointer :: on_a_sphere
  real (kind=RKIND), pointer :: sphere_radius
@@ -185,6 +187,7 @@ subroutine init_atm_static(mesh, dims, configs)
  real (kind=RKIND), pointer :: missing_value
  integer, dimension(:), pointer :: lu_index
  integer, dimension(:), pointer :: soilcat_top
+ integer, dimension(:), pointer :: soilcol_idx
  integer, dimension(:), pointer :: landmask
  integer, dimension(:), pointer :: bdyMaskCell
  character(len=StrKIND), pointer :: mminlu
@@ -270,6 +273,7 @@ subroutine init_atm_static(mesh, dims, configs)
  call mpas_pool_get_array(mesh, 'isice_lu', isice_lu)
  call mpas_pool_get_array(mesh, 'iswater_lu', iswater_lu)
  call mpas_pool_get_array(mesh, 'soilcat_top', soilcat_top)
+ call mpas_pool_get_array(mesh, 'soilcol_idx', soilcol_idx)
  call mpas_pool_get_array(mesh, 'landmask', landmask)
  call mpas_pool_get_array(mesh, 'snoalb', snoalb)
  call mpas_pool_get_array(mesh, 'greenfrac', greenfrac)
@@ -443,17 +447,28 @@ subroutine init_atm_static(mesh, dims, configs)
                      supersample_fac=supersample_fac_30s)
  call mpas_log_write('--- end interpolate SOILCAT_TOP')
 
+!
+! Interpolate SOILCOL_IDX
+!
+ geog_sub_path = 'soilcolour_30s/'
+
+ call mpas_log_write('--- start interpolate SOILCOL_IDX')
+ call interp_soilcol(mesh, tree, trim(geog_data_path)//trim(geog_sub_path), iswater_slcol, &
+                     supersample_fac=supersample_fac_30s)
+ call mpas_log_write('--- end interpolate SOILCOL_IDX')
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! KLUDGE TO FIX SOIL TYPE OVER ANTARCTICA
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  where (lu_index == isice_lu) soilcat_top = 16
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! CORRECT INCONSISTENT SOIL AND LAND USE DATA
+! CORRECT INCONSISTENT SOIL TYPE, SOIL COLOUR AND LAND USE DATA
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  do iCell = 1,nCells
-    if (lu_index(iCell) == iswater_lu .or. &
-        soilcat_top(iCell) == iswater_soil) then
+    if (lu_index(iCell) == iswater_lu      .or. &
+        soilcat_top(iCell) == iswater_soil .or. &
+        soilcol_idx(iCell) == iswater_slcol ) then
         if (lu_index(iCell) /= iswater_lu) then
             call mpas_log_write('Turning lu_index into water at $i', intArgs=(/iCell/))
             lu_index(iCell) = iswater_lu
@@ -462,6 +477,10 @@ subroutine init_atm_static(mesh, dims, configs)
             call mpas_log_write('Turning soilcat_top into water at $i', intArgs=(/iCell/))
             soilcat_top(iCell) = iswater_soil
         end if
+        if (soilcol_idx(iCell) /= iswater_slcol) then
+            call mpas_log_write('Turning soilcol_idx into water at $i', intArgs=(/iCell/))
+            soilcol_idx(iCell) = iswater_slcol
+        end if
     end if
  end do
 
@@ -1936,6 +1955,89 @@ subroutine interp_soilcat(mesh, kdtree, geog_data_path, iswater_soil, supersampl
 
  end subroutine interp_soilcat
 
+ !***********************************************************************
+ !
+ !  routine interp_soilcol
+ !
+ !> \brief   Interpolate soil colour indices
+ !> \author  Marcos Longo
+ !> \date    15 October 2025
+ !> \details
+ !> This sub-routine is based on interp_soilcat. The procedure interpolates soil 
+ !>    colour indices by using the init_atm_map_static_data routine and then by
+ !>    accumulating the pixel values into each cell using 
+ !>    category_interp_accumulation.
+ !>
+ !> Variable mesh should be an mpas_pool that contains nCells and soilcol_idx;
+ !>     kdtree should be an initialised mpas_kd_type tree with dimensions 
+ !>    (xCell, yCell, zCell); and geog_data_path should be the path to the 
+ !>    soil colour index data set.  The values used by this data set to flag
+ !>    pixels with no data (e.g., pixels fully covered with water or ice).  Some
+ !>    data sets may not have the no data flag, but we harmonise this data set
+ !>    with the land use class mapping.
+ !-----------------------------------------------------------------------
+ subroutine interp_soilcol(mesh, kdtree, geog_data_path, iswater_slcol, supersample_fac)
+
+     implicit none
+
+     ! Input variables
+     type (mpas_pool_type), intent(inout) :: mesh
+     type (mpas_kd_type), pointer, intent(in) :: kdtree
+     character (len=*), intent(in) :: geog_data_path
+     integer, intent(out) :: iswater_slcol
+     integer, intent(in), optional :: supersample_fac
+
+     ! Local variables
+     type (mpas_geotile_mgr_type) :: mgr
+     integer, pointer :: nCells
+     integer, pointer :: iswater_slcol_ptr
+
+     real (kind=RKIND), pointer :: scalefactor
+
+     integer :: iCell
+     integer :: ierr
+
+     ierr = mgr % init(trim(geog_data_path))
+     if (ierr /= 0) then
+         call mpas_log_write("Error occured initalising interpolation for "//trim(geog_data_path), messageType=MPAS_LOG_CRIT)
+         return
+     end if
+
+     call mpas_pool_get_dimension(mesh, 'nCells', nCells)
+     call mpas_pool_get_array(mesh, 'soilcol_idx', soilcol_idx)
+     call mpas_pool_get_config(mgr % pool, 'scale_factor', scalefactor      )
+     call mpas_pool_get_config(mgr % pool, 'category_min', category_min     )
+     call mpas_pool_get_config(mgr % pool, 'category_max', category_max     )
+     call mpas_pool_get_config(mgr % pool, 'islcolwater' , iswater_slcol_ptr)
+
+     iswater_slcol = iswater_slcol_ptr
+
+     allocate(ncat(category_min:category_max, nCells))
+     ncat(:,:) = 0
+
+     call init_atm_map_static_data(mesh, mgr, kdtree, categorical_interp_criteria, categorical_interp_accumulation, &
+                                   supersample_fac=supersample_fac)
+
+     do iCell = 1, nCells
+        !     Variable ncat will tally the classes occurring in the grid cell. We select the class that has the 
+        ! highest count as the value representing the grid cell. We use maxloc for this, however, maxloc returns
+        ! the index corresponding to the maximum value starting from 1. We shift the result to make sure the 
+        ! category value correctly maps between category_min and category_max. 
+        soilcol_idx(iCell) = maxloc(ncat(:,iCell), dim=1) - 1 + category_min
+     end do
+     deallocate(ncat)
+
+     ierr = mgr % finalize()
+     if (ierr /= 0) then
+         call mpas_log_write("Error occured finalising interpolation for "//trim(geog_data_path), messageType=MPAS_LOG_CRIT)
+         return
+     end if
+
+     nullify(category_min)
+     nullify(category_max)
+
+ end subroutine interp_soilcol
+
  !***********************************************************************
  !
  !  routine derive_landmask
diff --git a/src/framework/mpas_io.F b/src/framework/mpas_io.F
index 09514a3..bf5a898 100644
--- a/src/framework/mpas_io.F
+++ b/src/framework/mpas_io.F
@@ -868,7 +868,7 @@ subroutine MPAS_io_inq_var(handle, fieldname, fieldtype, ndims, dimnames, dimsiz
             if (present(ierr)) ierr = MPAS_IO_ERR_BACKEND
             deallocate(new_fieldlist_node % fieldhandle)
             deallocate(new_fieldlist_node)
-            call mpas_log_write('Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN)
+            call mpas_log_write('PIO_SUPPORT: Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN)
             return
          end if
 !call mpas_log_write('Inquired about variable ID $i.', intArgs=(/new_fieldlist_node % fieldhandle % fieldid/) )
@@ -907,7 +907,7 @@ subroutine MPAS_io_inq_var(handle, fieldname, fieldtype, ndims, dimnames, dimsiz
             if (present(ierr)) ierr = MPAS_IO_ERR_BACKEND
             deallocate(new_fieldlist_node % fieldhandle)
             deallocate(new_fieldlist_node)
-            call mpas_log_write('Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN)
+            call mpas_log_write('SMIOL_SUPPORT_A: Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN)
             return
          end if
          new_fieldlist_node % fieldhandle % field_type = smiol_type
@@ -955,7 +955,7 @@ subroutine MPAS_io_inq_var(handle, fieldname, fieldtype, ndims, dimnames, dimsiz
 
             deallocate(new_fieldlist_node % fieldhandle)
             deallocate(new_fieldlist_node)
-            call mpas_log_write('Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN)
+            call mpas_log_write('SMIOL_SUPPORT_B: Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN)
             return
          end if
          new_fieldlist_node% fieldhandle % ndims = smiol_ndims