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
      integer          :: increment
      integer          :: n,nmin,nmax
      integer          :: neu_shift(3)
      integer          :: ini_shift(3)

      ! Get extents
      extents = this%GetExtents(region)

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

      neu_shift = 0
      ini_shift = 0
      select case (side)
      case (BC_RIGHT)
        increment =  1
      case (BC_LEFT)
        increment = -1
      end select

      if (dir.eq.var%staggering) then

        ! Treatment for face-centered fields in dir-direction
        select case (side)
        case (BC_RIGHT)
          neu_shift(dir) =  0
          ini_shift(dir) =  0
          nmin = 1
          nmax = this%block%ngc
        case (BC_LEFT)
          neu_shift(dir) =  1
          ini_shift(dir) = -increment
          nmin = 0
          nmax = this%block%ngc
        end select

      else

        ! Treatment for cell-centered fields in dir-direction
        neu_shift(dir) =  0
        ini_shift(dir) = -increment
        nmin = 1
        nmax = this%block%ngc

      end if

      ! 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)

              shift = ini_shift

              do n=nmin,nmax
                shift(dir) = shift(dir) + increment
                var%cell(i+shift(1),j+shift(2),k+shift(3)) = var%cell(i+neu_shift(1),j+neu_shift(2),k+neu_shift(3))
              end do

            end do
          end do
        end do
      end associate

      return
    end subroutine bc_set_UpdateBoundaryNeumann