getmoms_notgc Subroutine

public subroutine getmoms_notgc(phinew, bparnew, dens, upar, tpar, tper, ntot, jpar, reset_coefficients)

Calculates moments at not guiding center coordinate

Arguments

Type IntentOptional 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

Contents

Source Code


Source Code

  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