evalmoment Function

   pure real(rk) function evalmoment(u, grid, order, normalize) result(res)
   !! \(m\)-th order moment of \( u(x,t) \):
   !! $$
   !! M_m(t) = \int _0^\infty x^m u(x,t)dx = {\sum_i \bar{u_i}(t) \Delta x_i x_i^m}
   !! $$
      real(rk), intent(in) :: u(:)
         !! cell-average number density, \( \bar{u_i} \)
      type(grid1), intent(in) :: grid
         !! `grid1` object
      integer, intent(in) :: order
         !! order of the moment
      logical, intent(in), optional :: normalize
         !! if `true`, the result will be normalized by the 0-th moment

      ! Check `order`
      if (order < 0) then
         error stop "Invalid 'order'. Valid range: order >=0."
      end if

      ! Compute moment
      res = sum(u*grid%width*grid%center**order)

      ! Normalize if required
      if (optval(normalize, .false.)) then
         res = res/sum(u*grid%width)
      end if

   end function evalmoment