bc_set_UpdateBoundaryNeumann Subroutine

    impure subroutine bc_set_UpdateBoundaryNeumann(this,region,var)
      !> Updates ghostcells to enforce Neumann BC.
      class(bc_set),        intent(in) :: this                                 !! Boundary conditions utility
      character(len=*),     intent(in)    :: region                            !! Region name
      type(eulerian_obj_r), intent(inout) :: var                               !! Eulerian variable
      ! Work variables
      type(extent_obj) :: extents
      integer          :: i,j,k
      integer          :: shift(3)
      integer          :: dir
      integer          :: side

      ! Get extents
      extents = this%GetExtents(region)

      ! Get direction and side of BC
      call this%GetSideDirByRegion(region,side,dir)

      shift = 0

      if (dir.eq.var%staggering) then
        ! Treatment for face-centered fields in dir-direction
        select case (side)
        case (BC_RIGHT)
          shift(dir) = 1
          ! Assuming zero gradient
          associate(lo => extents%lo, hi=>extents%hi)
            do k=lo(3),hi(3)
              do j=lo(2),hi(2)
                do i=lo(1),hi(1)
                  var%cell(i+shift(1),j+shift(2),k+shift(3)) = var%cell(i,j,k)
                end do
              end do
            end do
          end associate
        case (BC_LEFT)
          shift(dir) = 1
          ! Assuming zero gradient
          associate(lo => extents%lo, hi=>extents%hi)
            do k=lo(3),hi(3)
              do j=lo(2),hi(2)
                do i=lo(1),hi(1)
                  var%cell(i,j,k) = var%cell(i+shift(1),j+shift(2),k+shift(3))
                end do
              end do
            end do
          end associate
        end select


      else

        ! Treatment for cell-centered fields in dir-direction
        select case (side)
        case (BC_RIGHT)
          shift(dir) =  1
        case (BC_LEFT)
          shift(dir) = -1
        end select

        ! Assuming zero gradient
        associate(lo => extents%lo, hi=>extents%hi)
          do k=lo(3),hi(3)
            do j=lo(2),hi(2)
              do i=lo(1),hi(1)
                var%cell(i+shift(1),j+shift(2),k+shift(3)) = var%cell(i,j,k)
              end do
            end do
          end do
        end associate

      end if

      return
    end subroutine bc_set_UpdateBoundaryNeumann