Calculate the flux surface average term for the adiabatic response.
If using adiabatic electrons and the option
dist_fn_knobs
adiabatic_option = "iphi00=2"
/
the field solve should subtract the flux surface average
from the perturbed electron distribution function,
i.e. include the term
f1e = q( phi -
This function calculates
Joseph Parker, STFC
We might want to rename this routine to reflect that it only calculates a very specific flux surface average and not a general average of any field shaped array (e.g. as gamtot3 is only allocated/set in a special case and the gamtot factor is not required in the general flux surface average).
antot/gamtot is the expression for phi when beta = 0 and this reflects what is averaged in this routine. As such we appear to be assuming beta = 0 here. This is probably fine as we should only use the result for adiabatic electrons (where one would expect to be ignoring EM perturbations) but perhaps there should be a consistency check? This is something we should still flag and could rectify with a little effort.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| complex, | intent(out), | dimension (ntheta0, naky) | :: | fl_avg | ||
| complex, | intent(in), | dimension (-ntgrid:ntgrid,ntheta0,naky) | :: | antot |
subroutine calculate_flux_surface_average (fl_avg,antot)
use kt_grids, only: naky, ntheta0, aky, kwork_filter
use run_parameters, only: tite
use theta_grid, only: ntgrid, theta, jacob
use integration, only: trapezoidal_integration
implicit none
complex, dimension (ntheta0, naky), intent (out) :: fl_avg
complex, dimension (-ntgrid:ntgrid,ntheta0,naky), intent (in) :: antot
integer :: ik, it
fl_avg = 0.
if (.not. allocated(awgt)) then
allocate (awgt(ntheta0, naky))
awgt = 0.
do ik = 1, naky
if (aky(ik) > epsilon(0.0)) cycle
do it = 1, ntheta0
! We might not want this cycle given that this is a one-off setup of awgt
! and kwork_filter may change during the run (e.g. between initialisation
! and time advance). I suspect kwork_filter will be false the first
! time we call this routine, but there is no promise.
if(kwork_filter(it,ik)) cycle
awgt(it,ik) = 1.0/trapezoidal_integration(theta, jacob*gamtot3(:,it,ik))
end do
end do
endif
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(ik, it) &
!$OMP SHARED(naky, ntheta0, kwork_filter, fl_avg, tite, antot, &
!$OMP theta, jacob, gamtot, awgt) &
!$OMP SCHEDULE(static)
do ik = 1, naky
! Note awg == 0 for aky > 0 so we might be able to cycle/limit the loop.
do it = 1, ntheta0
if(kwork_filter(it,ik)) cycle
fl_avg(it,ik) = tite*trapezoidal_integration(theta, jacob*antot(:,it,ik)/gamtot(:,it,ik))*awgt(it,ik)
end do
end do
!$OMP END PARALLEL DO
end subroutine calculate_flux_surface_average