op_obj_BuildLaplacian Subroutine

    pure subroutine op_obj_BuildLaplacian(this,mat,stm)
      !> Builds Laplacian operator using Morinishi's schemes.
      class(op_obj), intent(inout) :: this                                     !! Differential operators utility
      real(wp),      intent(out),   &
                       allocatable :: mat(:,:,:,:,:)                           !! Matrix that stores the laplacian operator
      integer,       intent(out)   :: stm                                      !! Stencil extent
      ! Work variables
      integer :: i,j,k,n,m

      if (allocated(mat)) deallocate(mat)

      ! Set stencil extent based on current scheme order
      stm = 2*this%st-1

      associate ( lo=>this%block%lo, hi =>this%block%hi  )

        allocate(mat(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),3,-stm:stm))

        mat = 0.0_wp
        do k=lo(3),hi(3)
          do j=lo(2),hi(2)
            do i=lo(1),hi(1)

              do n=-this%st+1,this%st
                do m=-this%st,this%st-1
                  mat(i,j,k,1,n+m)= mat(i,j,k,1,n+m) + this%c_d1dx1 (n,i)*this%c_d1dx1m(m,i+n)
                  mat(i,j,k,2,n+m)= mat(i,j,k,2,n+m) + this%c_d1dx2 (n,j)*this%c_d1dx2m(m,j+n)
                  mat(i,j,k,3,n+m)= mat(i,j,k,3,n+m) + this%c_d1dx3 (n,k)*this%c_d1dx3m(m,k+n)
                end do
              end do

            end do
          end do
        end do

      end associate

      return
    end subroutine op_obj_BuildLaplacian