nl_center Subroutine

private subroutine nl_center(gtmp, bd)

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

Arguments

Type IntentOptional Attributes Name
complex, intent(inout), dimension (-ntgrid:,:,g_lo%llim_proc:) :: gtmp
real, intent(in) :: bd

Contents

Source Code


Source Code

  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