diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
new file mode 100644
index 00000000..1d5461ef
--- /dev/null
+++ b/.github/workflows/build.yml
@@ -0,0 +1,215 @@
+name: Build
+
+on: [push, pull_request]
+
+defaults:
+  run:
+    shell: bash
+
+jobs:
+  build:
+    name: ${{ matrix.compiler }}-${{ matrix.version }} (${{ matrix.os }})
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ macos-14, macos-15, macos-15-intel, macos-26, ubuntu-24.04 ]
+        compiler: [ gfortran ]
+        version: [ 13, 14, 15 ] # Julienne supports GFortran 13+
+        extra_flags: [ -g -O3 ]
+
+        include:
+
+          # --- LLVM flang coverage ---
+
+          - os: macos-14
+            compiler: flang
+            version: 21
+          - os: macos-15
+            compiler: flang
+            version: 21
+
+          # https://hub.docker.com/r/snowstep/llvm/tags
+          - os: ubuntu-24.04
+            compiler: flang
+            version: latest
+            container: snowstep/llvm:noble
+          - os: ubuntu-22.04
+            compiler: flang
+            version: latest
+            container: snowstep/llvm:jammy
+
+          # https://hub.docker.com/r/phhargrove/llvm-flang/tags
+          - os: ubuntu-24.04
+            compiler: flang
+            version: 21
+            network: smp
+            container: phhargrove/llvm-flang:21.1.0-latest
+          - os: ubuntu-24.04
+            compiler: flang
+            version: 20
+            container: phhargrove/llvm-flang:20.1.0-latest
+          - os: ubuntu-24.04
+            compiler: flang
+            version: 19
+            extra_flags: -g -mmlir -allow-assumed-rank -O3
+            container: phhargrove/llvm-flang:19.1.1-latest
+
+          # --- Intel coverage ---
+
+          # Intel provides two docker packages with Intel Fortran.
+          # The first is very compact, but lacks Intel C++:
+          #   https://hub.docker.com/r/intel/fortran-essentials/tags
+          #   This can be used with GNU C/C++
+          # The second includes Intel C++ and tons of other stuff,
+          # hence is much larger and can lead to slower container load:
+          #   https://hub.docker.com/r/intel/oneapi-hpckit/tags
+
+          # Julienne requires 2025.2.0 or later
+          - os: ubuntu-24.04
+            compiler: ifx-gcc
+            version: latest
+            container: intel/fortran-essentials:latest
+
+          - os: ubuntu-24.04
+            compiler: ifx
+            version: 2025.2.2
+            container: intel/oneapi-hpckit:2025.2.2-0-devel-ubuntu24.04
+
+          - os: ubuntu-24.04
+            compiler: ifx-gcc
+            version: 2025.2.0
+            container: intel/fortran-essentials:2025.2.0-0-devel-ubuntu24.04
+
+    container:
+      image: ${{ matrix.container }}
+
+    env:
+      COMPILER_VERSION: ${{ matrix.version }}
+      FC: ${{ matrix.compiler }}
+      FFLAGS: ${{ matrix.extra_flags }}
+      FPM_FLAGS:  --profile release --verbose
+
+    steps:
+    - name: Checkout code
+      uses: actions/checkout@v5
+
+    - name: Install Ubuntu Native Dependencies
+      if: ${{ contains(matrix.os, 'ubuntu') && matrix.container == '' && matrix.compiler == 'gfortran' }}
+      run: |
+        set -x
+        sudo apt-get update
+        sudo apt list -a 'gfortran-*'
+        sudo apt install -y build-essential
+        if [[ ${COMPILER_VERSION} < 15 ]] ; then \
+          sudo apt install -y gfortran-${COMPILER_VERSION} ; \
+        else \
+          curl -L https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh -o install-homebrew.sh ; \
+          chmod +x install-homebrew.sh ; \
+          env CI=1 ./install-homebrew.sh ; \
+          HOMEBREW_PREFIX="/home/linuxbrew/.linuxbrew" ; \
+          ${HOMEBREW_PREFIX}/bin/brew install -v gcc@${COMPILER_VERSION} binutils ; \
+          ls -al ${HOMEBREW_PREFIX}/bin ; \
+          echo "PATH=${HOMEBREW_PREFIX}/bin:${PATH}" >> "$GITHUB_ENV" ; \
+          : Homebrew GCC@15 needs binutils 2.44+ ; \
+          HOMEBREW_BINUTILS=$(ls -d ${HOMEBREW_PREFIX}/Cellar/binutils/2.*/bin ) ; \
+          ls -al ${HOMEBREW_BINUTILS} ; \
+          echo "FFLAGS=$FFLAGS -B ${HOMEBREW_BINUTILS}" >> "$GITHUB_ENV" ; \
+          echo "CFLAGS=$CFLAGS -B ${HOMEBREW_BINUTILS}" >> "$GITHUB_ENV" ; \
+          echo "CXXFLAGS=$CXXFLAGS -B ${HOMEBREW_BINUTILS}" >> "$GITHUB_ENV" ; \
+        fi
+
+    - name: Install Ubuntu Container Dependencies
+      if: ${{ contains(matrix.os, 'ubuntu') && matrix.container != '' && !contains(matrix.container, 'phhargrove') }}
+      run: |
+        set -x
+        apt update
+        apt install -y build-essential pkg-config make git curl
+
+    - name: Install macOS Dependencies
+      if: contains(matrix.os, 'macos')
+      run: |
+        set -x
+        brew update
+        # fpm binary distribution for macOS requires gfortran shared libraries from gcc@12
+        brew install gcc@12
+
+    - name: Install LLVM flang on macOS
+      if: contains(matrix.os, 'macos') && matrix.compiler == 'flang'
+      run: |
+        set -x
+        brew install llvm@${COMPILER_VERSION} flang
+        # workaround issue #228: clang cannot find homebrew flang's C header
+        for p in /opt/homebrew /usr/local $(brew --prefix) ; do find $p/Cellar/flang -name ISO_Fortran_binding.h 2>/dev/null || true ; done
+        echo "CFLAGS=-I$(dirname $(find $(brew --prefix)/Cellar/flang -name ISO_Fortran_binding.h | head -1)) ${CFLAGS}" >> "$GITHUB_ENV"
+        # Prepend homebrew clang to PATH:
+        echo "PATH=$(brew --prefix)/opt/llvm/bin:${PATH}" >> "$GITHUB_ENV"
+
+    - name: Setup Compilers
+      run: |
+        set -x
+        if test "$FC" = "flang" ; then \
+          echo "FPM_FC=flang-new" >> "$GITHUB_ENV" ; \
+        elif [[ "$FC" =~ "ifx" ]] ; then \
+          echo "FPM_FC=ifx" >> "$GITHUB_ENV" ; \
+          echo "FFLAGS=-fpp -coarray $FFLAGS" >> "$GITHUB_ENV" ; \
+          : ls -al /opt/intel/oneapi/compiler/2025.*/bin/ ; \
+          if type -p icpx ; then \
+            echo "FPM_CC=icx" >> "$GITHUB_ENV" ; \
+            echo "FPM_CXX=icpx" >> "$GITHUB_ENV" ; \
+          else \
+            echo "FPM_CC=gcc" >> "$GITHUB_ENV" ; \
+            echo "FPM_CXX=g++" >> "$GITHUB_ENV" ; \
+          fi ; \
+        elif test "$FC" = "lfortran" ; then \
+          echo "FPM_FC=lfortran" >> "$GITHUB_ENV" ; \
+          echo "FFLAGS=--cpp $FFLAGS" >> "$GITHUB_ENV" ; \
+        else \
+          echo "FPM_FC=gfortran-${COMPILER_VERSION}" >> "$GITHUB_ENV" ; \
+          echo "FFLAGS=-ffree-line-length-0 $FFLAGS" >> "$GITHUB_ENV" ; \
+          echo "FPM_CXXFLAGS=-std=c++14" >> "$GITHUB_ENV" ; : needed for Armadillo on macos-14 ; \
+        fi
+        if [[ "${{ matrix.container }}" =~ "snowstep/llvm" ]] ; then \
+          echo "LD_LIBRARY_PATH=/usr/lib/llvm-22/lib:$LD_LIBRARY_PATH" >> "$GITHUB_ENV" ; \
+        fi
+
+    - name: Setup FPM
+      uses: fortran-lang/setup-fpm@main
+      with:
+        github-token: ${{ secrets.GITHUB_TOKEN }}
+        fpm-version: latest
+
+    - name: Build FPM
+      if: false
+      run: |
+        set -x
+        curl --retry 5 -LOsS https://github.com/fortran-lang/fpm/releases/download/v0.11.0/fpm-0.11.0.F90
+        mkdir fpm-temp
+        gfortran-14 -o fpm-temp/fpm fpm-0.11.0.F90
+        echo "PATH=`pwd`/fpm-temp:${PATH}" >> "$GITHUB_ENV"
+
+    - name: Version info
+      run: |
+        if test -d .git && type -p git > /dev/null ; then \
+          git config --global --add safe.directory `pwd` ; \
+          echo == Commit Info == ; \
+          git log -n 1 ; echo ; echo ; \
+        fi
+        echo == Platform version info == 
+        uname -a
+        if test -r /etc/os-release ; then grep -e NAME -e VERSION /etc/os-release  ; fi
+        if test -x /usr/bin/sw_vers ; then /usr/bin/sw_vers ; fi
+        echo
+        echo PATH="$PATH"
+        for tool in ${FPM_FC} ${FPM_CC} ${FPM_CXX} fpm ; do
+          ( echo ; set -x ; w=$(which $tool) ; ls -al $w ; ls -alhL $w ; $tool --version )
+        done
+
+
+    - name: Build and Test
+      run: |
+        set -x
+        echo "FPM_FLAGS=${FPM_FLAGS}"
+        echo "FFLAGS=$FFLAGS"
+        ls -lt
+        fpm test ${FPM_FLAGS} --flag "$FFLAGS"
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 9e4d58a3..cb015c0c 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -75,7 +75,7 @@ jobs:
       - name: Install dependencies
         run: |
           brew update && (brew list cmake || brew install cmake)
-          brew install openblas eigen libomp lapack
+          brew install openblas eigen@3 libomp lapack
           sudo ln -sf /opt/homebrew/bin/gfortran-14 /opt/homebrew/bin/gfortran
         
       - name: Export Environment Variables
@@ -94,6 +94,8 @@ jobs:
           export CPPFLAGS+=" -I/opt/homebrew/opt/lapack/include"
           export PKG_CONFIG_PATH+=" /opt/homebrew/opt/lapack/lib/pkgconfig"
 
+          echo "Eigen3_DIR=/opt/homebrew/opt/eigen@3/share/eigen3/cmake" >> "$GITHUB_ENV"
+
       - name: Create build directory
         run: mkdir -p build
 
diff --git a/include/mole-language-support.F90 b/include/mole-language-support.F90
index 50171f15..9b332563 100644
--- a/include/mole-language-support.F90
+++ b/include/mole-language-support.F90
@@ -4,17 +4,12 @@
 #ifndef MOLE_LANGUAGE_SUPPORT
 #define MOLE_LANGUAGE_SUPPORT
 
-
 #ifdef __GNUC__
 #  define GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
 #else
 #  define GCC_VERSION 0
 #endif
 
-#if __GNUC__ && ( __GNUC__ < 14 || (__GNUC__ == 14 && __GNUC_MINOR__ < 3) )
-#define GCC_GE_MINIMUM
-#endif
-
 #ifndef HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT
 #  if defined(_CRAYFTN) || defined(__INTEL_COMPILER) || defined(NAGFOR) || defined(__flang__)
 #    define HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT 1
@@ -23,4 +18,12 @@
 #  endif
 #endif
 
+#ifndef HAVE_LOCALITY_SPECIFIER_SUPPORT
+#  if defined(NAGFOR) || defined(__flang__) || defined(__INTEL_COMPILER) || defined(_CRAYFTN) || (GCC_VERSION >= 150100)
+#    define HAVE_LOCALITY_SPECIFIER_SUPPORT 1
+#  else
+#    define HAVE_LOCALITY_SPECIFIER_SUPPORT 0
+#  endif
+#endif
+
 #endif
diff --git a/src/fortran/face_values_m.f90 b/src/fortran/face_values_m.f90
deleted file mode 100644
index 22999f89..00000000
--- a/src/fortran/face_values_m.f90
+++ /dev/null
@@ -1,114 +0,0 @@
-module face_values_m
-  !! Define an abstraction for face-centered values with a corresonding mimetic gradient operator
-  implicit none
-
-  private
-  public :: face_values_t
-  public :: gradient_t
-  public :: gradient_operator_t
-
-  type mimetic_matrix_t
-    !! Encapsulate a mimetic matrix with a corresponding matrix-vector product operator
-    private
-    double precision, allocatable :: upper_(:,:), inner_(:), lower_(:,:)
-  contains
-    generic :: operator(.x.) => matvec
-    procedure, private, non_overridable :: matvec
-  end type
-
-  type gradient_operator_t
-    !! Encapsulate kth-order mimetic gradient operator on dx-sized cells
-    private
-    integer k_
-    double precision dx_
-    type(mimetic_matrix_t) mimetic_matrix_
-  end type
-
-  type gradient_t
-    !! Encapsulate gradient values produced only by .grad. (private data, no constructors)
-    private
-    double precision, allocatable :: g_(:)
-  contains
-    procedure values
-  end type
-
-  interface 
-
-     pure module function values(self) result(gradients)
-       implicit none
-       class(gradient_t), intent(in) :: self
-       double precision, allocatable :: gradients(:)
-     end function
-
-  end interface
-
-  type face_values_t
-    !! Face-centered values 
-    private
-    double precision, allocatable :: f_(:)
-    type(gradient_operator_t) gradient_operator_
-  contains
-    generic :: operator(.grad.) => grad
-    procedure, non_overridable, private :: grad
-  end type
-
-  interface face_values_t
-
-    pure module function construct(f, k, dx) result(face_values)
-      !! Result is a collection of face-centered values with a mimetic gradient operator
-      implicit none
-      double precision, intent(in) :: f(:) !! face-centered values
-      double precision, intent(in) :: dx !! face spacing (cell width)
-      integer, intent(in) :: k !! order of accuracy
-      type(face_values_t) face_values
-    end function
-
-  end interface
-
-  interface
-
-    pure module function grad(self) result(grad_f)
-      !! Result is mimetic gradient of f
-      implicit none
-      class(face_values_t), intent(in) :: self
-      type(gradient_t) grad_f !! discrete gradient approximation
-    end function
-
-  end interface
-
-  interface gradient_operator_t
-
-    pure module function construct_from_parameters(k, dx) result(gradient_operator)
-      !! Construct a mimetic gradient operator
-      implicit none
-      integer, intent(in) :: k !! order of accuracy
-      double precision, intent(in) :: dx !! step size
-      type(gradient_operator_t) gradient_operator
-    end function
-
-  end interface
-
-  interface mimetic_matrix_t
-
-    pure module function construct_from_components(upper, inner, lower) result(mimetic_matrix)
-      !! Construct discrete operator from coefficient matrix
-      implicit none
-      double precision, intent(in) :: upper(:,:), inner(:), lower(:,:)
-      type(mimetic_matrix_t) mimetic_matrix
-    end function
-
-  end interface
-
-  interface
-
-    pure module function matvec(self, vector) result(gradient)
-      !! Apply a matrix operator to a vector
-      implicit none
-      class(mimetic_matrix_t), intent(in) :: self
-      type(face_values_t), intent(in) :: vector
-      type(gradient_t) gradient
-    end function
-
-  end interface
-
-end module face_values_m
\ No newline at end of file
diff --git a/src/fortran/face_values_s.F90 b/src/fortran/face_values_s.F90
deleted file mode 100644
index cc4a0a59..00000000
--- a/src/fortran/face_values_s.F90
+++ /dev/null
@@ -1,15 +0,0 @@
-submodule(face_values_m) face_values_s
-  implicit none
-
-contains
-
-  module procedure construct
-    face_values%f_ = f
-    face_values%gradient_operator_ = gradient_operator_t(k, dx)
-  end procedure
-
-  module procedure grad
-    grad_f = self%gradient_operator_%mimetic_matrix_ .x. self
-  end procedure
-
-end submodule face_values_s
\ No newline at end of file
diff --git a/src/fortran/gradient_1D_m.f90 b/src/fortran/gradient_1D_m.f90
new file mode 100644
index 00000000..62ed6a80
--- /dev/null
+++ b/src/fortran/gradient_1D_m.f90
@@ -0,0 +1,49 @@
+module gradient_1D_m
+  !! Define an abstraction for the collection of points used to compute gradidents:
+  !! cell centers plus oundaries.
+  implicit none
+
+  private
+  public :: gradient_1D_t
+
+  type gradient_1D_t
+    !! Encapsulate gradient_1D values produced only by .grad. (no other constructors)
+    private
+    double precision, allocatable :: vector_1D_(:) !! gradient_1D values at cell faces (nodes in 1D)
+    double precision x_min_ !! domain lower boundary
+    double precision x_max_ !! domain upper boundary
+    integer cells_ !! number of grid cells spanning the domain
+  contains
+    procedure values
+    procedure faces
+  end type
+
+  interface gradient_1D_t
+
+      pure module function construct_from_components(face_centered_values, x_min, x_max, cells) result(gradient_1D)
+        !! Result is an object storing gradient_1Ds at cell faces
+        implicit none
+        double precision, intent(in) :: face_centered_values(:), x_min, x_max
+        integer, intent(in) :: cells
+        type(gradient_1D_t) gradient_1D
+      end function
+
+  end interface
+
+  interface
+
+    pure module function faces(self) result(x)
+      implicit none
+      class(gradient_1D_t), intent(in) :: self
+      double precision, allocatable :: x(:)
+    end function
+
+     pure module function values(self) result(gradients)
+       implicit none
+       class(gradient_1D_t), intent(in) :: self
+       double precision, allocatable :: gradients(:)
+     end function
+
+  end interface
+
+end module gradient_1D_m
\ No newline at end of file
diff --git a/src/fortran/gradient_1D_s.F90 b/src/fortran/gradient_1D_s.F90
new file mode 100644
index 00000000..9cff5f90
--- /dev/null
+++ b/src/fortran/gradient_1D_s.F90
@@ -0,0 +1,25 @@
+submodule(gradient_1D_m) gradient_1D_s
+  implicit none
+
+contains
+
+  module procedure construct_from_components
+    gradient_1D%vector_1D_ = face_centered_values
+    gradient_1D%x_min_ = x_min
+    gradient_1D%x_max_ = x_max
+    gradient_1D%cells_ = cells
+  end procedure
+
+  module procedure values
+    gradients =  self%vector_1D_
+  end procedure
+
+  module procedure faces
+    integer cell
+    x = [ self%x_min_ &
+         ,self%x_min_ + [(cell*(self%x_max_ - self%x_min_)/self%cells_, cell = 1, self%cells_-1)] &
+         ,self%x_max_ &
+        ]
+  end procedure
+
+end submodule gradient_1D_s
\ No newline at end of file
diff --git a/src/fortran/gradient_operator_1D_s.F90 b/src/fortran/gradient_operator_1D_s.F90
new file mode 100644
index 00000000..f5f41ab2
--- /dev/null
+++ b/src/fortran/gradient_operator_1D_s.F90
@@ -0,0 +1,111 @@
+#include "julienne-assert-macros.h"
+#include "mole-language-support.F90"
+
+submodule(scalar_1D_m) gradient_operator_1D_s
+  use julienne_m, only : call_julienne_assert_, string_t
+#if ASSERTIONS
+  use julienne_m, only : operator(.isAtLeast.)
+#endif
+  implicit none
+
+contains
+
+  module procedure construct_from_parameters
+
+    call_julienne_assert(m .isAtLeast. 2*k)
+
+    gradient_operator_1D%mimetic_matrix_1D_ = mimetic_matrix_1D_t( &
+       corbino_castillo_A( k, dx) & 
+      ,corbino_castillo_M( k, dx) & 
+      ,corbino_castillo_Ap(k, dx) &
+    )
+    gradient_operator_1D%k_  = k
+    gradient_operator_1D%dx_ = dx
+    gradient_operator_1D%m_ = m
+  end procedure
+
+  pure function corbino_castillo_A(k, dx) result(matrix_block)
+    integer, intent(in) :: k
+    double precision, intent(in) :: dx
+    double precision, allocatable :: matrix_block(:,:)
+
+    order_of_accuracy: &
+    select case(k)
+    case(2)
+      matrix_block = reshape([-8D0/3D0, 3D0, -1D0/3D0] , shape=[1,3]) / dx
+    case(4)
+      matrix_block = reshape([ &
+         -352D0/105D0,  35D0/ 8D0, -35D0/24D0, 21D0/40D0, -5D0/ 56D0 &
+        ,  16D0/105D0, -31D0/24D0,  29D0/24D0, -3D0/40D0,  1D0/168D0 &
+      ], shape=[2,5], order=[2,1]) / dx
+    case default
+      associate(string_k => string_t(k))
+        error stop "corbino_castillo_A: unsupported order of accuracy: " // string_k%string()
+      end associate
+    end select order_of_accuracy
+
+  end function
+
+  pure function corbino_castillo_M(k, dx) result(row)
+    integer, intent(in) :: k
+    double precision, intent(in) :: dx
+    double precision, allocatable :: row(:)
+
+    order_of_accuracy: &
+    select case(k)
+    case(2)
+      row = [-1D0, 1D0]/ dx        
+    case(4)
+      row = [1D0/24D0, -9D0/8D0, 9D0/8D0, -1D0/24D0] / dx        
+    case default
+      associate(string_k => string_t(k))
+        error stop "corbino_castillo_A: unsupported order of accuracy: " // string_k%string()
+      end associate
+    end select order_of_accuracy
+
+  end function
+
+#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT
+
+  pure function corbino_castillo_Ap(k, dx) result(matrix_block)
+    integer, intent(in) :: k
+    double precision, intent(in) :: dx
+    double precision, allocatable :: matrix_block(:,:)
+
+    associate(A => corbino_castillo_A(k, dx))
+      allocate(matrix_block , mold=A)
+      reverse_elements_within_rows_and_flip_sign: &
+      do concurrent(integer :: row = 1:size(matrix_block,1)) default(none) shared(matrix_block, A)
+        matrix_block(row,:) = -A(row,size(A,2):1:-1)
+      end do reverse_elements_within_rows_and_flip_sign
+      reverse_elements_within_columns: &
+      do concurrent(integer :: column = 1 : size(matrix_block,2)) default(none) shared(matrix_block)
+        matrix_block(:,column) = matrix_block(size(matrix_block,1):1:-1,column)
+      end do reverse_elements_within_columns
+    end associate
+  end function
+ 
+#else
+
+  pure function corbino_castillo_Ap(k, dx) result(matrix_block)
+    integer, intent(in) :: k
+    double precision, intent(in) :: dx
+    double precision, allocatable :: matrix_block(:,:)
+    integer row, column
+
+    associate(A => corbino_castillo_A(k, dx))
+      allocate(matrix_block , mold=A)
+      reverse_elements_within_rows_and_flip_sign: &
+      do concurrent(row = 1:size(matrix_block,1))
+        matrix_block(row,:) = -A(row,size(A,2):1:-1)
+      end do reverse_elements_within_rows_and_flip_sign
+      reverse_elements_within_columns: &
+      do concurrent(column = 1 : size(matrix_block,2))
+        matrix_block(:,column) = matrix_block(size(matrix_block,1):1:-1,column)
+      end do reverse_elements_within_columns
+    end associate
+  end function
+ 
+#endif
+
+end submodule gradient_operator_1D_s
\ No newline at end of file
diff --git a/src/fortran/gradient_operator_s.F90 b/src/fortran/gradient_operator_s.F90
deleted file mode 100644
index 36053a35..00000000
--- a/src/fortran/gradient_operator_s.F90
+++ /dev/null
@@ -1,104 +0,0 @@
-#include "julienne-assert-macros.h"
-#include "mole-language-support.F90"
-
-submodule(face_values_m) gradient_operator_s
-  use julienne_m, only : call_julienne_assert_, string_t
-#if ASSERTIONS
-  use julienne_m, only : operator(.isAtLeast.)
-#endif
-  implicit none
-
-contains
-
-  module procedure construct_from_parameters
-
-    call_julienne_assert(m .isAtLeast. 2*k)
-
-    gradient_operator%mimetic_matrix_ = mimetic_matrix_t( &
-       corbino_castillo_A( k, dx) & 
-      ,corbino_castillo_M( k, dx) & 
-      ,corbino_castillo_Ap(k, dx) &
-    )
-    gradient_operator%k_  = k
-    gradient_operator%dx_ = dx
-  end procedure
-
-  pure function corbino_castillo_A(k, dx) result(rows)
-    integer, intent(in) :: k
-    double precision, intent(in) :: dx
-    double precision, allocatable :: rows(:,:)
-
-    order_of_accuracy: &
-    select case(k)
-    case(2)
-      rows = reshape([-8D0/3D0, 3D0, -1D0/3D0] , shape=[1,3]) / dx
-    case(4)
-      rows = transpose(reshape([ & 
-         [-352D0/105D0,  35D0/ 8D0, -35D0/24D0, -21D0/40D0, -5D0/ 56D0] &
-        ,[ -16D0/105D0, -31D0/24D0,  29D0/24D0,  -3D0/40D0,  1D0/168D0] &
-      ], shape=[5,2]))
-    case default
-      associate(string_k => string_t(k))
-        error stop "corbino_castillo_A: unsupported order of accuracy: " // string_k%string()
-      end associate
-    end select order_of_accuracy
-
-  end function
-
-  pure function corbino_castillo_M(k, dx) result(row)
-    integer, intent(in) :: k
-    double precision, intent(in) :: dx
-    double precision, allocatable :: row(:)
-
-    order_of_accuracy: &
-    select case(k)
-    case(2)
-      row = [-1D0, 1D0]/ dx        
-    case(4)
-      row = [1D0/24D0, -9D0/8D0, 9D0/8D0, -1D0/24D0] / dx        
-    case default
-      associate(string_k => string_t(k))
-        error stop "corbino_castillo_A: unsupported order of accuracy: " // string_k%string()
-      end associate
-    end select order_of_accuracy
-
-  end function
-
-#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT
-
-  pure function corbino_castillo_Ap(k, dx) result(rows)
-    integer, intent(in) :: k
-    double precision, intent(in) :: dx
-    double precision, allocatable :: rows(:,:)
-
-    associate(A => corbino_castillo_A(k, dx))
-      allocate(rows , mold=A)
-      reverse_and_flip_sign: &
-      do concurrent(integer :: row = 1:size(rows,1)) default(none) shared(rows, A) 
-        rows(row,:) = -A(row,size(A,2):1:-1)
-      end do reverse_and_flip_sign
-    end associate
-  end function
-
-#else
-
-  pure function corbino_castillo_Ap(k, dx) result(rows)
-    integer, intent(in) :: k
-    double precision, intent(in) :: dx
-    double precision, allocatable :: rows(:,:)
-
-    associate(A => corbino_castillo_A(k, dx))
-      allocate(rows , mold=A)
-      block
-        integer row
-        reverse_and_flip_sign: &
-        do concurrent(row = 1:size(rows,1)) default(none) shared(rows, A) 
-          rows(row,:) = -A(row,size(A,2):1)
-        end do reverse_and_flip_sign
-      end block
-    end associate
-  end function
-
-#endif
-
-end submodule gradient_operator_s
\ No newline at end of file
diff --git a/src/fortran/gradient_s.F90 b/src/fortran/gradient_s.F90
deleted file mode 100644
index fc289f04..00000000
--- a/src/fortran/gradient_s.F90
+++ /dev/null
@@ -1,10 +0,0 @@
-submodule(face_values_m) gradient_s
-  implicit none
-
-contains
-
-  module procedure values
-    gradients =  self%g_
-  end procedure
-
-end submodule gradient_s
\ No newline at end of file
diff --git a/src/fortran/mimetic_matrix_1D_s.F90 b/src/fortran/mimetic_matrix_1D_s.F90
new file mode 100644
index 00000000..1138d8bc
--- /dev/null
+++ b/src/fortran/mimetic_matrix_1D_s.F90
@@ -0,0 +1,87 @@
+#include "mole-language-support.F90"
+#include "julienne-assert-macros.h"
+
+submodule(scalar_1D_m) mimetic_matrix_1D_s
+  use julienne_m, only : call_julienne_assert_, string_t, operator(.equalsExpected.), operator(.csv.)
+  implicit none
+
+contains
+
+  module procedure construct_from_components
+    mimetic_matrix_1D%upper_ = upper
+    mimetic_matrix_1D%inner_ = inner
+    mimetic_matrix_1D%lower_ = lower
+  end procedure
+
+#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT && HAVE_LOCALITY_SPECIFIER_SUPPORT
+
+  module procedure matvec
+
+    double precision, allocatable :: product_inner(:)
+
+    associate(upper => size(self%upper_,1), lower => size(self%lower_,1))
+      associate(inner_rows => size(scalar_1D%scalar_1D_) - (upper + lower + 1))
+
+        allocate(product_inner(inner_rows))
+
+        do concurrent(integer :: row = 1 : inner_rows) default(none) shared(product_inner, self, scalar_1D)
+          product_inner(row) = dot_product(self%inner_, scalar_1D%scalar_1D_(row + 1 : row + size(self%inner_)))
+        end do
+
+        matvec_product = [ &
+           matmul(self%upper_, scalar_1D%scalar_1D_(1 : size(self%upper_,2))) &
+          ,product_inner &
+          ,matmul(self%lower_, scalar_1D%scalar_1D_(size(scalar_1D%scalar_1D_) - size(self%lower_,2) + 1 : )) &
+        ]
+      end associate
+    end associate
+  end procedure
+
+#else
+
+  module procedure matvec
+
+    integer row
+    double precision, allocatable :: product_inner(:)
+
+    associate(upper => size(self%upper_,1), lower => size(self%lower_,1))
+      associate(inner_rows => size(scalar_1D%scalar_1D_) - (upper + lower + 1))
+
+        allocate(product_inner(inner_rows))
+
+        do concurrent(row = 1 : inner_rows)
+          product_inner(row) = dot_product(self%inner_, scalar_1D%scalar_1D_(row + 1 : row + size(self%inner_)))
+        end do
+
+        matvec_product = [ &
+           matmul(self%upper_, scalar_1D%scalar_1D_(1 : size(self%upper_,2))) &
+          ,product_inner &
+          ,matmul(self%lower_, scalar_1D%scalar_1D_(size(scalar_1D%scalar_1D_) - size(self%lower_,2) + 1 : )) &
+        ]
+      end associate
+    end associate
+  end procedure
+
+#endif
+
+  module procedure to_file_t
+    type(string_t), allocatable :: lines(:)
+    integer, parameter :: inner_rows = 1
+    integer row
+
+    associate(upper_rows => size(self%upper_,1), lower_rows => size(self%lower_,1))
+      allocate(lines(upper_rows + inner_rows + lower_rows))
+      do row = 1, upper_rows
+        lines(row) = .csv. string_t(self%upper_(row,:))
+      end do
+      lines(upper_rows + inner_rows) = .csv. string_t(self%inner_)
+      do row = 1, lower_rows
+        lines(upper_rows + inner_rows + row) = .csv. string_t(self%lower_(row,:))
+      end do
+    end associate
+
+    file = file_t(lines)
+
+  end procedure
+
+end submodule mimetic_matrix_1D_s
\ No newline at end of file
diff --git a/src/fortran/mimetic_matrix_s.F90 b/src/fortran/mimetic_matrix_s.F90
deleted file mode 100644
index ebe8165e..00000000
--- a/src/fortran/mimetic_matrix_s.F90
+++ /dev/null
@@ -1,74 +0,0 @@
-#include "mole-language-support.F90"
-#include "julienne-assert-macros.h"
-
-submodule(face_values_m) mimetic_matrix_s
-  use julienne_m, only : call_julienne_assert_
-  implicit none
-
-contains
-
-  module procedure construct_from_components
-    mimetic_matrix%upper_ = upper
-    mimetic_matrix%inner_ = inner
-    mimetic_matrix%lower_ = lower
-  end procedure
-
-#if HAVE_DO_CONCURRENT_TYPE_SPEC_SUPPORT
-
-  module procedure matvec
-
-    double precision, allocatable :: product_inner(:)
-
-    associate(upper_rows  => size(self%upper_,1))
-      associate(inner_rows => (size(vector%f_) - 1) - 2*upper_rows) ! inner_rows = matrix rows - (upper rows + lower rows)
-
-        allocate(product_inner(inner_rows))
-
-        associate(inner_bandwidth => size(self%inner_))
-          do concurrent(integer :: row = 1 : inner_rows) default(none) &
-            shared(product_inner, self, vector, upper_rows, inner_bandwidth)
-            product_inner(row) = dot_product(self%inner_, vector%f_(upper_rows + row : upper_rows + row + inner_bandwidth - 1))
-          end do
-        end associate
-
-        associate(upper_bandwidth => size(self%upper_,2), lower_bandwidth => size(self%lower_,2))
-          gradient%g_ = [ &
-             matmul(self%upper_, vector%f_(1 : upper_bandwidth)) &
-            ,product_inner &
-            ,matmul(self%lower_, vector%f_(size(vector%f_) - lower_bandwidth + 1 : )) &
-          ]
-        end associate
-      end associate
-    end associate
-  end procedure
-
-#else
-
-  module procedure matvec
-
-    double precision, allocatable :: product_inner(:)
-
-    associate(upper_rows  => size(self%upper_,1), inner_bandwidth => size(self%inner_))
-      associate(inner_rows => (size(vector%f_) - 1) - 2*upper_rows) ! inner_rows = matrix rows - (upper rows + lower rows)
-
-        allocate(product_inner(inner_rows))
-
-        block
-          integer row
-          do concurrent(row = 1 : inner_rows) default(none) shared(product_inner, self, vector, upper_rows, inner_bandwidth)
-            product_inner(row) = dot_product(self%inner_, vector%f_(upper_rows + row : upper_rows + inner_bandwidth))
-          end do
-        end block
-
-        gradient%g_ = [ &
-           matmul(self%upper_, vector%f_(1 : upper_rows)) &
-          ,product_inner &
-          ,matmul(self%lower_, vector%f_(upper_rows + inner_rows + 1 : )) &
-        ]
-      end associate
-    end associate
-  end procedure
-
-#endif
-
-end submodule mimetic_matrix_s
diff --git a/src/fortran/mole_m.f90 b/src/fortran/mole_m.f90
index e81b61c1..06fbda92 100644
--- a/src/fortran/mole_m.f90
+++ b/src/fortran/mole_m.f90
@@ -1,5 +1,6 @@
 module mole_m 
   !! MOLE Fortran public entities
-  use face_values_m, only : face_values_t, gradient_t
+  use scalar_1D_m, only : scalar_1D_t, scalar_1D_initializer_i
+  use gradient_1D_m, only : gradient_1D_t
   implicit none
 end module mole_m
diff --git a/src/fortran/scalar_1D_m.f90 b/src/fortran/scalar_1D_m.f90
new file mode 100644
index 00000000..dc1e5c25
--- /dev/null
+++ b/src/fortran/scalar_1D_m.f90
@@ -0,0 +1,132 @@
+module scalar_1D_m
+  !! Define an abstraction for the collection of points used to compute gradidents:
+  !! cell centers plus oundaries.
+  use julienne_m, only : file_t
+  use gradient_1D_m, only : gradient_1D_t
+  implicit none
+
+  private
+  public :: scalar_1D_t
+  public :: gradient_operator_1D_t
+  public :: scalar_1D_initializer_i
+
+  abstract interface
+
+    pure function scalar_1D_initializer_i(x) result(f)
+      implicit none
+      double precision, intent(in) :: x(:)
+      double precision, allocatable :: f(:)
+    end function
+
+  end interface
+
+  type mimetic_matrix_1D_t
+    !! Encapsulate a mimetic matrix with a corresponding matrix-vector product operator
+    private
+    double precision, allocatable :: upper_(:,:), inner_(:), lower_(:,:)
+  contains
+    procedure to_file_t
+  end type
+
+  type gradient_operator_1D_t
+    !! Encapsulate kth-order mimetic gradient operator on dx-sized cells
+    private
+    integer k_, m_
+    double precision dx_
+    type(mimetic_matrix_1D_t) mimetic_matrix_1D_
+  end type
+
+  type scalar_1D_t
+    !! Encapsulate information at cell centers and boundaries
+    private
+    double precision, allocatable :: scalar_1D_(:)
+    double precision x_min_, x_max_
+    integer cells_
+    type(gradient_operator_1D_t) gradient_operator_1D_
+  contains
+    procedure grid
+    generic :: operator(.grad.) => grad
+    procedure, non_overridable, private :: grad
+  end type
+
+  interface
+
+     pure module function to_file_t(self) result(file)
+       implicit none
+       class(mimetic_matrix_1D_t), intent(in) :: self
+       type(file_t) file
+     end function
+
+  end interface
+
+  interface scalar_1D_t
+
+    pure module function construct_from_function(initializer, order, cells, x_min, x_max) result(scalar_1D)
+      !! Result is a collection of cell-centered-extended values with a corresponding mimetic gradient operator
+      implicit none
+      procedure(scalar_1D_initializer_i), pointer :: initializer 
+      integer, intent(in) :: order !! order of accuracy
+      integer, intent(in) :: cells !! number of grid cells spanning the domain
+      double precision, intent(in) :: x_min !! grid location minimum
+      double precision, intent(in) :: x_max !! grid location maximum
+      type(scalar_1D_t) scalar_1D
+    end function
+
+  end interface
+
+  interface
+
+    pure module function grid(self) result(x)
+      !! Result is array of cell-centers-extended grid locations (cell centers + boundaries) 
+      !! as described in Corbino & Castillo (2020) https://doi.org/10.1016/j.cam.2019.06.042
+      implicit none
+      class(scalar_1D_t), intent(in) :: self
+      double precision, allocatable :: x(:)
+    end function
+
+    pure module function grad(self) result(grad_f)
+      !! Result is mimetic gradient of f
+      implicit none
+      class(scalar_1D_t), intent(in) :: self
+      type(gradient_1D_t) grad_f !! discrete gradient approximation
+    end function
+
+  end interface
+
+  interface gradient_operator_1D_t
+
+    pure module function construct_from_parameters(k, dx, m) result(gradient_operator_1D)
+      !! Construct a mimetic gradient operator
+      implicit none
+      integer, intent(in) :: k !! order of accuracy
+      double precision, intent(in) :: dx !! step siz
+      integer, intent(in) :: m !! number of grid cells
+      type(gradient_operator_1D_t) gradient_operator_1D
+    end function
+
+  end interface
+
+  interface mimetic_matrix_1D_t
+
+    pure module function construct_from_components(upper, inner, lower) result(mimetic_matrix_1D)
+      !! Construct discrete operator from coefficient matrix
+      implicit none
+      double precision, intent(in) :: upper(:,:), inner(:), lower(:,:)
+      type(mimetic_matrix_1D_t) mimetic_matrix_1D
+    end function
+
+  end interface
+
+  interface
+
+    pure module function matvec(self, scalar_1D) result(matvec_product)
+      !! Apply a mimetic matrix operator to a vector encapsulated in a scalar_1D_t object
+      implicit none
+      class(mimetic_matrix_1D_t), intent(in) :: self
+      type(scalar_1D_t), intent(in) :: scalar_1D
+      double precision, allocatable :: matvec_product(:)
+    end function
+
+  end interface
+
+end module scalar_1D_m
\ No newline at end of file
diff --git a/src/fortran/scalar_1D_s.F90 b/src/fortran/scalar_1D_s.F90
new file mode 100644
index 00000000..c7dde20f
--- /dev/null
+++ b/src/fortran/scalar_1D_s.F90
@@ -0,0 +1,47 @@
+#include "julienne-assert-macros.h"
+
+submodule(scalar_1D_m) scalar_1D_s
+  use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.), operator(.greaterThan.), operator(.isAtLeast.)
+  implicit none
+
+contains
+
+  pure module function construct_from_function(initializer, order, cells, x_min, x_max) result(scalar_1D)
+    implicit none
+    procedure(scalar_1D_initializer_i), pointer :: initializer
+    integer, intent(in) :: order !! order of accuracy
+    integer, intent(in) :: cells !! number of grid cells spanning the domain
+    double precision, intent(in) :: x_min !! grid location minimum
+    double precision, intent(in) :: x_max !! grid location maximum
+    type(scalar_1D_t) scalar_1D
+
+    call_julienne_assert(x_max .greaterThan. x_min)
+    call_julienne_assert(cells .isAtLeast. 2*order)
+
+    scalar_1D%x_min_ = x_min
+    scalar_1D%x_max_ = x_max
+    scalar_1D%cells_ = cells
+    scalar_1D%gradient_operator_1D_ = gradient_operator_1D_t(k=order, dx=(x_max - x_min)/cells, m=cells)
+    scalar_1D%scalar_1D_ = initializer(grid_(x_min, x_max, cells))
+  end function
+
+  pure function grid_(x_min, x_max, cells) result(x)
+    double precision, intent(in) :: x_min, x_max
+    integer, intent(in) :: cells
+    double precision, allocatable :: x(:)
+    integer cell
+
+    associate(dx => (x_max - x_min)/cells)
+      x = [x_min, x_min + dx/2. + [((cell-1)*dx, cell = 1, cells)], x_max]
+    end associate
+  end function
+
+  module procedure grid
+      x = grid_(self%x_min_, self%x_max_, self%cells_)
+  end procedure
+
+  module procedure grad
+    grad_f = gradient_1D_t(matvec(self%gradient_operator_1D_%mimetic_matrix_1D_, self), self%x_min_, self%x_max_, self%cells_)
+  end procedure
+
+end submodule scalar_1D_s
\ No newline at end of file
diff --git a/test/driver.f90 b/test/driver.f90
index 703733f4..724b7fca 100644
--- a/test/driver.f90
+++ b/test/driver.f90
@@ -1,12 +1,10 @@
-! Copyright (c) 2024-2025, The Regents of the University of California and Sourcery Institute
-! Terms of use are as specified in https://github.com/BerkeleyLab/julienne/blob/3.1.2/LICENSE.txt
 program test_suite_driver
   use julienne_m, only : test_fixture_t, test_harness_t
-  use gradient_operator_test_m, only : gradient_operator_test_t
+  use gradient_operator_1D_test_m, only : gradient_operator_1D_test_t
   implicit none
 
   associate(test_harness => test_harness_t([ &
-     test_fixture_t(gradient_operator_test_t()) &
+     test_fixture_t(gradient_operator_1D_test_t()) &
   ]))
     call test_harness%report_results
   end associate
diff --git a/test/gradient_operator_1D_test_m.F90 b/test/gradient_operator_1D_test_m.F90
new file mode 100644
index 00000000..1ef85f7b
--- /dev/null
+++ b/test/gradient_operator_1D_test_m.F90
@@ -0,0 +1,202 @@
+#include "language-support.F90"
+  !! include Julienne preprocessor  macros
+
+module gradient_operator_1D_test_m
+  use julienne_m, only : &
+     string_t &
+    ,test_t &
+    ,test_description_t &
+    ,test_diagnosis_t &
+    ,test_result_t &
+    ,operator(//) &
+    ,operator(.all.) &
+    ,operator(.also.) &
+    ,operator(.approximates.) &
+    ,operator(.csv.) &
+    ,operator(.within.)
+  use mole_m, only : scalar_1D_t, gradient_1D_t, scalar_1D_initializer_i
+#if ! HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
+  use julienne_m, only : diagnosis_function_i
+#endif
+  implicit none
+
+  type, extends(test_t) :: gradient_operator_1D_test_t
+  contains
+    procedure, nopass :: subject
+    procedure, nopass :: results
+  end type
+
+  double precision, parameter :: tight_tolerance = 1D-14, loose_tolerance = 1D-12, rough_tolerance = 1D-02, crude_tolerance = 5D-02
+
+contains
+
+  pure function subject() result(test_subject)
+    character(len=:), allocatable :: test_subject
+    test_subject = 'A 1D mimetic gradient operator'
+  end function
+
+#if HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
+
+  function results() result(test_results)
+    type(gradient_operator_1D_test_t) gradient_operator_1D_test
+    type(test_result_t), allocatable :: test_results(:)
+    test_results = gradient_operator_1D_test%run([ & 
+       test_description_t('computing 2nd- & 4th-order 1D gradients of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const) &
+      ,test_description_t('computing 2nd- & 4th-order 1D gradients of a line within tolerance ' // string_t(loose_tolerance), check_grad_line) &
+      ,test_description_t('computing 2nd- & 4th-order 1D gradients of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola) &
+      ,test_description_t('computing 2nd-order 1D gradients of a sinusoid with a convergence rate of 2 within tolerance ' // string_t(crude_tolerance), check_2nd_order_grad_sinusoid) &
+      ,test_description_t('computing 4th-order 1D gradients of a sinusoid with a convergence rate of 4 within tolerance ' // string_t(crude_tolerance), check_4th_order_grad_sinusoid) &
+    ])
+  end function
+
+#else
+
+  function results() result(test_results)
+    type(gradient_operator_1D_test_t) gradient_operator_1D_test
+    type(test_result_t), allocatable :: test_results(:)
+    procedure(diagnosis_function_i), pointer :: &
+       check_grad_const_ptr => check_grad_const &
+      ,check_grad_line_ptr => check_grad_line &
+      ,check_grad_parabola_ptr => check_grad_parabola
+
+    test_results = gradient_operator_1D_test%run([ &
+       test_description_t('computing 2nd & 4th-order 1D gradients of a constant within tolerance ' // string_t(tight_tolerance), check_grad_const_ptr) &
+      ,test_description_t('computing 2nd & 4th-order 1D gradients of a line within tolerance ' // string_t(loose_tolerance), check_grad_line_ptr) &
+      ,test_description_t('computing 2nd & 4th-order 1D gradients of a parabola within tolerance ' // string_t(loose_tolerance), check_grad_parabola_ptr) &
+    ])
+  end function
+
+#endif
+
+  pure function const(x) result(y)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: y(:)
+    integer i
+    y = [(5D0, i=1,size(x))]
+  end function
+
+  function check_grad_const() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    type(gradient_1D_t) grad_f
+    double precision, parameter :: df_dx = 0.
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => const
+
+    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=4D0)
+    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (2nd-order d(line)/dx)"
+
+    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=8, x_min=0D0, x_max=8D0)
+    test_diagnosis = test_diagnosis .also. (.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance)) // " (4th-order d(line)/dx)"
+  end function
+
+  pure function line(x) result(y)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: y(:)
+    y = 14*x + 3
+  end function
+
+  function check_grad_line() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    type(gradient_1D_t) grad_f
+    double precision, parameter :: df_dx = 14D0
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => line
+
+    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=2, cells=4, x_min=0D0, x_max=4D0)
+    test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (2nd-order d(line)/dx)"
+
+    grad_f = .grad. scalar_1D_t(scalar_1D_initializer, order=4, cells=8, x_min=0D0, x_max=8D0)
+    test_diagnosis = test_diagnosis .also. (.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance)) // " (4th-order d(line)/dx)"
+
+  end function
+
+  pure function parabola(x) result(y)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: y(:)
+    y = 7*x**2 + 3*x + 5
+  end function
+
+  function check_grad_parabola() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    type(scalar_1D_t) quadratic
+    type(gradient_1D_t) grad_f
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => parabola
+
+    quadratic = scalar_1D_t(scalar_1D_initializer , order=2, cells=4, x_min=0D0, x_max=4D0)
+    grad_f = .grad. quadratic
+
+    associate(x => grad_f%faces())
+      associate(df_dx => 14*x + 3)
+        test_diagnosis = .all. (grad_f%values() .approximates. df_dx .within. loose_tolerance) // " (2nd-order d(parabola)/dx)"
+      end associate
+    end associate
+
+    quadratic = scalar_1D_t(scalar_1D_initializer , order=4, cells=8, x_min=0D0, x_max=8D0)
+    grad_f = .grad. quadratic
+
+    associate(x => grad_f%faces())
+      associate(df_dx => 14*x + 3)
+        test_diagnosis = test_diagnosis .also. (.all. (grad_f%values() .approximates. df_dx .within. loose_tolerance)) // " (4th-order d(parabola)/dx)"
+      end associate
+    end associate
+  end function
+
+  pure function sinusoid(x) result(y)
+    double precision, intent(in) :: x(:)
+    double precision, allocatable :: y(:)
+    y = sin(x) + cos(x)
+  end function
+
+  function check_2nd_order_grad_sinusoid() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    type(scalar_1D_t) coarse, fine
+    type(gradient_1D_t) grad_coarse, grad_fine
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => sinusoid
+    double precision, parameter :: pi = 3.141592653589793D0
+    integer, parameter :: order_desired = 2, coarse_cells=100, fine_cells=1000
+
+    coarse = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi)
+    fine   = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0D0, x_max=2*pi)
+
+    grad_coarse = .grad. coarse
+    grad_fine   = .grad. fine
+
+    associate(x_coarse => grad_coarse%faces(), x_fine => grad_fine%faces())
+      associate(df_dx_coarse => cos(x_coarse) - sin(x_coarse), df_dx_fine => cos(x_fine) - sin(x_fine), grad_coarse_values => grad_coarse%values(), grad_fine_values => grad_fine%values())
+        test_diagnosis = .all. (grad_coarse_values .approximates. df_dx_coarse .within. rough_tolerance) // " (2nd-order d(sinusoid)/dx point-wise errors)"
+        test_diagnosis = test_diagnosis .also. (.all. (grad_fine_values .approximates. df_dx_fine .within. rough_tolerance)) // " (2nd-order d(sinusoid)/dx point-wise)"
+        associate(error_coarse_max => maxval(abs(grad_coarse_values - df_dx_coarse)), error_fine_max => maxval(abs(grad_fine_values - df_dx_fine)))
+          associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells))
+            test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. crude_tolerance)  // " (2nd-order d(sinusoid)/dx order of accuracy)"
+          end associate
+        end associate
+      end associate
+    end associate
+  end function
+
+  function check_4th_order_grad_sinusoid() result(test_diagnosis)
+    type(test_diagnosis_t) test_diagnosis
+    type(scalar_1D_t) coarse, fine
+    type(gradient_1D_t) grad_coarse, grad_fine
+    procedure(scalar_1D_initializer_i), pointer :: scalar_1D_initializer => sinusoid
+    double precision, parameter :: pi = 3.141592653589793D0
+    integer, parameter :: order_desired = 4, coarse_cells=100, fine_cells=1000
+
+    coarse = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=coarse_cells, x_min=0D0, x_max=2*pi)
+    fine   = scalar_1D_t(scalar_1D_initializer , order=order_desired, cells=fine_cells  , x_min=0D0, x_max=2*pi)
+
+    grad_coarse = .grad. coarse
+    grad_fine   = .grad. fine
+
+    associate(x_coarse => grad_coarse%faces(), x_fine => grad_fine%faces())
+      associate(df_dx_coarse => cos(x_coarse) - sin(x_coarse), df_dx_fine => cos(x_fine) - sin(x_fine), grad_coarse_values => grad_coarse%values(), grad_fine_values => grad_fine%values())
+        test_diagnosis = .all. (grad_coarse_values .approximates. df_dx_coarse .within. rough_tolerance) // " (4th-order d(sinusoid)/dx point-wise errors)"
+        test_diagnosis = test_diagnosis .also. (.all. (grad_fine_values .approximates. df_dx_fine .within. rough_tolerance)) // " (4th-order d(sinusoid)/dx point-wise)"
+        associate(error_coarse_max => maxval(abs(grad_coarse_values - df_dx_coarse)), error_fine_max => maxval(abs(grad_fine_values - df_dx_fine)))
+          associate(order_actual => log(error_coarse_max/error_fine_max)/log(dble(fine_cells)/coarse_cells))
+            test_diagnosis = test_diagnosis .also. (order_actual .approximates. dble(order_desired) .within. crude_tolerance)  // " (4th-order d(sinusoid)/dx order of accuracy)"
+          end associate
+        end associate
+      end associate
+    end associate
+  end function
+
+end module
\ No newline at end of file
diff --git a/test/gradient_operator_test_m.F90 b/test/gradient_operator_test_m.F90
deleted file mode 100644
index 2df4dd6a..00000000
--- a/test/gradient_operator_test_m.F90
+++ /dev/null
@@ -1,72 +0,0 @@
-! Copyright (c) 2024-2025, The Regents of the University of California and Sourcery Institute
-! Terms of use are as specified in https://github.com/BerkeleyLab/julienne/blob/3.1.2/LICENSE.txt
-
-#include "language-support.F90"
-  !! include Julienne preprocessor  macros
-
-module gradient_operator_test_m
-  use julienne_m, only : &
-     test_t, test_description_t, test_diagnosis_t, test_result_t, operator(.all.), operator(.approximates.), operator(.within.)
-  use mole_m, only : face_values_t
-#if ! HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
-  use julienne_m, only : diagnosis_function_i
-#endif
-  implicit none
-
-  type, extends(test_t) :: gradient_operator_test_t
-  contains
-    procedure, nopass :: subject
-    procedure, nopass :: results
-  end type
-
-contains
-
-  pure function subject() result(test_subject)
-    character(len=:), allocatable :: test_subject
-    test_subject = 'A 1D mimetic gradient operator'
-  end function
-
-#if HAVE_PROCEDURE_ACTUAL_FOR_POINTER_DUMMY
-
-  function results() result(test_results)
-    type(gradient_operator_test_t) gradient_operator_test
-    type(test_result_t), allocatable :: test_results(:)
-    test_results = gradient_operator_test%run([ & 
-       test_description_t('computing the gradient of a linear function', check_gradient_of_line) &
-    ])
-  end function
-
-#else
-
-  function results() result(test_results)
-    type(gradient_operator_test_t) gradient_operator_test
-    type(test_result_t), allocatable :: test_results(:)
-    procedure(diagnosis_function_i), pointer :: &
-      check_gradient_of_line_ptr => check_gradient_of_line
-    test_results = gradient_operator_test%run([ &
-       test_description_t('computing the gradient of a linear function', check_gradient_of_line_ptr) &
-    ])
-  end function
-
-#endif
-
-  function check_gradient_of_line() result(test_diagnosis)
-    type(test_diagnosis_t) test_diagnosis
-    double precision, parameter :: dx = 1D0, df_dx_exact = 1D0
-    double precision, parameter :: x(*) = [0D0,.5D0, 1.5D0, 2.5D0, 3.5D0, 4D0]*dx
-    double precision, parameter :: tolerance = 1D-15
-
-    associate(f => face_values_t(linear(x), k=2, dx=dx))
-      associate(grad_f => .grad. f)
-        test_diagnosis = .all. (grad_f%values() .approximates. df_dx_exact .within. tolerance)
-      end associate
-    end associate
-
-  contains
-    double precision elemental function linear(x)
-      double precision, intent(in) :: x
-      linear = x
-    end function
-  end function
-
-end module