xigridset Subroutine

private subroutine xigridset()

Uses

FIXME : Add documentation

Arguments

None

Contents

Source Code


Source Code

  subroutine xigridset
    use theta_grid, only: ntgrid, bmag
    implicit none

    integer :: ig, je, ixi, il

    if (.not. allocated(xi)) allocate (xi(2*nlambda, -ntgrid:ntgrid))
    if (.not. allocated(ixi_to_il)) allocate (ixi_to_il(2*nlambda, -ntgrid:ntgrid))
    if (.not. allocated(ixi_to_isgn)) allocate (ixi_to_isgn(2*nlambda, -ntgrid:ntgrid))
    if (.not. allocated(wxi)) allocate (wxi(2*nlambda, -ntgrid:ntgrid))

    ! define array 'sgn' that returns sign of vpa associated with isgn=1,2
    sgn(1) = 1
    sgn(2) = -1

    ! xi is vpa / v and goes from 1 --> -1
    xi = 0.0
    do ig = -ntgrid, ntgrid
       xi(:jend(ig), ig) = sqrt(max(1.0 - al(:jend(ig))*spread(bmag(ig),1,jend(ig)),0.0))
       xi(jend(ig)+1:2*jend(ig)-1, ig) = -xi(jend(ig)-1:1:-1, ig)
    end do
    
    ! get mapping from ixi (which runs between 1 and 2*nlambda) and il (runs between 1 and nlambda)
    do ig = -ntgrid, ntgrid
       je = jend(ig)
       ! if no trapped particles
       if (je == 0) then
          do ixi = 1, 2*nlambda
             if (ixi > nlambda) then
                ixi_to_isgn(ixi, ig) = 2
                ixi_to_il(ixi, ig) = 2*nlambda - ixi + 1
             else
                ixi_to_isgn(ixi, ig) = 1
                ixi_to_il(ixi, ig) = ixi
             end if
          end do
       else
!CMR, 1/11/2013:
! Sketch of how ixi=>il mapping is arranged
!===============================================================================================
!  ixi=   1, ... , je-1, je, je+1, ... , 2je-1, || 2je, 2je+1, ... , nl+je, nl+je+1, ... ,  2nl
!   il=   1, ... , je-1, je, je-1, ... ,     1, ||  je,  je+1, ... ,    nl,      nl, ... , je+1
! isgn=   1, ... ,    1,  2,    2, ... ,     2, ||   1,     1, ... ,     1,       2, ... ,    2
!         particles passing through             ||  je,   + forbidden trapped velocity space
!         nb only need one isigma for je, as v||=0 at the bounce point
!===============================================================================================

          do ixi = 1, 2*nlambda
             if (ixi >= nlambda + je + 1) then
                ixi_to_isgn(ixi, ig) = 2
                ixi_to_il(ixi, ig) = 2*nlambda + je + 1 - ixi
             else if (ixi >= 2*je) then
                ixi_to_isgn(ixi, ig) = 1
                ixi_to_il(ixi, ig) = ixi - je
             else if (ixi >= je) then
                ixi_to_isgn(ixi, ig) = 2
                ixi_to_il(ixi, ig) = 2*je - ixi
             else
                ixi_to_isgn(ixi, ig) = 1
                ixi_to_il(ixi, ig) = ixi
             end if
          end do
       end if
    end do

    wxi = 0.0
    do ig = -ntgrid, ntgrid
       do ixi = 1, 2*nlambda
          il = ixi_to_il(ixi, ig)
          wxi(ixi, ig) = wl(ig, il)
       end do
    end do
  end subroutine xigridset