Calculate the phi0
array, used to normalise eigenfunctions and
moments.
function get_phi0() result(phi0)
use fields_arrays, only: phi, apar
use kt_grids, only: ntheta0, naky
use dist_fn, only: def_parity, even
use run_parameters, only: fapar
use theta_grid, only: theta
implicit none
complex, dimension (ntheta0, naky) :: phi0
!Should this actually use igomega instead of 0?
!What if fphi==0? --> See where statements below
phi0 = phi(0,:,:)
!This looks like a hack for the case where we know we've forced phi(theta=0) to be 0
!this could probably be better addressed by the use of igomega above
if (def_parity .and. fapar > 0 .and. (.not. even)) phi0 = apar(0, :, :)
!Address locations where phi0=0 by using next point
where (abs(phi0) < 10.0*epsilon(0.0))
phi0 = phi(1,:,:)/(theta(1)-theta(0))
end where
!Check again if any locations are 0, this could be true if fphi (or fapar)
!is zero.
where (abs(phi0) < 10.0*epsilon(0.0))
phi0 = 1.0
end where
end function get_phi0