pressure_gradient_cell Subroutine

   subroutine pressure_gradient_cell(mesh, bc, p, cell_id, gradp)
      type(mesh_t), intent(in) :: mesh
      type(bc_set_t), intent(in) :: bc
      real(rk), intent(in) :: p(:)
      integer, intent(in) :: cell_id
      real(rk), intent(out) :: gradp(3)

      integer :: lf, f, nb, patch_id
      real(rk) :: nvec(3)
      real(rk) :: pf, dist
      logical :: is_dirichlet

      gradp = zero

      do lf = 1, mesh%ncell_faces(cell_id)
         f = mesh%cell_faces(lf, cell_id)
         nvec = outward_normal(mesh, f, cell_id)

         nb = face_effective_neighbor(mesh, bc, f, cell_id)

         if (nb > 0) then
            pf = face_linear_scalar(mesh, bc, f, cell_id, nb, p(cell_id), p(nb))
         else
            if (boundary_pressure_type(mesh, bc, f) == bc_dirichlet) then
               call boundary_pressure(mesh, bc, f, p(cell_id), pf, is_dirichlet)
            else
               patch_id = mesh%faces(f)%patch
               if (patch_id > 0) then
                  dist = face_normal_distance(mesh, bc, f, cell_id, 0)
                  pf = p(cell_id) + bc%patches(patch_id)%dpdn * dist
               else
                  pf = p(cell_id)
               end if
            end if
         end if

         gradp = gradp + pf * nvec * mesh%faces(f)%area
      end do

      gradp = gradp / mesh%cells(cell_id)%volume
   end subroutine pressure_gradient_cell