anchor.diff:              ng_local = 1_ip
           end if
   
  +        ! Check if x is increasing and evenly spaced
  +        gsi%even_x_spacing = .true.
  +        dx = (x(size(x)) - x(1) ) / (size(x) - 1.0_dp)
           do i = 1, (size(x, kind=ip)-1)
               if(x(i+1) <= x(i)) then
                   write(log_output_unit, *) "x must be increasing", x
                   call generic_stop
               end if
  -        end do
  +            if( abs(x(i+1) - x(i) - dx) > 10*maxval(spacing([dx, x(i:i+1)])) ) gsi%even_x_spacing = .false.
  +        end do
  +
  +        ! Check if y is increasing and evenly spaced
  +        gsi%even_y_spacing = .true.
  +        dy = (y(size(y)) - y(1) ) / (size(y) - 1.0_dp)
           do i = 1, (size(y, kind=ip)-1)
               if(y(i+1) <= y(i)) then
                   write(log_output_unit, *) "y must be increasing", y
                   call generic_stop
               end if
  +            if( abs(y(i+1) - y(i) - dy) > 10*maxval(spacing([dy, y(i:i+1)])) ) gsi%even_y_spacing = .false.
           end do
   
           ! Grid size
  
positive.code:         do i = 1, size(xs, kind=ip)
  
              ! Get x-indices for bilinear interpolation: i0 < i1
              call nearest_index_sorted(gsi%nx(1), gsi%x, xs(i), i0)
              if(xs(i) > gsi%x(i0)) then
                  i1 = min(i0+1, gsi%nx(1))
              else
  
negative.code:         end do
          !print*, ys
  
          do testcase = 1, 2
  
              ! Setup the gsi object
              start_time = 35.0_dp
              dt = 15.0_dp
              gsi%grid_update_context = c_loc(cxy(1))
              gsi%get_grid_at_time => for_test_get_grid_at_time
              call gsi%initialise(xs, ys, start_time, dt, ng)
  
              if(testcase == 1) then
                  ! The x/y spacing is even, which should be detected
                  if(gsi%even_x_spacing .and. gsi%even_y_spacing) then
                      print*, 'PASS'
                  else
                      print*, 'FAIL'
                  end if
              else
                  ! Forcibly set even_x_spacing/even_y_spacing to .false., so that
                  ! we exercise the test code for uneven spacing. 
                  ! This is not recommended in general usage!
                  gsi%even_x_spacing = .false.
                  gsi%even_y_spacing = .false.
              end if
  
              !
              ! Setup the test data -- off the grid in this case
              !
              tx = x0 - 1.0_dp
              ty = y0 - 1.0_dp
              time = start_time + 5.0_dp
              call gsi%interpolate(time, tx, ty, tv)
              do k = 1, ng
                  if(all(abs(tv(:,k) - k*(time + cxy(1) * x0 + cxy(2) * y0)) < errtol)) then
                      print*, 'PASS'
                  else
                      print*, 'FAIL', abs(tv(:,k) - k*(time + cxy(1) * x0 + cxy(2) * y0))
                  end if
              end do
  
              !
              ! Setup the test data -- on the grid in this case
              !
              tx = x0 + dx * ([1.0 , 3.0, 6.3, 4.0, 8.0] - 1.0)
              ty = y0 + dy * ([20.0, 4.0, 5.2, 1.0, 2.0] - 1.0)
              time = start_time + 5.0_dp
              call gsi%interpolate(time, tx, ty, tv)
              do k = 1, ng
                  if(all(abs(tv(:,k) - k*(time + cxy(1) * tx + cxy(2) * ty)) < errtol)) then
                      print*, 'PASS'
                  else
                      print*, 'FAIL', abs(tv(:,k) - k*(time + cxy(1) * tx + cxy(2) * ty))
                  end if
              end do
  
              ! As before, with time large enough to trigger a grid update
              time = start_time + 5.0_dp + dt
              call gsi%interpolate(time, tx, ty, tv) ! This will cause g1 => g2, and g2 to be updated using for_test_get_grid_at_time
              do k = 1, ng
                  if(all(abs(tv(:,k) - k*(time + cxy(1) * tx + cxy(2) * ty)) < errtol)) then
                      print*, 'PASS'
                  else
                      print*, 'FAIL', abs(tv(:,k) - k*(time + cxy(1) * tx + cxy(2) * ty))
                  end if
              end do
  
              ! Finally check we can interpolate the grids in time
              time = 0.6_dp*gsi%g1_time + 0.4_dp*gsi%g2_time
              if(abs(time - (0.6_dp*50.0_dp + 0.4_dp*65.0_dp)) < errtol) then
                  print*, 'PASS'
              else
                  print*, 'FAIL', time, time - (0.6_dp*50.0_dp + 0.4_dp*65.0_dp), errtol
              end if
              call gsi%interpolate_grid_in_time(time, test_grid)
              if(all(abs(test_grid - (0.6_dp * gsi%g1 + 0.4_dp * gsi%g2)) < errtol)) then
                  print*, 'PASS'
              else
                  print*, 'FAIL'
              end if
  
              call gsi%finalise
  
          end do
  
          ! Check that uneven x/y coordinates are detected
          xs(1) = xs(1) - 0.1_dp
          ys(size(ys)) = ys(size(ys)) + 0.4_dp
          gsi%grid_update_context = c_loc(cxy(1))
          gsi%get_grid_at_time => for_test_get_grid_at_time
          call gsi%initialise(xs, ys, start_time, dt, ng)
          ! The x/y spacing is uneven, which should be detected
          if(gsi%even_x_spacing .or. gsi%even_y_spacing) then
              print*, 'FAIL'
          else
              print*, 'PASS'
          end if
  
      end subroutine
  
  