Calculates moments at not guiding center coordinate
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
complex, | intent(in), | dimension(-ntgrid:,:,:) | :: | phinew | ||
complex, | intent(in), | dimension(-ntgrid:,:,:) | :: | bparnew | ||
complex, | intent(out) | :: | dens(-ntgrid:,:,:,:) | |||
complex, | intent(out) | :: | upar(-ntgrid:,:,:,:) | |||
complex, | intent(out) | :: | tpar(-ntgrid:,:,:,:) | |||
complex, | intent(out) | :: | tper(-ntgrid:,:,:,:) | |||
complex, | intent(out), | optional | :: | ntot(-ntgrid:,:,:,:) | ||
complex, | intent(out), | optional | :: | jpar(-ntgrid:,:,:) | ||
logical, | intent(in), | optional | :: | reset_coefficients |
subroutine getmoms_notgc (phinew, bparnew, dens, upar, tpar, tper, ntot, jpar, reset_coefficients)
use dist_fn_arrays, only: vpa, vperp2, aj0, aj1, gnew, g_work
use gs2_layouts, only: g_lo, is_idx, ik_idx, it_idx, ie_idx
use species, only: nspec, spec, nonmaxw_corr
use theta_grid, only: ntgrid
use kt_grids, only: nakx => ntheta0, naky
use le_grids, only: integrate_moment
use optionals, only: get_option_with_default
implicit none
logical, parameter :: full_arr=moment_to_allprocs
complex, intent (out) :: &
& dens(-ntgrid:,:,:,:), upar(-ntgrid:,:,:,:), &
& tpar(-ntgrid:,:,:,:), tper(-ntgrid:,:,:,:)
complex, intent (out), optional :: ntot(-ntgrid:,:,:,:)
complex, intent (out), optional :: jpar(-ntgrid:,:,:)
complex, dimension(-ntgrid:,:,:), intent(in) :: phinew, bparnew
logical, intent(in), optional :: reset_coefficients
real, dimension(:, :, :, :), allocatable, save :: mom_coeff
real, dimension(:, :, :), allocatable, save :: mom_coeff_npara, mom_coeff_nperp
real, dimension(:, :, :), allocatable, save :: mom_coeff_tpara, mom_coeff_tperp
real, dimension(:, :, :), allocatable, save :: mom_shift_para, mom_shift_perp
integer :: isgn, iglo, is, ie
real :: a, b, tpar2, tper2
integer :: it, ik, ig
! returns moment integrals to PE 0
if (get_option_with_default(reset_coefficients, .false.)) then
if (allocated(mom_coeff)) deallocate(mom_coeff)
if (allocated(mom_coeff_npara)) deallocate(mom_coeff_npara)
if (allocated(mom_coeff_nperp)) deallocate(mom_coeff_nperp)
if (allocated(mom_coeff_tpara)) deallocate(mom_coeff_tpara)
if (allocated(mom_coeff_tperp)) deallocate(mom_coeff_tperp)
if (allocated(mom_shift_para)) deallocate(mom_shift_para)
if (allocated(mom_shift_perp)) deallocate(mom_shift_perp)
end if
! Note we mark mom_coeff* as 'save' so that the result of this call is cached
! and not recalculated each time. This can help save overhead and should be
! fine. The result of this call should only depend on the velocity grids, the
! wavenumbers and the species parameters (and the derived Bessel functions).
! These should not change during a single standard run, but if we wish to enable
! more exotic behaviour we might want to remove the save from above declarations and
! recalculate these coefficients each time. For now just add an optional flag to
! request recalculation.
call init_mom_coeff(mom_coeff, mom_coeff_npara, mom_coeff_nperp, mom_coeff_tpara, &
mom_coeff_tperp, mom_shift_para, mom_shift_perp)
! not guiding center n_total
if(present(ntot)) then
do iglo = g_lo%llim_proc, g_lo%ulim_proc
is = is_idx(g_lo,iglo)
ik = ik_idx(g_lo,iglo)
it = it_idx(g_lo,iglo)
ie = ie_idx(g_lo,iglo)
do isgn = 1, 2
g_work(:,isgn,iglo) = aj0(:,iglo) * gnew(:,isgn,iglo)
end do
do isgn = 1, 2
g_work(:,isgn,iglo) = g_work(:,isgn,iglo) + phinew(:,it,ik) &
& *(aj0(:,iglo)**2-1.0) * spec(is)%zt * nonmaxw_corr(ie,is)
end do
do isgn = 1, 2
g_work(:,isgn,iglo) = g_work(:,isgn,iglo) &
& + 2.*vperp2(:,iglo)*aj1(:,iglo)*aj0(:,iglo) &
& * bparnew(:,it,ik) * nonmaxw_corr(ie,is)
end do
end do
call integrate_moment (g_work, ntot, moment_to_allprocs,full_arr)
endif
! not guiding center density
do iglo = g_lo%llim_proc, g_lo%ulim_proc
do isgn = 1, 2
g_work(:,isgn,iglo) = aj0(:,iglo)*gnew(:,isgn,iglo)
end do
end do
call integrate_moment (g_work, dens, moment_to_allprocs,full_arr)
! not guiding center upar
do iglo = g_lo%llim_proc, g_lo%ulim_proc
do isgn = 1, 2
g_work(:,isgn,iglo) = aj0(:,iglo)*vpa(:,isgn,iglo)*gnew(:,isgn,iglo)
end do
end do
call integrate_moment (g_work, upar, moment_to_allprocs,full_arr)
! not guiding center tpar
do iglo = g_lo%llim_proc, g_lo%ulim_proc
do isgn = 1, 2
g_work(:,isgn,iglo) = 2.*aj0(:,iglo)*vpa(:,isgn,iglo)**2*gnew(:,isgn,iglo)
end do
end do
call integrate_moment (g_work, tpar, moment_to_allprocs,full_arr)
! not guiding center tperp
do iglo = g_lo%llim_proc, g_lo%ulim_proc
do isgn = 1, 2
g_work(:,isgn,iglo) = 2.*vperp2(:,iglo)*aj1(:,iglo)*gnew(:,isgn,iglo)
end do
end do
call integrate_moment (g_work, tper, moment_to_allprocs,full_arr)
do ig=-ntgrid,ntgrid
do it=1,nakx
do ik=1,naky
do is=1,nspec
tpar2=tpar(ig,it,ik,is) &
& - dens(ig,it,ik,is)*mom_coeff_npara(it,ik,is)
tper2=tper(ig,it,ik,is) &
& - dens(ig,it,ik,is)*mom_coeff_nperp(it,ik,is)
a=mom_coeff_tperp(it,ik,is)
b=mom_coeff_tpara(it,ik,is)
tpar(ig,it,ik,is)=( tpar2-a*tper2)/(1.-a*b)
tper(ig,it,ik,is)=(-b*tpar2+ tper2)/(1.-a*b)
end do
end do
end do
end do
do is=1,nspec
dens(:,:,:,is)=dens(:,:,:,is)*spec(is)%dens
upar(:,:,:,is)=upar(:,:,:,is)*spec(is)%stm
tpar(:,:,:,is)=tpar(:,:,:,is)*spec(is)%temp
tper(:,:,:,is)=tper(:,:,:,is)*spec(is)%temp
end do
if(present(jpar)) then
jpar(:,:,:)=cmplx(0.,0.)
do is=1,nspec
jpar(:,:,:)=jpar(:,:,:)+spec(is)%z*spec(is)%dens*upar(:,:,:,is)
end do
endif
end subroutine getmoms_notgc