npa_f Subroutine

  1. fidasim.f90
  2. libfida
  3. npa_f

public subroutine npa_f()

Calculate NPA flux using a fast-ion distribution function F(E,p,r,z)

Arguments

None

Calls

proc~~npa_f~~CallsGraph proc~npa_f npa_f proc~npa_gyro_range npa_gyro_range proc~npa_f->proc~npa_gyro_range proc~store_npa store_npa proc~npa_f->proc~store_npa proc~get_nlaunch get_nlaunch proc~npa_f->proc~get_nlaunch proc~get_plasma get_plasma proc~npa_f->proc~get_plasma proc~get_indices get_indices proc~npa_f->proc~get_indices proc~hit_npa_detector hit_npa_detector proc~npa_f->proc~hit_npa_detector proc~gyro_trajectory gyro_trajectory proc~npa_f->proc~gyro_trajectory proc~gyro_surface gyro_surface proc~npa_f->proc~gyro_surface proc~attenuate attenuate proc~npa_f->proc~attenuate proc~get_beam_cx_rate get_beam_cx_rate proc~npa_f->proc~get_beam_cx_rate proc~mc_fastion mc_fastion proc~npa_f->proc~mc_fastion proc~approx_eq approx_eq proc~npa_gyro_range->proc~approx_eq proc~gyro_range gyro_range proc~npa_gyro_range->proc~gyro_range proc~get_fields get_fields proc~store_npa->proc~get_fields proc~xyz_to_uvw xyz_to_uvw proc~store_npa->proc~xyz_to_uvw proc~in_plasma in_plasma proc~get_plasma->proc~in_plasma proc~line_plane_intersect line_plane_intersect proc~hit_npa_detector->proc~line_plane_intersect proc~in_boundary in_boundary proc~hit_npa_detector->proc~in_boundary proc~attenuate->proc~get_plasma proc~colrad colrad proc~attenuate->proc~colrad proc~get_beam_cx_rate->proc~get_plasma proc~bt_cx_rates bt_cx_rates proc~get_beam_cx_rate->proc~bt_cx_rates proc~bb_cx_rates bb_cx_rates proc~get_beam_cx_rate->proc~bb_cx_rates proc~mc_fastion->proc~get_fields interface~randind randind proc~mc_fastion->interface~randind proc~randu randu proc~mc_fastion->proc~randu proc~get_distribution get_distribution proc~mc_fastion->proc~get_distribution proc~get_fields->proc~in_plasma proc~calc_perp_vectors calc_perp_vectors proc~get_fields->proc~calc_perp_vectors interface~interpol_coeff interpol_coeff proc~bt_cx_rates->interface~interpol_coeff omp_get_thread_num omp_get_thread_num proc~randu->omp_get_thread_num proc~rng_uniform rng_uniform proc~randu->proc~rng_uniform proc~get_distribution->proc~xyz_to_uvw interface~interpol interpol proc~get_distribution->interface~interpol proc~in_plasma->proc~xyz_to_uvw proc~in_plasma->interface~interpol_coeff proc~eigen eigen proc~colrad->proc~eigen proc~get_rate_matrix get_rate_matrix proc~colrad->proc~get_rate_matrix proc~matinv matinv proc~colrad->proc~matinv proc~gyro_range->proc~line_plane_intersect proc~gyro_range->proc~in_boundary proc~boundary_edge boundary_edge proc~gyro_range->proc~boundary_edge proc~line_gyro_surface_intersect line_gyro_surface_intersect proc~gyro_range->proc~line_gyro_surface_intersect proc~gyro_surface_coordinates gyro_surface_coordinates proc~gyro_range->proc~gyro_surface_coordinates proc~in_gyro_surface in_gyro_surface proc~gyro_range->proc~in_gyro_surface proc~bb_cx_rates->interface~interpol_coeff proc~hqr2 hqr2 proc~eigen->proc~hqr2 proc~balback balback proc~eigen->proc~balback proc~elmtrans elmtrans proc~eigen->proc~elmtrans proc~elmhes elmhes proc~eigen->proc~elmhes proc~balance balance proc~eigen->proc~balance proc~interpol1d_coeff interpol1D_coeff interface~interpol_coeff->proc~interpol1d_coeff proc~interpol2d_coeff_arr interpol2D_coeff_arr interface~interpol_coeff->proc~interpol2d_coeff_arr proc~interpol2d_coeff interpol2D_coeff interface~interpol_coeff->proc~interpol2d_coeff proc~interpol1d_coeff_arr interpol1D_coeff_arr interface~interpol_coeff->proc~interpol1d_coeff_arr proc~get_rate_matrix->interface~interpol_coeff proc~interpol2d_arr interpol2D_arr interface~interpol->proc~interpol2d_arr proc~interpol2d_2d_arr interpol2D_2D_arr interface~interpol->proc~interpol2d_2d_arr proc~interpol1d_arr interpol1D_arr interface~interpol->proc~interpol1d_arr proc~ludcmp ludcmp proc~matinv->proc~ludcmp proc~interpol2d_arr->interface~interpol_coeff proc~hqrvec hqrvec proc~hqr2->proc~hqrvec proc~outerprod outerprod proc~ludcmp->proc~outerprod proc~interpol2d_2d_arr->interface~interpol_coeff proc~interpol1d_arr->interface~interpol_coeff proc~interpol2d_coeff_arr->proc~interpol2d_coeff proc~interpol1d_coeff_arr->proc~interpol1d_coeff

Called by

proc~~npa_f~~CalledByGraph proc~npa_f npa_f program~fidasim fidasim program~fidasim->proc~npa_f

Contents

Source Code

npa_f

Source Code

subroutine npa_f
    !+ Calculate NPA flux using a fast-ion distribution function F(E,p,r,z)
    integer :: i,j,k,det,ip
    integer(Int64) :: iion
    real(Float64), dimension(3) :: rg,ri,rf,vi
    integer, dimension(3) :: ind,pind
    real(Float64) :: denf
    integer, dimension(3,beam_grid%ngrid) :: pcell
    type(LocalProfiles) :: plasma
    type(LocalEMFields) :: fields
    type(GyroSurface) :: gs
    real(Float64), dimension(2,4) :: gyrange
    integer, dimension(4) :: neut_types=[1,2,3,4]
    real(Float64), dimension(nlevs) :: rates
    real(Float64), dimension(nlevs) :: states
    real(Float64) :: flux, theta, dtheta, eb, ptch

    integer :: inpa,pcnt,ichan,nrange,ir
    real(Float64) :: papprox_tot, maxcnt, cnt, inv_maxcnt
    real(Float64), dimension(beam_grid%nx,beam_grid%ny,beam_grid%nz) :: papprox, nlaunch

    papprox=0.d0
    papprox_tot=0.d0

    pcnt=1
    do k=1,beam_grid%nz
        do j=1,beam_grid%ny
            do i=1,beam_grid%nx
                ind =[i,j,k]
                call get_plasma(plasma,ind=ind)
                papprox(i,j,k)=(sum(neut%dens(:,nbif_type,i,j,k)) + &
                                sum(neut%dens(:,nbih_type,i,j,k)) + &
                                sum(neut%dens(:,nbit_type,i,j,k)) + &
                                sum(neut%dens(:,halo_type,i,j,k)))* &
                                plasma%denf

                if(papprox(i,j,k).gt.0) then
                    pcell(:,pcnt)= ind
                    pcnt = pcnt + 1
                endif

                if(plasma%in_plasma) papprox_tot=papprox_tot+papprox(i,j,k)
            enddo
        enddo
    enddo
    pcnt = pcnt - 1
    maxcnt=real(pcnt)
    inv_maxcnt = 100.0/maxcnt

    call get_nlaunch(inputs%n_npa,papprox,papprox_tot,nlaunch)
    if(inputs%verbose.ge.1) then
        write(*,'(T6,"# of markers: ",i12)') int(sum(nlaunch),Int64)
    endif

    !! Loop over all cells that can contribute to NPA signal
    cnt=0.d0
    loop_over_cells: do ip = 1, int(pcnt)
        i = pcell(1,ip)
        j = pcell(2,ip)
        k = pcell(3,ip)
        ind = [i, j, k]
        !$OMP PARALLEL DO schedule(guided) private(iion,ichan,fields,nrange,gyrange, &
        !$OMP& pind,vi,ri,rf,det,plasma,rates,states,flux,denf,eb,ptch,gs,ir,theta,dtheta)
        loop_over_fast_ions: do iion=1,int(nlaunch(i, j, k),Int64)
            !! Sample fast ion distribution for energy and pitch
            call mc_fastion(ind, fields, eb, ptch, denf)
            if(denf.eq.0.0) cycle loop_over_fast_ions

            call gyro_surface(fields, eb, ptch, gs)

            detector_loop: do ichan=1,npa_chords%nchan
                call npa_gyro_range(ichan, gs, gyrange, nrange)
                if(nrange.eq.0) cycle detector_loop
                gyro_range_loop: do ir=1,nrange
                    dtheta = gyrange(2,ir)
                    theta = gyrange(1,ir) + 0.5*dtheta
                    call gyro_trajectory(gs, theta, ri, vi)

                    !! Check if particle hits a NPA detector
                    call hit_npa_detector(ri, vi ,det, rf, ichan)
                    if(det.ne.ichan) then
                        if (inputs%verbose.ge.0)then
                            write(*,*) "NPA_F: Missed Detector ",ichan
                        endif
                        cycle gyro_range_loop
                    endif

                    !! Get beam grid indices at ri
                    call get_indices(ri,pind)

                    !! Calculate CX probability with beam and halo neutrals
                    call get_beam_cx_rate(pind,ri,vi,beam_ion,neut_types,rates)
                    if(sum(rates).le.0.) cycle gyro_range_loop

                    !! Weight CX rates by ion source density
                    states=rates*denf

                    !! Attenuate states as the particle move through plasma
                    call attenuate(ri,rf,vi,states)

                    !! Store NPA Flux
                    flux = (dtheta/(2*pi))*sum(states)*beam_grid%dv/nlaunch(i,j,k)
                    call store_npa(det,ri,rf,vi,flux)
                enddo gyro_range_loop
            enddo detector_loop
        enddo loop_over_fast_ions
        !$OMP END PARALLEL DO
        cnt=cnt+1
        if (inputs%verbose.ge.2)then
            WRITE(*,'(f7.2,"% completed",a,$)') cnt*inv_maxcnt,char(13)
        endif
    enddo loop_over_cells
    if(inputs%verbose.ge.1) then
        write(*,'(T4,"Number of NPA particles that hit a detector: ",i8)') npa%npart
    endif

end subroutine npa_f