diff --git a/.github/actions/testing-setup/action.yml b/.github/actions/testing-setup/action.yml
index 8a3264b140..6ba149d927 100644
--- a/.github/actions/testing-setup/action.yml
+++ b/.github/actions/testing-setup/action.yml
@@ -5,10 +5,6 @@ inputs:
     description: 'If true, will build the symmetric MOM6 executable'
     required: false
     default: 'true'
-  install_python:
-    description: 'If true, will install the local python env needed for .testing'
-    required: false
-    default: 'true'
 runs:
   using: 'composite'
   steps:
@@ -54,14 +50,6 @@ runs:
         test ${{ inputs.build_symmetric }} == true && make build/symmetric/MOM6 -j
         echo "::endgroup::"
 
-    - name: Install local python venv for generating input data
-      shell: bash
-      run: |
-        echo "::group::Create local python env for input data generation"
-        cd .testing
-        test ${{ inputs.install_python }} == true && make work/local-env
-        echo "::endgroup::"
-
     - name: Set flags
       shell: bash
       run: |
diff --git a/.github/workflows/coupled-api.yml b/.github/workflows/coupled-api.yml
index 443755c7f4..2c9fa32720 100644
--- a/.github/workflows/coupled-api.yml
+++ b/.github/workflows/coupled-api.yml
@@ -20,7 +20,6 @@ jobs:
     - uses: ./.github/actions/testing-setup
       with:
         build_symmetric: 'false'
-        install_python: 'false'
 
     - name: Compile MOM6 for the GFDL coupled driver
       shell: bash
diff --git a/.testing/Makefile b/.testing/Makefile
index 530a552181..3e5c174239 100644
--- a/.testing/Makefile
+++ b/.testing/Makefile
@@ -11,7 +11,7 @@
 #       Delete the MOM6 test executables and dependency builds (FMS)
 #
 #   make clean.build
-#   	Delete only the MOM6 test executables
+#       Delete only the MOM6 test executables
 #
 #
 # Configuration:
@@ -204,34 +204,11 @@ endif
 FMS_SOURCE = $(call SOURCE,deps/fms/src)
 
 
-#---
-# Python preprocessing environment configuration
-
-HAS_NUMPY = $(shell python -c "import numpy" 2> /dev/null && echo "yes")
-HAS_NETCDF4 = $(shell python -c "import netCDF4" 2> /dev/null && echo "yes")
-
-USE_VENV =
-ifneq ($(HAS_NUMPY), yes)
-  USE_VENV = yes
-endif
-ifneq ($(HAS_NETCDF4), yes)
-  USE_VENV = yes
-endif
-
-# When disabled, activation is a null operation (`true`)
-VENV_PATH =
-VENV_ACTIVATE = true
-ifeq ($(USE_VENV), yes)
-  VENV_PATH = work/local-env
-  VENV_ACTIVATE = . $(VENV_PATH)/bin/activate
-endif
-
-
 #---
 # Rules
 
 .PHONY: all build.regressions build.prof
-all: $(foreach b,$(BUILDS),build/$(b)) $(VENV_PATH)
+all: $(foreach b,$(BUILDS),build/$(b))
 build.regressions: $(foreach b,symmetric target,build/$(b)/MOM6)
 build.prof: $(foreach b,opt opt_target,build/$(b)/MOM6)
 
@@ -382,8 +359,8 @@ deps/Makefile: ../ac/deps/Makefile
 # broken the ability to compile.  This is not a means to build a complete
 # coupled executable.
 # TODO:
-# 	- Avoid re-building FMS and MOM6 src by re-using existing object/mod files
-# 	- Use autoconf rather than mkmf templates
+#   - Avoid re-building FMS and MOM6 src by re-using existing object/mod files
+#   - Use autoconf rather than mkmf templates
 MK_TEMPLATE ?= ../../deps/mkmf/templates/ncrc-gnu.mk
 
 # NUOPC driver
@@ -402,21 +379,6 @@ build/mct/mom_ocean_model_mct.o: build/mct/Makefile
 check_mom6_api_mct: build/mct/mom_ocean_model_mct.o
 
 
-#---
-# Python preprocessing
-
-# NOTE: Some less mature environments (e.g. Arm64 Ubuntu) require explicit
-#   installation of numpy before netCDF4, as well as wheel and cython support.
-work/local-env:
-	python3 -m venv $@
-	. $@/bin/activate \
-	  && python3 -m pip install --upgrade pip \
-	  && pip3 install wheel \
-	  && pip3 install cython \
-	  && pip3 install numpy \
-	  && pip3 install netCDF4
-
-
 #---
 # Testing
 
@@ -555,6 +517,20 @@ $(foreach c,$(CONFIGS),$(eval $(call CONFIG_DIM_RULE,$(c))))
 	@echo -e "$(PASS): Diagnostics $*.regression.diag agree."
 
 
+#---
+# Preprocessing
+# NOTE: This only support tc4, but can be generalized over all tests.
+.PHONY: preproc
+preproc: tc4/Makefile
+	cd tc4 && make
+
+tc4/Makefile: tc4/configure tc4/Makefile.in
+	cd $(@D) && ./configure || (cat config.log && false)
+
+tc4/configure: tc4/configure.ac
+	cd $(@D) && autoreconf -if
+
+
 #---
 # Test run output files
 
@@ -567,17 +543,10 @@ $(foreach c,$(CONFIGS),$(eval $(call CONFIG_DIM_RULE,$(c))))
 # $(6): Number of MPI ranks
 
 define STAT_RULE
-work/%/$(1)/ocean.stats work/%/$(1)/chksum_diag: build/$(2)/MOM6 $(VENV_PATH)
+work/%/$(1)/ocean.stats work/%/$(1)/chksum_diag: build/$(2)/MOM6 | preproc
 	@echo "Running test $$*.$(1)..."
 	mkdir -p $$(@D)
 	cp -RL $$*/* $$(@D)
-	if [ -f $$(@D)/Makefile ]; then \
-	  $$(VENV_ACTIVATE) \
-	    && cd $$(@D) \
-	    && $(MAKE); \
-	else \
-	  cd $$(@D); \
-	fi
 	mkdir -p $$(@D)/RESTART
 	echo -e "$(4)" > $$(@D)/MOM_override
 	rm -f results/$$*/std.$(1).{out,err}
@@ -615,7 +584,7 @@ report.cov: run.cov codecov
 	  || { \
 	    cat build/cov/codecov.err ; \
 	    echo -e "${RED}Failed to upload report.${RESET}" ; \
-		if [ "$(REQUIRE_COVERAGE_UPLOAD)" = true ] ; then false ; fi ; \
+	    if [ "$(REQUIRE_COVERAGE_UPLOAD)" = true ] ; then false ; fi ; \
 	  }
 
 # Define $(,) as comma escape character
@@ -643,17 +612,10 @@ $(eval $(call STAT_RULE,cov,cov,true,,,1))
 #  2. Convert DAYMAX from TIMEUNIT to seconds
 #  3. Apply seconds to `ocean_solo_nml` inside input.nml.
 # NOTE: Assumes that runtime set by DAYMAX, will fail if set by input.nml
-work/%/restart/ocean.stats: build/symmetric/MOM6 $(VENV_PATH)
+work/%/restart/ocean.stats: build/symmetric/MOM6 | preproc
 	rm -rf $(@D)
 	mkdir -p $(@D)
 	cp -RL $*/* $(@D)
-	if [ -f $(@D)/Makefile ]; then \
-	  $(VENV_ACTIVATE) \
-	    && cd work/$*/restart \
-	    && $(MAKE); \
-	else \
-	  cd work/$*/restart; \
-	fi
 	mkdir -p $(@D)/RESTART
 	# Set the half-period
 	cd $(@D) \
@@ -754,7 +716,7 @@ report.cov.unit: build/unit/MOM_file_parser_tests.F90.gcov codecov
 	  || { \
 	    cat build/unit/codecov.err ; \
 	    echo -e "${RED}Failed to upload report.${RESET}" ; \
-		if [ "$(REQUIRE_COVERAGE_UPLOAD)" = true ] ; then false ; fi ; \
+	    if [ "$(REQUIRE_COVERAGE_UPLOAD)" = true ] ; then false ; fi ; \
 	  }
 
 
@@ -833,6 +795,13 @@ clean.build:
 
 
 .PHONY: clean.stats
-clean.stats:
+clean.stats: clean.preproc
 	@[ $$(basename $$(pwd)) = .testing ]
 	rm -rf work results
+
+
+.PHONY: clean.preproc
+clean.preproc:
+	@if [ -f tc4/Makefile ] ; then \
+	  cd tc4 && make clean ; \
+	fi
diff --git a/.testing/tc4/.gitignore b/.testing/tc4/.gitignore
index 29f62fb208..4f9cc2826f 100644
--- a/.testing/tc4/.gitignore
+++ b/.testing/tc4/.gitignore
@@ -1,4 +1,15 @@
+# Autoconf
+aclocal.m4
+autom4te.cache/
+config.log
+config.status
+configure~
+
+# Output
+gen_grid
 ocean_hgrid.nc
+topog.nc
+
+gen_data
 sponge.nc
 temp_salt_ic.nc
-topog.nc
diff --git a/.testing/tc4/Makefile b/.testing/tc4/Makefile
deleted file mode 100644
index a9aa395b9c..0000000000
--- a/.testing/tc4/Makefile
+++ /dev/null
@@ -1,8 +0,0 @@
-OUT=ocean_hgrid.nc sponge.nc temp_salt_ic.nc topog.nc
-
-$(OUT):
-	python build_grid.py
-	python build_data.py
-
-clean:
-	rm -rf $(OUT)
diff --git a/.testing/tc4/Makefile.in b/.testing/tc4/Makefile.in
new file mode 100644
index 0000000000..249d86b0b6
--- /dev/null
+++ b/.testing/tc4/Makefile.in
@@ -0,0 +1,38 @@
+FC = @FC@
+LD = @LD@
+FCFLAGS = @FCFLAGS@
+LDFLAGS = @LDFLAGS@
+LIBS = @LIBS@
+
+OUT = topog.nc ocean_hgrid.nc temp_salt_ic.nc sponge.nc
+
+all: $(OUT)
+
+ocean_hgrid.nc topog.nc: gen_grid
+	./gen_grid
+
+temp_salt_ic.nc sponge.nc: gen_data
+	./gen_data
+
+gen_grid: gen_grid.F90
+	$(FC) $(FCFLAGS) $(LDFLAGS) -o $@ $^ $(LIBS)
+
+gen_data: gen_data.F90
+	$(FC) $(FCFLAGS) $(LDFLAGS) -o $@ $^ $(LIBS)
+
+.PHONY: clean
+clean:
+	rm -rf $(OUT) gen_grid gen_data
+
+.PHONY: distclean
+distclean: clean
+	rm -f config.log
+	rm -f config.status
+	rm -f Makefile
+
+.PHONY: ac-clean
+ac-clean: distclean
+	rm -f aclocal.m4
+	rm -rf autom4te.cache
+	rm -f configure
+	rm -f configure~
diff --git a/.testing/tc4/build_data.py b/.testing/tc4/build_data.py
deleted file mode 100644
index e060d05cb1..0000000000
--- a/.testing/tc4/build_data.py
+++ /dev/null
@@ -1,80 +0,0 @@
-import netCDF4 as nc
-import numpy as np
-
-x = nc.Dataset('ocean_hgrid.nc').variables['x'][1::2, 1::2]
-y = nc.Dataset('ocean_hgrid.nc').variables['y'][1::2, 1::2]
-zbot = nc.Dataset('topog.nc').variables['depth'][:]
-zbot0 = zbot.max()
-
-
-def t_fc(x, y, z, radius=5.0, tmag=1.0):
-    """a radially symmetric anomaly in the center of the domain.
-    units are meters and degC.
-    """
-    ny, nx = x.shape
-    nz = z.shape[0]
-
-    x0 = x[int(ny/2), int(nx/2)]
-    y0 = y[int(ny/2), int(nx/2)]
-
-    tl = np.zeros((nz, ny, nx))
-    zb = z[-1]
-    if len(z) > 1:
-        zd = z / zb
-    else:
-        zd = [0.]
-    for k in np.arange(len(zd)):
-        r = np.sqrt((x - x0)**2 + (y - y0)**2)
-        tl[k, :] += (1.0 - np.minimum(r / radius, 1.0)) * tmag * (1.0 - zd[k])
-    return tl
-
-
-ny, nx = x.shape
-nz = 3
-z = (np.arange(nz) * zbot0) / nz
-
-temp = t_fc(x, y, z)
-salt = np.zeros(temp.shape)+35.0
-fl = nc.Dataset('temp_salt_ic.nc', 'w', format='NETCDF3_CLASSIC')
-fl.createDimension('lon', nx)
-fl.createDimension('lat', ny)
-fl.createDimension('depth', nz)
-fl.createDimension('Time', None)
-zv = fl.createVariable('depth', 'f8', ('depth'))
-lonv = fl.createVariable('lon', 'f8', ('lon'))
-latv = fl.createVariable('lat', 'f8', ('lat'))
-timev = fl.createVariable('Time', 'f8', ('Time'))
-timev.calendar = 'noleap'
-timev.units = 'days since 0001-01-01 00:00:00.0'
-timev.modulo = ' '
-tv = fl.createVariable('ptemp', 'f8', ('Time', 'depth', 'lat', 'lon'),
-                       fill_value=-1.e20)
-sv = fl.createVariable('salt', 'f8', ('Time', 'depth', 'lat', 'lon'),
-                       fill_value=-1.e20)
-tv[:] = temp[np.newaxis, :]
-sv[:] = salt[np.newaxis, :]
-zv[:] = z
-lonv[:] = x[0, :]
-latv[:] = y[:, 0]
-timev[0] = 0.
-fl.sync()
-fl.close()
-
-
-# Make Sponge forcing file
-dampTime = 20.0     # days
-secDays = 8.64e4
-fl = nc.Dataset('sponge.nc', 'w', format='NETCDF3_CLASSIC')
-fl.createDimension('lon', nx)
-fl.createDimension('lat', ny)
-lonv = fl.createVariable('lon', 'f8', ('lon'))
-latv = fl.createVariable('lat', 'f8', ('lat'))
-spv = fl.createVariable('Idamp', 'f8', ('lat', 'lon'), fill_value=-1.e20)
-Idamp = np.zeros((ny, nx))
-if dampTime > 0.:
-    Idamp = 0.0 + 1.0 / (dampTime * secDays)
-spv[:] = Idamp
-lonv[:] = x[0, :]
-latv[:] = y[:, 0]
-fl.sync()
-fl.close()
diff --git a/.testing/tc4/build_grid.py b/.testing/tc4/build_grid.py
deleted file mode 100644
index 7f1be74efd..0000000000
--- a/.testing/tc4/build_grid.py
+++ /dev/null
@@ -1,76 +0,0 @@
-import netCDF4 as nc
-from netCDF4 import stringtochar
-import numpy as np
-
-nx, ny = 14, 10     # Grid size
-depth0 = 100.       # Uniform depth
-ds = 0.01           # grid resolution at the equator in degrees
-Re = 6.378e6        # Radius of earth
-
-topo_ = np.zeros((ny, nx)) + depth0
-f_topo = nc.Dataset('topog.nc', 'w', format='NETCDF3_CLASSIC')
-ny, nx = topo_.shape
-f_topo.createDimension('ny', ny)
-f_topo.createDimension('nx', nx)
-f_topo.createDimension('ntiles', 1)
-f_topo.createVariable('depth', 'f8', ('ny', 'nx'))
-f_topo.createVariable('h2', 'f8', ('ny', 'nx'))
-f_topo.variables['depth'][:] = topo_
-f_topo.sync()
-f_topo.close()
-
-x_ = np.arange(0, 2*nx + 1) * ds  # units are degrees E
-y_ = np.arange(0, 2*ny + 1) * ds  # units are degrees N
-x, y = np.meshgrid(x_, y_)
-
-dx = np.zeros((2*ny + 1, 2*nx))
-dy = np.zeros((2*ny, 2*nx + 1))
-rad_deg = np.pi / 180.
-dx[:] = (rad_deg * Re * (x[:, 1:] - x[:, 0:-1])
-         * np.cos(0.5*rad_deg*(y[:, 0:-1] + y[:, 1:])))
-dy[:] = rad_deg * Re * (y[1:, :] - y[0:-1, :])
-
-f_sg = nc.Dataset('ocean_hgrid.nc', 'w', format='NETCDF3_CLASSIC')
-f_sg.createDimension('ny', 2*ny)
-f_sg.createDimension('nx', 2*nx)
-f_sg.createDimension('nyp', 2*ny + 1)
-f_sg.createDimension('nxp', 2*nx + 1)
-f_sg.createDimension('string', 5)
-f_sg.createVariable('y', 'f8', ('nyp', 'nxp'))
-f_sg.createVariable('x', 'f8', ('nyp', 'nxp'))
-dyv = f_sg.createVariable('dy', 'f8', ('ny', 'nxp'))
-dxv = f_sg.createVariable('dx', 'f8', ('nyp', 'nx'))
-areav = f_sg.createVariable('area', 'f8', ('ny', 'nx'))
-dxv.units = 'm'
-dyv.units = 'm'
-areav.units = 'm2'
-f_sg.createVariable('angle_dx', 'f8', ('nyp', 'nxp'))
-f_sg.createVariable('tile', 'S1', ('string'))
-f_sg.variables['y'].units = 'degrees'
-f_sg.variables['x'].units = 'degrees'
-f_sg.variables['dy'].units = 'meters'
-f_sg.variables['dx'].units = 'meters'
-f_sg.variables['area'].units = 'm2'
-f_sg.variables['angle_dx'].units = 'degrees'
-f_sg.variables['y'][:] = y
-f_sg.variables['x'][:] = x
-f_sg.variables['dx'][:] = dx
-f_sg.variables['dy'][:] = dy
-
-# Compute the area bounded by lines of constant
-# latitude-longitud on a sphere in m2.
-dlon = x_[1:] - x_[:-1]
-dlon = np.tile(dlon[np.newaxis, :], (2*ny, 1))
-y1_ = y_[:-1]
-y1_ = y1_[:, np.newaxis]*rad_deg
-y2_ = y_[1:]
-y2_ = y2_[:, np.newaxis]*rad_deg
-y1_ = np.tile(y1_, (1, 2*nx))
-y2_ = np.tile(y2_, (1, 2*nx))
-area = rad_deg * Re * Re * (np.sin(y2_) - np.sin(y1_)) * dlon
-f_sg.variables['area'][:] = area
-f_sg.variables['angle_dx'][:] = 0.
-str_ = stringtochar(np.array(['tile1'], dtype='S5'))
-f_sg.variables['tile'][:] = str_
-f_sg.sync()
-f_sg.close()
diff --git a/.testing/tc4/configure.ac b/.testing/tc4/configure.ac
new file mode 100644
index 0000000000..c431ad65ef
--- /dev/null
+++ b/.testing/tc4/configure.ac
@@ -0,0 +1,71 @@
+# tc4 preprocessor configuration
+AC_PREREQ([2.63])
+AC_INIT([], [])
+
+# Validate srdcir and configure input
+AC_CONFIG_SRCDIR([gen_grid.F90])
+AC_CONFIG_MACRO_DIR([../../ac/m4])
+
+
+# Explicitly assume free-form Fortran
+AC_LANG([Fortran])
+AC_FC_SRCEXT([f90])
+
+# We do not need MPI, but we want to emulate the executable used in MOM6
+AX_MPI([], [AC_MSG_ERROR([Could not find MPI launcher.])])
+AC_SUBST([FC], [$MPIFC])
+AC_SUBST([LD], [$MPILD])
+
+
+# netCDF configuration
+
+# Search for the Fortran netCDF module.
+AX_FC_CHECK_MODULE([netcdf], [], [
+  AS_UNSET([ax_fc_cv_mod_netcdf])
+  AC_PATH_PROG([NF_CONFIG], [nf-config])
+  AS_IF([test -n "$NF_CONFIG"], [
+    AC_SUBST([FCFLAGS], ["$FCFLAGS -I$($NF_CONFIG --includedir)"])
+    ], [AC_MSG_ERROR([Could not find nf-config.])]
+  )
+  AX_FC_CHECK_MODULE([netcdf], [], [
+    AC_MSG_ERROR([Could not find netcdf module.])
+  ])
+])
+
+# Confirm that the Fortran compiler can link the netCDF C library
+AX_FC_CHECK_C_LIB([netcdf], [nc_create], [], [
+  AS_UNSET([ax_fc_cv_c_lib_netcdf_nc_create])
+  AC_PATH_PROG([NC_CONFIG], [nc-config])
+  AS_IF([test -n "$NC_CONFIG"], [
+    AC_SUBST([LDFLAGS], ["$LDFLAGS -L$($NC_CONFIG --libdir)"])
+  ], [
+    AC_MSG_ERROR([Could not find nc-config.])
+  ])
+  AX_FC_CHECK_C_LIB([netcdf], [nc_create], [], [
+    AC_MSG_ERROR([Could not find netCDF C library.])
+  ])
+])
+
+# Confirm that the Fortran compiler can link to the netCDF Fortran library.
+# NOTE:
+#   - We test nf_create, rather than nf90_create, since AX_FC_CHECK_LIB can
+#     not currently probe the Fortran 90 interfaces.
+#   - nf-config does not have --libdir, so we parse the --flibs output.
+AX_FC_CHECK_LIB([netcdff], [nf_create], [], [], [
+  AS_UNSET([ax_fc_cv_lib_netcdff_nf_create])
+  AC_PATH_PROG([NF_CONFIG], [nf-config])
+  AS_IF([test -n "$NF_CONFIG"], [
+    AC_SUBST([LDFLAGS],
+      ["$LDFLAGS $($NF_CONFIG --flibs | xargs -n1 | grep "^-L" | sort -u | xargs)"]
+    )
+  ], [
+    AC_MSG_ERROR([Could not find nf-config.])
+  ])
+  AX_FC_CHECK_LIB([netcdff], [nf_create], [], [], [
+    AC_MSG_ERROR([Could not find netCDF Fortran library.])
+  ])
+])
+
+
+AC_CONFIG_FILES([Makefile])
+AC_OUTPUT
diff --git a/.testing/tc4/gen_data.F90 b/.testing/tc4/gen_data.F90
new file mode 100644
index 0000000000..8f44aa1465
--- /dev/null
+++ b/.testing/tc4/gen_data.F90
@@ -0,0 +1,189 @@
+use netcdf
+implicit none
+
+integer, parameter :: dp = selected_real_kind(10, 100)
+  !! Double precision (8-byte)
+
+integer, parameter :: nz = 3
+  !! Number of vertical layers
+real(kind=dp), parameter :: salt0 = 35._dp
+  !! Background salinity
+real(kind=dp), parameter :: dampTime = 20._dp
+  !! Sponge damping timescale [days]
+real(kind=dp), parameter :: secs_per_day = 86400._dp
+  !! Seconds per calendar day
+
+integer :: ncid
+
+integer :: x_id, y_id
+integer :: lon_dimid, lat_dimid, depth_dimid, time_dimid
+integer :: lon_id, lat_id, depth_id, time_id, temp_id, salt_id, idamp_id
+integer :: field_dimids(2)
+integer :: nx, ny
+
+integer :: i, rc
+
+real(kind=dp), allocatable :: x(:,:), y(:,:), z(:)
+  !! Temperature grid positions
+real(kind=dp), allocatable :: zbot(:,:)
+  !! Bottom topography
+real(kind=dp) :: zbot0
+  !! Maximum topographic depth
+real(kind=dp), allocatable :: temp(:,:,:), salt(:,:,:)
+  !! Initial temperature and salinity fields
+real(kind=dp), allocatable :: Idamp(:,:)
+  !! Sponge dampening rate
+
+! Read the domain grid
+rc = nf90_open('ocean_hgrid.nc', NF90_NOWRITE, ncid)
+
+rc = nf90_inq_varid(ncid, 'x', x_id)
+rc = nf90_inq_varid(ncid, 'y', y_id)
+
+rc = nf90_inquire_variable(ncid, x_id, dimids=field_dimids)
+rc = nf90_inquire_dimension(ncid, field_dimids(1), len=nx)
+rc = nf90_inquire_dimension(ncid, field_dimids(2), len=ny)
+
+! Extract center ("T") points of supergrid
+nx = nx / 2
+ny = ny / 2
+allocate(x(nx, ny), y(nx, ny))
+rc = nf90_get_var(ncid, x_id, x, start=[2,2], stride=[2,2])
+rc = nf90_get_var(ncid, y_id, y, start=[2,2], stride=[2,2])
+
+rc = nf90_close(ncid)
+
+
+! Read the topographic domain
+rc = nf90_open('topog.nc', NF90_NOWRITE, ncid)
+
+rc = nf90_inq_varid(ncid, 'depth', depth_id)
+rc = nf90_inquire_variable(ncid, depth_id, dimids=field_dimids)
+rc = nf90_inquire_dimension(ncid, field_dimids(1), len=nx)
+rc = nf90_inquire_dimension(ncid, field_dimids(2), len=ny)
+
+allocate(zbot(nx, ny))
+rc = nf90_get_var(ncid, depth_id, zbot)
+rc = nf90_close(ncid)
+
+
+! Construct the vertical axis
+allocate(z(nz))
+z = [(i, i=0,nz-1)] * maxval(zbot) / nz
+
+allocate(temp(nx, ny, nz), salt(nx, ny, nz))
+call t_fc(x, y, z, temp)
+salt(:,:,:) = salt0
+
+
+! Write T/S initial state
+rc = nf90_create('temp_salt_ic.nc', NF90_CLOBBER, ncid)
+
+rc = nf90_def_dim(ncid, 'lon', nx, lon_dimid)
+rc = nf90_def_dim(ncid, 'lat', ny, lat_dimid)
+rc = nf90_def_dim(ncid, 'depth', nz, depth_dimid)
+rc = nf90_def_dim(ncid, 'Time', NF90_UNLIMITED, time_dimid)
+
+rc = nf90_def_var(ncid, 'depth', NF90_DOUBLE, [depth_dimid], depth_id)
+rc = nf90_def_var(ncid, 'lon', NF90_DOUBLE, [lon_dimid], lon_id)
+rc = nf90_def_var(ncid, 'lat', NF90_DOUBLE, [lat_dimid], lat_id)
+rc = nf90_def_var(ncid, 'Time', NF90_DOUBLE, [time_dimid], time_id)
+
+rc = nf90_put_att(ncid, time_id, 'calendar', 'noleap')
+rc = nf90_put_att(ncid, time_id, 'units', 'days since 0001-01-01 00:00:00.0')
+! NOTE: nf90_put_att() truncates empty strings, so use nf90_put_att_any()
+rc = nf90_put_att_any(ncid, time_id, 'modulo', NF90_CHAR, 1, ' ')
+
+rc = nf90_def_var(ncid, 'ptemp', NF90_DOUBLE, &
+    [lon_dimid, lat_dimid, depth_dimid, time_dimid], temp_id)
+rc = nf90_def_var_fill(ncid, temp_id, 0, -1e20_dp)
+
+rc = nf90_def_var(ncid, 'salt', NF90_DOUBLE, &
+    [lon_dimid, lat_dimid, depth_dimid, time_dimid], salt_id)
+rc = nf90_def_var_fill(ncid, salt_id, 0, -1e20_dp)
+
+rc = nf90_enddef(ncid)
+
+rc = nf90_put_var(ncid, lon_id, x(:,1))
+rc = nf90_put_var(ncid, lat_id, y(1,:))
+rc = nf90_put_var(ncid, depth_id, z)
+rc = nf90_put_var(ncid, time_id, 0.)
+rc = nf90_put_var(ncid, temp_id, temp)
+rc = nf90_put_var(ncid, salt_id, salt)
+
+rc = nf90_close(ncid)
+
+
+! Sponge file
+rc = nf90_create('sponge.nc', NF90_CLOBBER, ncid)
+
+rc = nf90_def_dim(ncid, 'lon', nx, lon_dimid)
+rc = nf90_def_dim(ncid, 'lat', ny, lat_dimid)
+
+rc = nf90_def_var(ncid, 'lon', NF90_DOUBLE, lon_id)
+rc = nf90_def_var(ncid, 'lat', NF90_DOUBLE, lat_id)
+rc = nf90_def_var(ncid, 'Idamp', NF90_DOUBLE, [lon_dimid, lat_dimid], Idamp_id)
+rc = nf90_def_var_fill(ncid, Idamp_id, 0, -1e20_dp)
+
+rc = nf90_enddef(ncid)
+
+allocate(Idamp(nx, ny))
+Idamp = 0.
+if (dampTime > 0.) &
+  Idamp(:,:) = 1. / (dampTime * secs_per_day)
+
+rc = nf90_put_var(ncid, Idamp_id, Idamp)
+rc = nf90_put_var(ncid, lon_id, x(:,1))
+rc = nf90_put_var(ncid, lat_id, y(1,:))
+
+rc = nf90_close(ncid)
+
+contains
+
+subroutine t_fc(x, y, z, tl, radius, tmag)
+  real(kind=dp), intent(in) :: x(:,:), y(:,:), z(:)
+    !! Grid positions
+  real(kind=dp), intent(inout) :: tl(:,:,:)
+    !! Temperature field on the model grid
+  real(kind=dp), intent(in), optional :: radius
+    !! Temperature anomaly radius
+  real(kind=dp), intent(in), optional :: tmag
+    !! Temperature anomaly maximum
+
+  real(kind=dp) :: t_rad, t_max
+    !! Temperature field parameters (radius, max value)
+  real(kind=dp) :: x0, y0
+    !! Center of anomaly (currently midpoint of domain)
+  real(kind=dp), allocatable :: r(:,:), zd(:)
+    !! Radial and vertical extent of anomaly
+  integer :: k, nz
+    !! Vertical level indexing
+
+  t_rad = 5._dp
+  if (present(radius)) t_rad = radius
+
+  t_max = 1._dp
+  if (present(tmag)) t_max = tmag
+
+  ! Reduce supergrid size to T/S grid
+  allocate(zd, source=z)
+
+  x0 = x(1 + size(x, 1)/2, 1 + size(x, 2)/2)
+  y0 = y(1 + size(y, 1)/2, 1 + size(y, 2)/2)
+
+  tl(:,:,:) = 0.
+  nz = size(z)
+  if (nz > 1) then
+    zd(:) = z(:) / z(nz)
+  else
+    zd(:) = 0.
+  endif
+
+  allocate(r, source=x)
+  r(:,:) = hypot(x(:,:) - x0, y(:,:) - y0)
+  do k = 1, nz
+    tl(:,:,k) = (1. - min(r(:,:) / t_rad, 1.)) * t_max * (1. - zd(k))
+  enddo
+end subroutine t_fc
+
+end
diff --git a/.testing/tc4/gen_grid.F90 b/.testing/tc4/gen_grid.F90
new file mode 100644
index 0000000000..e76a681924
--- /dev/null
+++ b/.testing/tc4/gen_grid.F90
@@ -0,0 +1,108 @@
+use netcdf
+
+implicit none
+
+integer, parameter :: dp = selected_real_kind(10, 100)
+  !! Double precision (8-byte)
+
+integer, parameter :: nx = 14, ny = 10
+  !! Grid size
+real(kind=dp), parameter :: depth0 = 100._dp
+  !! Uniform depth
+real(kind=dp), parameter :: ds = 0.01_dp
+  !! Grid resolution at the equator in degrees
+real(kind=dp), parameter :: Re = 6.378e6_dp
+  !! Radius of earth
+real(kind=dp), parameter :: rad_per_deg = (4. * atan(1._dp)) / 180._dp
+  !! Degress to radians (= pi/180.)
+
+integer :: ncid
+integer :: nx_id, ny_id, nxp_id, nyp_id, ntile_id, string_id
+integer :: depth_id, h2_id
+integer :: x_id, y_id, dx_id, dy_id, area_id, angle_id, tile_id
+
+! Fields on model grid
+real(kind=dp) :: depth(nx, ny)
+
+! Grid fields (defined on supergrid)
+real(kind=dp) :: xg(0:2*nx), yg(0:2*ny)
+real(kind=dp) :: x(0:2*nx, 0:2*ny), y(0:2*nx, 0:2*ny)
+real(kind=dp) :: dx(0:2*nx-1, 0:2*ny)
+real(kind=dp) :: dy(0:2*nx, 0:2*ny-1)
+real(kind=dp) :: area(0:2*nx-1, 0:2*ny-1)
+real(kind=dp) :: angle_dx(0:2*nx, 0:2*ny)
+
+integer :: i, j, rc
+
+
+! Topography
+rc = nf90_create('topog.nc', NF90_CLOBBER, ncid)
+
+rc = nf90_def_dim(ncid, 'ny', ny, ny_id)
+rc = nf90_def_dim(ncid, 'nx', nx, nx_id)
+rc = nf90_def_dim(ncid, 'ntiles', 1, ntile_id)
+
+rc = nf90_def_var(ncid, 'depth', NF90_DOUBLE, [nx_id, ny_id], depth_id)
+rc = nf90_def_var(ncid, 'h2', NF90_DOUBLE, [nx_id, ny_id], h2_id)
+
+rc = nf90_enddef(ncid)
+
+depth(:,:) = depth0
+rc = nf90_put_var(ncid, depth_id, depth)
+
+rc = nf90_close(ncid)
+
+
+! Horizontal grid
+rc = nf90_create('ocean_hgrid.nc', NF90_CLOBBER, ncid)
+
+rc = nf90_def_dim(ncid, 'ny', 2*ny, ny_id)
+rc = nf90_def_dim(ncid, 'nx', 2*nx, nx_id)
+rc = nf90_def_dim(ncid, 'nyp', 2*ny+1, nyp_id)
+rc = nf90_def_dim(ncid, 'nxp', 2*nx+1, nxp_id)
+rc = nf90_def_dim(ncid, 'string', 5, string_id)
+
+rc = nf90_def_var(ncid, 'y', NF90_DOUBLE, [nxp_id, nyp_id], y_id)
+rc = nf90_def_var(ncid, 'x', NF90_DOUBLE, [nxp_id, nyp_id], x_id)
+rc = nf90_def_var(ncid, 'dy', NF90_DOUBLE, [nxp_id, ny_id], dy_id)
+rc = nf90_def_var(ncid, 'dx', NF90_DOUBLE, [nx_id, nyp_id], dx_id)
+rc = nf90_def_var(ncid, 'area', NF90_DOUBLE, [nx_id, ny_id], area_id)
+rc = nf90_def_var(ncid, 'angle_dx', NF90_DOUBLE, [nxp_id, nyp_id], angle_id)
+rc = nf90_def_var(ncid, 'tile', NF90_CHAR, string_id, tile_id)
+
+rc = nf90_put_att(ncid, y_id, 'units', 'degrees')
+rc = nf90_put_att(ncid, x_id, 'units', 'degrees')
+rc = nf90_put_att(ncid, dy_id, 'units', 'meters')
+rc = nf90_put_att(ncid, dx_id, 'units', 'meters')
+rc = nf90_put_att(ncid, area_id, 'units', 'm2')
+rc = nf90_put_att(ncid, angle_id, 'units', 'degrees')
+
+rc = nf90_enddef(ncid)
+
+xg = ds * [(i, i=0, 2*nx)]
+yg = ds * [(j, j=0, 2*ny)]
+
+! NOTE: sin() and cos() are compiler-dependent
+
+x(:,:) = spread(xg(:), 2, 2*ny+1)
+y(:,:) = spread(yg(:), 1, 2*nx+1)
+dx(:,:) = rad_per_deg * Re * (x(1:,:) - x(:2*nx-1,:)) &
+    * cos(0.5 * rad_per_deg * (y(1:,:) + y(:2*nx-1,:)))
+dy(:,:) = rad_per_deg * Re * (y(:,1:) - y(:,:2*ny-1))
+
+area(:,:) = rad_per_deg * Re * Re &
+  * spread(sin(rad_per_deg * yg(1:)) - sin(rad_per_deg * yg(:2*ny-1)), 1, 2*nx) &
+  * spread(xg(1:) - xg(:2*nx-1), 2, 2*ny)
+
+angle_dx(:,:) = 0.
+
+rc = nf90_put_var(ncid, x_id, x)
+rc = nf90_put_var(ncid, y_id, y)
+rc = nf90_put_var(ncid, dx_id, dx)
+rc = nf90_put_var(ncid, dy_id, dy)
+rc = nf90_put_var(ncid, area_id, area)
+rc = nf90_put_var(ncid, angle_id, angle_dx)
+rc = nf90_put_var(ncid, tile_id, 'tile1')
+
+rc = nf90_close(ncid)
+end