anchor.diff:      !! are extrapolated outside of that subset region
       !!
       !! ---------------------------------
  -    subroutine subset_array(input, output, grid)
  -        implicit none
  -        real,         intent(in)    :: input(:,:)
  -        real,         intent(inout) :: output(:,:)
  +    subroutine extend_array(input, output, grid, extrapolate)
  +        implicit none
  +        real,         intent(in)    :: input(:,:,:)
  +        real,         intent(inout) :: output(:,:,:)
           type(grid_t), intent(in)    :: grid
  +        logical,      intent(in),   optional :: extrapolate
   
           ! loop counter
           integer :: i
  
positive.code:         ! this should only be necessary for border images
          if (grid%ims < 1) then
              do i=1,xs_out-1
                  output(i,:) = output(xs_out,:) + (output(xs_out,:) - output(xs_out+1,:)) * (xs_out - i)
              enddo
          endif
  
          if (grid%ime > nx) then
              do i=xe_out+1,nxo
                  output(i,:) = output(xe_out,:) + (output(xe_out,:) - output(xe_out-1,:)) * (i - xe_out)
              enddo
          endif
  
          if (grid%jms < 1) then
              do i=1,ys_out-1
                  output(:,i) = output(:,ys_out) + (output(:,ys_out) - output(:,ys_out+1)) * (ys_out - i)
              enddo
          endif
  
          if (grid%jme > ny) then
              do i=ye_out+1,nyo
                  output(:,i) = output(:,ye_out) + (output(:,ye_out) - output(:,ye_out-1)) * (i - ye_out)
              enddo
          endif
  
  
negative.code:         call this%v_grid%set_grid_dimensions(       nx_global, ny_global, nz_global, ny_extra = 1)
  
          call this%grid2d%set_grid_dimensions(       nx_global, ny_global, 0)
  
          call this%u_grid2d%set_grid_dimensions(     nx_global, ny_global, 0, nx_extra = 1)
          this%u_grid2d%ims = this%u_grid2d%ims - nsmooth
          this%u_grid2d%ime = this%u_grid2d%ime + nsmooth
  