FIXME : Add documentation
returns results to PE 0 [or to all processors if 'all' is present in input arg list] NOTE: Takes f = f(y, z, sigma, lambda, E, species) and returns int f, where the integral is over all velocity space
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real, | intent(in), | dimension (-ntgrid:,:,p_lo%llim_proc:) | :: | g | ||
| real, | intent(out), | dimension (-ntgrid:,:,:) | :: | total | ||
| integer, | intent(in), | optional | :: | all |
subroutine integrate_moment_r33 (g, total, all)
use mp, only: nproc
use gs2_layouts, only: p_lo, is_idx, ik_idx, ie_idx, il_idx
use mp, only: sum_reduce, sum_allreduce
use theta_grid, only: ntgrid
use array_utils, only: zero_array
implicit none
real, dimension (-ntgrid:,:,p_lo%llim_proc:), intent (in) :: g
real, dimension (-ntgrid:,:,:), intent (out) :: total
integer, optional, intent(in) :: all
integer :: is, il, ie, ik, iplo
!Ensure zero to start
call zero_array(total)
!Do local velocity space integral
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(iplo, is, ik, ie, il) &
!$OMP SHARED(p_lo, w, wl, g) &
!$OMP REDUCTION(+ : total) &
!$OMP SCHEDULE(static)
do iplo = p_lo%llim_proc, p_lo%ulim_proc
ik = ik_idx(p_lo,iplo)
ie = ie_idx(p_lo,iplo)
is = is_idx(p_lo,iplo)
il = il_idx(p_lo,iplo)
total(:, ik, is) = total(:, ik, is) + &
w(ie,is)*wl(:,il)*(g(:,1,iplo)+g(:,2,iplo))
end do
!$OMP END PARALLEL DO
!Do we really need this if?
if (nproc > 1) then
!Complete distributed integral
if (present(all)) then
!Return result to all procs
call sum_allreduce (total)
else
!Only proc0 knows the result
call sum_reduce (total, 0)
end if
end if
end subroutine integrate_moment_r33