box_gc_for_fine_neighbor Subroutine

  subroutine box_gc_for_fine_neighbor(box, nb, di, nc, iv, gc)
    type(mg_box_t), intent(in) :: box
    integer, intent(in)     :: nb
    !< Direction of fine neighbor
    integer, intent(in)     :: di(3)
    !< Index offset of fine neighbor
    integer, intent(in)     :: nc, iv
    real(dp), intent(out)   :: gc(nc, nc)
    real(dp)                :: tmp(0:nc/2+1, 0:nc/2+1), grad(3-1)
    integer                 :: i, j, hnc

    hnc = nc/2

    ! First fill a temporary array with data next to the fine grid
    select case (nb)
    case (mg_neighb_lowx)
       tmp = box%cc(1, di(2):di(2)+hnc+1, di(3):di(3)+hnc+1, iv)
    case (mg_neighb_highx)
       tmp = box%cc(nc, di(2):di(2)+hnc+1, di(3):di(3)+hnc+1, iv)
    case (mg_neighb_lowy)
       tmp = box%cc(di(1):di(1)+hnc+1, 1, di(3):di(3)+hnc+1, iv)
    case (mg_neighb_highy)
       tmp = box%cc(di(1):di(1)+hnc+1, nc, di(3):di(3)+hnc+1, iv)
    case (mg_neighb_lowz)
       tmp = box%cc(di(1):di(1)+hnc+1, di(2):di(2)+hnc+1, 1, iv)
    case (mg_neighb_highz)
       tmp = box%cc(di(1):di(1)+hnc+1, di(2):di(2)+hnc+1, nc, iv)
    case default
       error stop
    end select

    ! Now interpolate the coarse grid data to obtain values 'straight' next to
    ! the fine grid points
    do j = 1, hnc
       do i = 1, hnc
          grad(1)          = 0.125_dp * (tmp(i+1, j) - tmp(i-1, j))
          grad(2)          = 0.125_dp * (tmp(i, j+1) - tmp(i, j-1))
          gc(2*i-1, 2*j-1) = tmp(i, j) - grad(1) - grad(2)
          gc(2*i, 2*j-1)   = tmp(i, j) + grad(1) - grad(2)
          gc(2*i-1, 2*j)   = tmp(i, j) - grad(1) + grad(2)
          gc(2*i, 2*j)     = tmp(i, j) + grad(1) + grad(2)
       end do
    end do
  end subroutine box_gc_for_fine_neighbor