calculate_simple_quasilinear_flux_metric_by_k Function

  function calculate_simple_quasilinear_flux_metric_by_k(growth_rates) result(ql_metric)
    use run_parameters, only: has_phi, has_apar, has_bpar
    use kt_grids, only: naky, ntheta0, aky
    use fields_arrays, only: phinew, aparnew, bparnew
    use warning_helpers, only: is_zero
    implicit none
    real, dimension(ntheta0, naky), intent(in) :: growth_rates
    real, dimension(ntheta0, naky) :: limited_growth_rates, average
    real, dimension(ntheta0, naky) :: ql_metric, normalising_field
    integer :: ik

    ! Initialise flux to zero
    ql_metric = 0.0

    ! Decide on the normalising field
    if (has_phi) then
       normalising_field = maxval(abs(phinew), dim = 1)
    else if (has_apar) then
       normalising_field = maxval(abs(aparnew), dim = 1)
    else if (has_bpar) then
       normalising_field = maxval(abs(bparnew), dim = 1)
    else
       ! If we get here then no fields are active so the QL flux
       ! is zero and we can exit
       return
    end if

    where (normalising_field > 0)
       normalising_field = 1.0 / normalising_field
    end where

    if (has_phi) then
       average = get_field_weighted_average_kperp2(phinew)
       where(average > 0.0)
          ql_metric = maxval(abs(phinew), dim = 1) * normalising_field / average
       end where
    end if

    if (has_apar) then
       average = get_field_weighted_average_kperp2(aparnew)
       where(average > 0.0)
          ql_metric = ql_metric + &
               maxval(abs(aparnew), dim = 1) * normalising_field / average
       end where
    end if

    if (has_bpar) then
       average = get_field_weighted_average_kperp2(bparnew)
       where(average > 0.0)
          ql_metric = ql_metric + &
               maxval(abs(bparnew), dim = 1) * normalising_field / average
       end where
    end if

    ! Limit the growth rate to positive values
    where (growth_rates > 0)
       limited_growth_rates = growth_rates
    elsewhere
       limited_growth_rates = 0.0
    end where

    ! MG: Avoid spurious contribution from the zonal mode
    do ik = 1, naky
       if (is_zero(aky(ik))) ql_metric(:, ik) = 0.0
    end do

    ql_metric = ql_metric * limited_growth_rates
  end function calculate_simple_quasilinear_flux_metric_by_k