anchor.diff:          ! 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)
  +                if (do_extrapolate) then
  +                    output(i,:) = output(xs_out,:) + (output(xs_out,:) - output(xs_out+1,:)) * (xs_out - i)
  +                else
  +                    output(i,:) = output(xs_out,:)
  +                endif
               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)
  +                if (do_extrapolate) then
  +                    output(i,:) = output(xe_out,:) + (output(xe_out,:) - output(xe_out-1,:)) * (i - xe_out)
  +                else
  +                    output(i,:) = output(xe_out,:)
  +                endif
               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)
  +                if (do_extrapolate) then
  +                    output(:,i) = output(:,ys_out) + (output(:,ys_out) - output(:,ys_out+1)) * (ys_out - i)
  +                else
  +                    output(:,i) = output(:,ys_out)
  +                endif
               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)
  +                if (do_extrapolate) then
  +                    output(:,i) = output(:,ye_out) + (output(:,ye_out) - output(:,ye_out-1)) * (i - ye_out)
  +                else
  +                    output(:,i) = output(:,ye_out)
  +                endif
               enddo
           endif
   
  
positive.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
  
negative.code:     !! are extrapolated outside of that subset region
      !!
      !! ---------------------------------
      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
  