flux_dist Subroutine

public subroutine flux_dist(phi, pflux_dist, vflux_par_dist, vflux_perp_dist, qflux_dist)

Diagnose the poloidal distribution of the particle, angular momentum, and energy fluxes

JRB

Arguments

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

Contents

Source Code


Source Code

  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