fluxavg Function

  pure real function fluxavg(f, theta, jacob)
    implicit none
    real, dimension(-nth:nth), intent(in) :: f, theta, jacob
    real, dimension(-nth:nth) :: delth
    delth(-nth+1:nth-1) = 0.5 * (theta(-nth+2:nth) - theta(-nth:nth-2))
    delth(-nth) = 0.5 * (theta(-nth+1) - theta(-nth))
    delth(nth) = 0.5 * (theta(nth) - theta(nth-1))

    ! We might want to write this to use trapezoidal_integration instead
    ! of this sum form. Currently only used in [[f_trap]], which itself
    ! is only used here for one piece of information written to screen.
    ! As it doesn't get into the simulated system we probably don't mind
    ! too much beyond readability concerns.
    fluxavg = sum(f * jacob * delth) / sum(jacob * delth)
  end function fluxavg