cdifs_obj_PrepareSolverOperatorsDIV Subroutine

    subroutine cdifs_obj_PrepareSolverOperatorsDIV(this)
      !> Generates the divergence operator for the velocity
      ! 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
      logical          :: found
      type(extent_obj) :: extents
      integer          :: dir
      character(len=2) :: varname(3)
      integer          :: ierr

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

      ! Define divergence operator
      associate (lo => this%block%lo, hi => this%block%hi, &
        op => this%op, bcs => this%bcs)
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)

              do n=-op%st+1,op%st
                this%DIV(n,i,j,k,1) = this%DIV(n,i,j,k,1) + op%c_d1dx1 (n,i)
                this%DIV(n,i,j,k,2) = this%DIV(n,i,j,k,2) + op%c_d1dx2 (n,j)
                this%DIV(n,i,j,k,3) = this%DIV(n,i,j,k,3) + op%c_d1dx3 (n,k)
              end do

            end do
          end do
        end do

        ! Apply BC
        varname=['V1','V2','V3']
        do n=1,bcs%count
          do dir=1,3
            ! Velocity component in the this direction
            found = bcs%CheckBCExists(bcs%region(n)%name, varname(dir))
            if (found) then
              select case ( bcs%GetBCType(bcs%region(n)%name, varname(dir)) )
              case (BC_NEUMANN)
                ! Get extents
                extents = bcs%GetExtents(bcs%region(n)%name)

                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%DIV(:,i,j,k,dir) = 0.0_wp
                    end do
                  end do
                end do

              end select
            end if
          end do
        end do
      end associate

      return
    end subroutine cdifs_obj_PrepareSolverOperatorsDIV