diffusivity Function

  pure real function diffusivity()
    use kt_grids, only: ntheta0, naky, kperp2
    use theta_grid, only: grho
    use warning_helpers, only: is_zero
    real, dimension(ntheta0, naky) :: diffusivity_by_k
    integer  :: ik, it
    diffusivity_by_k = 0.0
    do ik=1,naky
       do it=1,ntheta0
          if (is_zero(kperp2(igomega,it,ik))) cycle
          diffusivity_by_k(it,ik) = &
               max(aimag(omegaavg(it,ik)), 0.0)/kperp2(igomega, it, ik)*2.0
       end do
    end do
    diffusivity = maxval(diffusivity_by_k) * grho(igomega)
  end function diffusivity