FIXME : Add documentation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! First get v0y0 $OMP PARALLEL DO DEFAULT(none) & $OMP PRIVATE(iglo, it, ik, ie, is, isgn) & $OMP SHARED(g_lo, kwork_filter, gtmp, vnmult, vnew_E, aj0, g) & $OMP SCHEDULE(static) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Get y1 = y0 - v0y0 * z0 / (1 + v0z0) $OMP PARALLEL DO DEFAULT(none) & $OMP PRIVATE(iglo, it, ik, is, isgn) & $OMP SHARED(g_lo, kwork_filter, g1, g, ieqzip, v0y0, bz0) & $OMP SCHEDULE(static) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Now get v1y1 $OMP PARALLEL DO DEFAULT(none) & $OMP PRIVATE(iglo, it, ik, ie, is, isgn) & $OMP SHARED(g_lo, kwork_filter, gtmp, vns, vpa, aj0, g1) & $OMP SCHEDULE(static) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Get y2 = y1 - v1y1 * s1 / (1 + v1s1) $OMP PARALLEL DO DEFAULT(none) & $OMP PRIVATE(iglo, it, ik, is, isgn) & $OMP SHARED(g_lo, kwork_filter, g1, ieqzip, v1y1, bs0) & $OMP SCHEDULE(static) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Now get v2y2 $OMP PARALLEL DO DEFAULT(none) & $OMP PRIVATE(iglo, it, ik, ie, il, is, isgn) & $OMP SHARED(g_lo, kwork_filter, gtmp, vns, energy, al, aj1, g1) & $OMP SCHEDULE(static) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Finally get x = y2 - v2y2 * w2 / (1 + v2w2) $OMP PARALLEL DO DEFAULT(none) & $OMP PRIVATE(iglo, it, ik, is, isgn) & $OMP SHARED(g_lo, kwork_filter, g, g1, ieqzip, v2y2, bw0) & $OMP SCHEDULE(static)
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
complex, | intent(inout), | dimension (-ntgrid:,:,g_lo%llim_proc:) | :: | g | ||
complex, | intent(inout), | dimension (-ntgrid:,:,g_lo%llim_proc:) | :: | g1 |
subroutine conserve_diffuse_standard_layout (g, g1)
use theta_grid, only: ntgrid
use species, only: nspec
use kt_grids, only: naky, ntheta0, kwork_filter
use gs2_layouts, only: g_lo, ik_idx, it_idx, ie_idx, il_idx, is_idx
use le_grids, only: energy => energy_maxw, al, integrate_moment, negrid
use dist_fn_arrays, only: aj0, aj1, vpa
use run_parameters, only: ieqzip
use array_utils, only: zero_array
implicit none
complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in out) :: g, g1
complex, dimension (:,:,:), allocatable :: gtmp
real, dimension (:,:,:), allocatable :: vns
integer :: isgn, iglo, ik, ie, il, is, it
complex, dimension (:,:,:,:), allocatable :: v0y0, v1y1, v2y2
logical, parameter :: all_procs = .true.
allocate (v0y0(-ntgrid:ntgrid, ntheta0, naky, nspec))
allocate (v1y1(-ntgrid:ntgrid, ntheta0, naky, nspec))
allocate (v2y2(-ntgrid:ntgrid, ntheta0, naky, nspec))
allocate (vns(naky, negrid, nspec))
allocate (gtmp(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc))
vns = vnmult(2)*delvnew
!This is needed to to ensure the it,ik values we don't set aren't included
!in the integral (can also be enforced in integrate_moment routine)
if(any(kwork_filter)) call zero_array(gtmp)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! First get v0y0
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(iglo, it, ik, ie, is, isgn) &
!$OMP SHARED(g_lo, kwork_filter, gtmp, vnmult, vnew_E, aj0, g) &
!$OMP SCHEDULE(static)
do iglo = g_lo%llim_proc, g_lo%ulim_proc
ik = ik_idx(g_lo,iglo)
it = it_idx(g_lo,iglo)
if(kwork_filter(it,ik))cycle
ie = ie_idx(g_lo,iglo)
is = is_idx(g_lo,iglo)
do isgn = 1, 2
! v0 = nu_E E J0 f0
gtmp(:, isgn, iglo) = vnmult(2) * vnew_E(ik, ie, is) * aj0(:, iglo) &
* g(:, isgn, iglo)
end do
end do
!$OMP END PARALLEL DO
call integrate_moment (gtmp, v0y0, all_procs) ! v0y0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Get y1 = y0 - v0y0 * z0 / (1 + v0z0)
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(iglo, it, ik, is, isgn) &
!$OMP SHARED(g_lo, kwork_filter, g1, g, ieqzip, v0y0, bz0) &
!$OMP SCHEDULE(static)
do iglo = g_lo%llim_proc, g_lo%ulim_proc
it = it_idx(g_lo,iglo)
ik = ik_idx(g_lo,iglo)
if(kwork_filter(it,ik)) cycle
if(ieqzip(it,ik)) cycle
is = is_idx(g_lo,iglo)
do isgn = 1, 2
g1(:, isgn, iglo) = g(:, isgn, iglo) - v0y0(:, it, ik, is) * bz0(:, isgn, iglo)
end do
end do
!$OMP END PARALLEL DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Now get v1y1
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(iglo, it, ik, ie, is, isgn) &
!$OMP SHARED(g_lo, kwork_filter, gtmp, vns, vpa, aj0, g1) &
!$OMP SCHEDULE(static)
do iglo = g_lo%llim_proc, g_lo%ulim_proc
ik = ik_idx(g_lo,iglo)
it = it_idx(g_lo,iglo)
if(kwork_filter(it,ik))cycle
ie = ie_idx(g_lo,iglo)
is = is_idx(g_lo,iglo)
do isgn = 1, 2
! v1 = (nus-nud) vpa J0 f0
gtmp(:, isgn, iglo) = vns(ik, ie, is) * vpa(:, isgn, iglo) * aj0(:, iglo) &
* g1(:, isgn, iglo)
end do
end do
!$OMP END PARALLEL DO
call integrate_moment (gtmp, v1y1, all_procs) ! v1y1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Get y2 = y1 - v1y1 * s1 / (1 + v1s1)
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(iglo, it, ik, is, isgn) &
!$OMP SHARED(g_lo, kwork_filter, g1, ieqzip, v1y1, bs0) &
!$OMP SCHEDULE(static)
do iglo = g_lo%llim_proc, g_lo%ulim_proc
it = it_idx(g_lo,iglo)
ik = ik_idx(g_lo,iglo)
if(kwork_filter(it,ik)) cycle
if(ieqzip(it,ik)) cycle
is = is_idx(g_lo,iglo)
do isgn = 1, 2
g1(:, isgn, iglo) = g1(:, isgn, iglo) - v1y1(:, it, ik, is) * bs0(:, isgn, iglo)
end do
end do
!$OMP END PARALLEL DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Now get v2y2
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(iglo, it, ik, ie, il, is, isgn) &
!$OMP SHARED(g_lo, kwork_filter, gtmp, vns, energy, al, aj1, g1) &
!$OMP SCHEDULE(static)
do iglo = g_lo%llim_proc, g_lo%ulim_proc
ik = ik_idx(g_lo,iglo)
it = it_idx(g_lo,iglo)
if(kwork_filter(it,ik))cycle
ie = ie_idx(g_lo,iglo)
il = il_idx(g_lo,iglo)
is = is_idx(g_lo,iglo)
do isgn = 1, 2
! v2 = (nus-nud) vperp J1 f0
gtmp(:, isgn, iglo) = vns(ik, ie, is) * energy(ie) * al(il) * aj1(:, iglo) &
* g1(:, isgn, iglo)
end do
end do
!$OMP END PARALLEL DO
call integrate_moment (gtmp, v2y2, all_procs) ! v2y2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Finally get x = y2 - v2y2 * w2 / (1 + v2w2)
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(iglo, it, ik, is, isgn) &
!$OMP SHARED(g_lo, kwork_filter, g, g1, ieqzip, v2y2, bw0) &
!$OMP SCHEDULE(static)
do iglo = g_lo%llim_proc, g_lo%ulim_proc
it = it_idx(g_lo,iglo)
ik = ik_idx(g_lo,iglo)
if(kwork_filter(it,ik)) cycle
if(ieqzip(it,ik)) cycle
is = is_idx(g_lo,iglo)
do isgn = 1, 2
g(:, isgn, iglo) = g1(:, isgn, iglo) - v2y2(:, it, ik, is) * bw0(:, isgn, iglo)
end do
end do
!$OMP END PARALLEL DO
deallocate (vns, v0y0, v1y1, v2y2, gtmp)
end subroutine conserve_diffuse_standard_layout