Diagnose the poloidal distribution of the particle, angular momentum, and energy fluxes
JRB
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
complex, | intent(in), | dimension (-ntgrid:,:,:) | :: | phi | ||
real, | intent(out), | dimension (-ntgrid:,:,:,:) | :: | pflux_dist | ||
real, | intent(out), | dimension (-ntgrid:,:,:,:) | :: | vflux_par_dist | ||
real, | intent(out), | dimension (-ntgrid:,:,:,:) | :: | vflux_perp_dist | ||
real, | intent(out), | dimension (-ntgrid:,:,:,:) | :: | qflux_dist |
subroutine flux_dist (phi, pflux_dist, vflux_par_dist, vflux_perp_dist, qflux_dist)
!CMR, 15/1/08:
! Implemented Clemente Angioni's fix for fluxes by replacing g with gnew
! so fields and distribution function are evaluated self-consistently in time.
! This fixed unphysical oscillations in non-ambipolar particle fluxes
!
use species, only: spec
use theta_grid, only: ntgrid, bmag
use theta_grid, only: qval, shat, gds21, gds22
use kt_grids, only: theta0, aky
use le_grids, only: energy
use dist_fn_arrays, only: gnew, aj0, vpac, aj1, vperp2, g_work
use gs2_layouts, only: g_lo, ie_idx, is_idx, it_idx, ik_idx
use run_parameters, only:has_phi
use constants, only: zi
use theta_grid, only: Rplot, Bpol, rhoc
use array_utils, only: zero_array
implicit none
complex, dimension (-ntgrid:,:,:), intent (in) :: phi
real, dimension (-ntgrid:,:,:,:), intent (out) :: pflux_dist, vflux_par_dist, vflux_perp_dist, qflux_dist
integer :: it, ik, is, isgn, ig
integer :: iglo
call zero_array(pflux_dist) ; call zero_array(vflux_par_dist) ; call zero_array(vflux_perp_dist) ; call zero_array(qflux_dist)
if (has_phi) then
do isgn = 1, 2
g_work(:,isgn,:) = gnew(:,isgn,:)*aj0
end do
call get_flux_dist (g_work, phi, pflux_dist)
do iglo = g_lo%llim_proc, g_lo%ulim_proc
g_work(:,:,iglo) = g_work(:,:,iglo)*energy(ie_idx(g_lo,iglo),is_idx(g_lo,iglo))
end do
call get_flux_dist (g_work, phi, qflux_dist)
do isgn = 1, 2
do ig = -ntgrid, ntgrid
! Not clear that this should use the grid-centre v|| (vpac)
! as this is not defined at ig = ntgrid and all other quantities
! are evaluated on the grid points.
g_work(ig,isgn,:) = gnew(ig,isgn,:)*aj0(ig,:)*vpac(ig,isgn,:)*Rplot(ig)*sqrt(1.0-Bpol(ig)**2/bmag(ig)**2)
end do
end do
call get_flux_dist (g_work, phi, vflux_par_dist)
do iglo = g_lo%llim_proc, g_lo%ulim_proc
it = it_idx(g_lo,iglo)
ik = ik_idx(g_lo,iglo)
is = is_idx(g_lo,iglo)
do isgn = 1, 2
g_work(:,isgn,iglo) = -zi*aky(ik)*gnew(:,isgn,iglo)*aj1(:,iglo) &
*rhoc*(gds21+theta0(it,ik)*gds22)*vperp2(:,iglo)*spec(is)%smz/(qval*shat*bmag**2)
end do
end do
call get_flux_dist (g_work, phi, vflux_perp_dist)
end if
end subroutine flux_dist