anchor.diff:      end type
   
       interface
  -        subroutine get_grid_at_time(time, x, y, grid)
  +        subroutine get_grid_at_time(time, x, y, grid, grid_update_context)
               !! A subroutine with this interface must be provided to the grid_spacetime_interpolator
               !! It populates the grid given the time, and the x/y coordinates associated with the grid.
  -            import dp
  +            import dp, c_ptr
               implicit none
               real(dp), intent(in) :: time
               real(dp), intent(in) :: x(:), y(:)
               real(dp), intent(inout) :: grid(:,:,:)
  +            type(c_ptr), intent(in) :: grid_update_context
           end subroutine
       end interface
   
       contains
  +
  +    subroutine swap_rank3_arrays(g1, g2)
  +        !! Swap g1 and g2, without any copy.
  +        real(dp), allocatable, intent(inout) :: g1(:,:,:), g2(:,:,:)
  +
  +        real(dp), allocatable :: tmp_g1(:,:,:)
  +
  +        call move_alloc(g1, tmp_g1)
  +        call move_alloc(g2, g1)
  +        call move_alloc(tmp_g1, g2)
  +
  +    end subroutine
   
       subroutine initialise_grid_spacetime_interpolator(gsi, x, y, start_time, dt, ng)
           !! Setup the grid_spacetime_interpolator
  
positive.code:         gsi%g2_time = start_time + dt
  
          ! Read the grids
          call gsi%get_grid_at_time(gsi%g1_time, x, y, gsi%g1)
          call gsi%get_grid_at_time(gsi%g2_time, x, y, gsi%g2)
  
      end subroutine
  
  
negative.code:         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
  
  