cdifs_obj_PrepareSolverOperatorsVLAP Subroutine

    subroutine cdifs_obj_PrepareSolverOperatorsVLAP(this)
      !> Generates the viscous Laplacian operator 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
      integer          :: stm
      logical          :: found
      type(extent_obj) :: extents
      integer          :: dir
      integer          :: bcdir
      integer          :: side
      integer          :: varstag(3)
      character(len=2) :: varname(3)
      integer          :: ierr

      associate (lo => this%block%lo, hi => this%block%hi, &
        op => this%op, bcs => this%bcs)
        ! Set stencil extent based on current scheme order
        stm = 2*op%st-1

        allocate(this%LAP(-stm:stm,lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),3,3), &
                 source = 0.0_wp, stat=ierr)
        if (ierr.ne.0) call this%parallel%Stop('Unable to allocate LAP array')
      end associate

      ! Define the viscous Laplacian
      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
                do m=-op%st,op%st-1

                  ! Laplacian for velocity component in x1-dir
                  this%LAP(n+m,i,j,k,1,1) = this%LAP(n+m,i,j,k,1,1) + op%c_d1dx1m(m,i)*op%c_d1dx1 (n,i+m)
                  this%LAP(n+m,i,j,k,2,1) = this%LAP(n+m,i,j,k,2,1) + op%c_d1dx2 (n,j)*op%c_d1dx2m(m,j+n)
                  this%LAP(n+m,i,j,k,3,1) = this%LAP(n+m,i,j,k,3,1) + op%c_d1dx3 (n,k)*op%c_d1dx3m(m,k+n)

                  ! Laplacian for velocity component in x2-dir
                  this%LAP(n+m,i,j,k,1,2) = this%LAP(n+m,i,j,k,1,2) + op%c_d1dx1 (n,i)*op%c_d1dx1m(m,i+n)
                  this%LAP(n+m,i,j,k,2,2) = this%LAP(n+m,i,j,k,2,2) + op%c_d1dx2m(m,j)*op%c_d1dx2 (n,j+m)
                  this%LAP(n+m,i,j,k,3,2) = this%LAP(n+m,i,j,k,3,2) + op%c_d1dx3 (n,k)*op%c_d1dx3m(m,k+n)

                  ! Laplacian for velocity component in x3-dir
                  this%LAP(n+m,i,j,k,1,3) = this%LAP(n+m,i,j,k,1,3) + op%c_d1dx1 (n,i)*op%c_d1dx1m(m,i+n)
                  this%LAP(n+m,i,j,k,2,3) = this%LAP(n+m,i,j,k,2,3) + op%c_d1dx2 (n,j)*op%c_d1dx2m(m,j+n)
                  this%LAP(n+m,i,j,k,3,3) = this%LAP(n+m,i,j,k,3,3) + op%c_d1dx3m(m,k)*op%c_d1dx3 (n,k+m)
                end do
              end do

            end do
          end do
        end do

        varname=['V1','V2','V3']
        varstag=[  1,   2,   3]
        ! Apply BC
        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_DIRICHLET, BC_NEUMANN)
                ! Get direction and side of BC
                call bcs%GetSideDirByRegion(bcs%region(n)%name,side,bcdir)

                ! If dir.eq.stagggering and side.eq.BC_LEFT
                if ( bcdir.eq.varstag(dir) .and. side.eq.0 ) then
                  ! 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%LAP(:,i,j,k,1,dir) = 0.0_wp
                        this%LAP(:,i,j,k,2,dir) = 0.0_wp
                        this%LAP(:,i,j,k,3,dir) = 0.0_wp
                      end do
                    end do
                  end do
                end if
              end select
            end if
          end do
        end do

      end associate

      return
    end subroutine cdifs_obj_PrepareSolverOperatorsVLAP