cdifs_obj_PrepareSolverOperatorsPGRAD Subroutine

    subroutine cdifs_obj_PrepareSolverOperatorsPGRAD(this)
      !> Generates the gradient operator for the pressure
      ! field, and adjusts at the boundaries.
      use leapBC
      implicit none
      class(cdifs_obj), intent(inout) :: this                                  !! CDIFS solver
      ! Work variables
      integer          :: i,j,k,n,m
      logical          :: found
      type(extent_obj) :: extents
      integer          :: shift(3)
      integer          :: dir
      integer          :: side
      integer          :: ierr

      ! Allocate array
      associate (lo => this%block%lo, hi => this%block%hi, &
        op => this%op, bcs => this%bcs)
        allocate(this%pGRAD(-op%st:op%st-1,lo(1)-op%st+1:hi(1)+op%st,lo(2)-op%st+1:hi(2)+op%st,lo(3)-op%st+1:hi(3)+op%st,3), &
                 source = 0.0_wp, stat=ierr)
        if (ierr.ne.0) call this%parallel%Stop('Unable to allocate pGRAD array')
      end associate

      ! Define gradients of pressure
      associate (lo => this%block%lo, hi => this%block%hi, &
        op => this%op, bcs => this%bcs)
        do k=lo(3)-op%st+1,hi(3)+op%st
          do j=lo(2)-op%st+1,hi(2)+op%st
            do i=lo(1)-op%st+1,hi(1)+op%st

              do m=-op%st,op%st-1
                this%pGRAD(m,i,j,k,1) = this%pGRAD(m,i,j,k,1) + op%c_d1dx1m(m,i)
                this%pGRAD(m,i,j,k,2) = this%pGRAD(m,i,j,k,2) + op%c_d1dx2m(m,j)
                this%pGRAD(m,i,j,k,3) = this%pGRAD(m,i,j,k,3) + op%c_d1dx3m(m,k)
              end do

            end do
          end do
        end do

        ! Apply Neumann BC
        do n=1,bcs%count
          ! Check Whether we have a BC for this variable on this region
          found = bcs%CheckBCExists(bcs%region(n)%name,'P')
          if (.not.found) cycle

          select case (bcs%GetBCType(bcs%region(n)%name,'P'))
          case (BC_NEUMANN)
           ! Get extents
           extents = bcs%GetExtents(bcs%region(n)%name)

           ! Get direction and side of BC
           call bcs%GetSideDirByRegion(bcs%region(n)%name,side,dir)

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

           ! Assuming zero gradient
           do k=extents%lo(3),extents%hi(3)
             do j=extents%lo(2),extents%hi(2)
               do i=extents%lo(1),extents%hi(1)
                 this%pGRAD(:,i+shift(1),j+shift(2),k+shift(3),dir) = 0.0_wp
               end do
             end do
           end do
          end select
        end do
      end associate

      return
    end subroutine cdifs_obj_PrepareSolverOperatorsPGRAD