compute_bj_coefficient Function

   real(rk) function compute_bj_coefficient(mesh, phi, grad_phi, c) result(alpha)
      type(mesh_t), intent(in) :: mesh
      real(rk), intent(in) :: phi(:)
      real(rk), intent(in) :: grad_phi(:,:)
      integer, intent(in) :: c

      integer :: lf, fid, nb
      real(rk) :: phi_cell, lo, hi, phi_f, diff, val
      
      phi_cell = phi(c)
      lo = phi_cell
      hi = phi_cell
      
      ! 1. Find the min/max among neighbors
      do lf = 1, mesh%ncell_faces(c)
         fid = mesh%cell_faces(lf, c)
         nb = adjacent_cell(mesh, fid, c)
         if (nb > 0) then
            lo = min(lo, phi(nb))
            hi = max(hi, phi(nb))
         end if
      end do
      
      ! 2. Compute the BJ limiter over all faces of cell c
      alpha = 1.0_rk
      do lf = 1, mesh%ncell_faces(c)
         fid = mesh%cell_faces(lf, c)
         phi_f = phi_cell + dot_product(grad_phi(:, c), mesh%faces(fid)%center - mesh%cells(c)%center)
         diff = phi_f - phi_cell
         if (diff > 1.0e-14_rk) then
            val = min(1.0_rk, (hi - phi_cell) / diff)
            alpha = min(alpha, val)
         else if (diff < -1.0e-14_rk) then
            val = min(1.0_rk, (lo - phi_cell) / diff)
            alpha = min(alpha, val)
         end if
      end do
      
      alpha = max(0.0_rk, min(1.0_rk, alpha))
   end function compute_bj_coefficient