op_obj_div Function

    impure function op_obj_div(this,name,in1,in2,in3) result(out)
      !> Computes the divergence of a vector (in1,in2,in3).
      ! This function takes in1,in2,in3 cell-centered (stag=0)
      ! and returns the divergence on cell centers (stag=0)
      class(op_obj),        intent(in) :: this                                 !! Differential operators utility
      character(len=*),     intent(in) :: name                                 !! Name to give this variable
      type(eulerian_obj_r), intent(in) :: in1                                  !! Component in 1-dir
      type(eulerian_obj_r), intent(in) :: in2                                  !! Component in 2-dir
      type(eulerian_obj_r), intent(in) :: in3                                  !! Component in 3-dir
      type(eulerian_obj_r)             :: out                                  !! Result
      ! Work variables
      integer :: i,j,k,n,m

      call out%Initialize(name,this%block,this%parallel, 0 )

      associate (lo => this%block%lo, hi => this%block%hi )
        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
                  out%cell(i,j,k) = out%cell(i,j,k)                               &
                    + this%c_d1dx1(n,i)*this%c_intrp1m(m,i+n)*in1%cell(i+n+m,j,k) &
                    + this%c_d1dx2(n,j)*this%c_intrp2m(m,j+n)*in2%cell(i,j+n+m,k) &
                    + this%c_d1dx3(n,k)*this%c_intrp3m(m,k+n)*in3%cell(i,j,k+n+m)
                end do
              end do
            end do
          end do
        end do
      end associate

      return
    end function op_obj_div