Takes input array evaluated at theta grid points and overwrites it with array evaluated at cell centers note that there is an extra factor of 2 in output array
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
complex, | intent(inout), | dimension (-ntgrid:,:,g_lo%llim_proc:) | :: | gtmp | ||
real, | intent(in) | :: | bd |
subroutine nl_center (gtmp, bd)
!
!CMR, 13/10/2014
! Fixing some (cosmetic) issues from prior to R3021.
! (1) forbid(-ntgrid:ntgrid) refers to grid points at cell boundaries
! gtmp(-ntgrid:ntgrid-1) refers to cell centers
! => applying forbid to output gtmp did not zero the correct things!!!
! (2) totally trapped particles are special, so return source without upwinding is fine
! BUT source for other trapped particles should NOT enter forbidden region.
! if ig or ig+1 forbidden (i.e. ig+1 or ig is bounce point) now return gtmp(ig,:,iglo)=0
!
use gs2_layouts, only: g_lo, il_idx
use theta_grid, only: ntgrid
use le_grids, only: forbid, grid_has_trapped_particles, jend, is_ttp
use mp, only: mp_abort
implicit none
real, intent (in) :: bd
complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in out) :: gtmp
integer :: iglo, il, ig
logical :: trapped
trapped = grid_has_trapped_particles()
! factor of one-half appears elsewhere
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(iglo, il, ig) &
!$OMP SHARED(g_lo, nl_forbid_force_zero, forbid, gtmp, ntgrid, is_ttp, jend, trapped, bd) &
!$OMP SCHEDULE(static)
do iglo = g_lo%llim_proc, g_lo%ulim_proc
il = il_idx(g_lo, iglo)
!
!CMR, 7/10/2014:
! Incoming gtmp SHOULD vanish in forbidden region,
! New logical in nonlinear_terms_knobs namelist: nl_forbid_force_zero
! nl_forbid_force_zero =.t. : force zeroing (default)
! nl_forbid_force_zero =.f. : NO forced zeroing
! ie assume forbidden gtmp is zero on entry
!
if ( nl_forbid_force_zero ) then
! force spurious gtmp outside trapped boundary to be zero
where (forbid(:,il))
gtmp(:,1,iglo) = 0.0
gtmp(:,2,iglo) = 0.0
end where
endif
do ig = -ntgrid, ntgrid-1
!
!CMR, 7/10/2014:
! loop sets gtmp to value at upwinded cell center RIGHT of theta(ig)
! except for ttp where upwinding makes no sense!
!
! Note this cycle doesn't skip when ig is in the forbidden region
! even if the pitch angle is greater than the totally trapped one
! at this location. This ensures that whilst we leave the (allowed)
! totally trapped source alone we still zero the NL source for forbidden
! points. Previously we cycled simply if the pitch angle index was larger
! than the totally trapped one. When there are multiple totally trapped
! pitch angles this may not be equivalent.
if (is_ttp(ig, il)) cycle
! The below might be better off using forbid instead of checking jend
if ( trapped .and. ( il > jend(ig) .or. il > jend(ig+1)) ) then
!
!CMR, 7/10/2014:
! if either ig or ig+1 is forbidden, no source possible in a cell RIGHT of theta(ig)
! => gtmp(ig,1:2,iglo)=0
!
gtmp(ig,1:2,iglo) = 0.0
else
!
!CMR, 7/10/2014:
! otherwise ig and ig+1 BOTH allowed, and upwinding in cell RIGHT of theta(ig) is fine
!
gtmp(ig,1,iglo) = (1.+bd)*gtmp(ig+1,1,iglo) + (1.-bd)*gtmp(ig,1,iglo)
gtmp(ig,2,iglo) = (1.-bd)*gtmp(ig+1,2,iglo) + (1.+bd)*gtmp(ig,2,iglo)
endif
end do
end do
!$OMP END PARALLEL DO
end subroutine nl_center