Calculates neutral beam deposition and spectra
subroutine ndmc
!+ Calculates neutral beam deposition and spectra
integer :: neut_type !! full half third energy
real(Float64) :: nlaunch !! nr. of markers
real(Float64) :: nneutrals !! # NBI particles
real(Float64), dimension(3) :: vnbi !! velocities(full..)
real(Float64), dimension(3) :: rnbi !! initial position
integer(Int64) :: jj, ii, kk, cnt
integer :: ntrack
type(ParticleTrack), dimension(beam_grid%ntrack) :: tracks
type(LocalProfiles) :: plasma
type(LocalEMFields) :: fields
real(Float64), dimension(nlevs) :: states, dens
real(Float64) :: photons, iflux,flux_tot
integer(Int32), dimension(3) :: ind
real(Float64), dimension(3) :: ri,ri_gc,r_gyro
real(Float64), dimension(1) :: randomu
integer, dimension(1) :: randi
logical :: err
if(inputs%verbose.ge.1) then
write(*,'(T6,"# of markers: ",i10)') inputs%n_nbi
if(inputs%calc_birth.ge.1) then
write(*,'(T6,"# of birth markers: 3 x",i10)') int(inputs%n_nbi*inputs%n_birth)
endif
endif
!! # of injected neutrals = NBI power/energy_per_particle
nneutrals=1.d6*nbi%pinj/ (1.d3*nbi%einj*e0 &
*( nbi%current_fractions(1) &
+ nbi%current_fractions(2)/2.d0 &
+ nbi%current_fractions(3)/3.d0 ) )
nlaunch=real(inputs%n_nbi)
!$OMP PARALLEL DO schedule(guided) &
!$OMP& private(vnbi,rnbi,tracks,ntrack,plasma,fields,randi,flux_tot, &
!$OMP& states,dens,iflux,photons,neut_type,jj,ii,kk,ind,err,ri,ri_gc,r_gyro)
loop_over_markers: do ii=istart,inputs%n_nbi,istep
energy_fractions: do neut_type=1,3
!! (type = 1: full energy, =2: half energy, =3: third energy
call mc_nbi(vnbi,neut_type,rnbi,err)
if(err) cycle energy_fractions
call track(rnbi,vnbi,tracks,ntrack)
if(ntrack.eq.0) cycle energy_fractions
!! Solve collisional radiative model along track
flux_tot = 0.d0
states=0.d0
states(1)=nneutrals*nbi%current_fractions(neut_type)/beam_grid%dv
loop_along_track: do jj=1,ntrack
iflux = sum(states)
ind = tracks(jj)%ind
call get_plasma(plasma,pos=tracks(jj)%pos)
call colrad(plasma,beam_ion,vnbi,tracks(jj)%time,states,dens,photons)
call store_neutrals(ind,neut_type,dens/nlaunch,vnbi)
tracks(jj)%flux = (iflux - sum(states))/nlaunch
flux_tot = flux_tot + tracks(jj)%flux*beam_grid%dv
if(inputs%calc_birth.ge.1) then
call store_births(ind,neut_type,tracks(jj)%flux)
endif
if((photons.gt.0.d0).and.(inputs%calc_bes.ge.1)) then
call store_bes_photons(tracks(jj)%pos,vnbi,photons/nlaunch,neut_type)
endif
enddo loop_along_track
if((inputs%calc_birth.ge.1).and.(flux_tot.gt.0.d0)) then
!! Sample according to deposited flux along neutral trajectory
!$OMP CRITICAL(ndmc_birth)
do kk=1,inputs%n_birth
call randind(tracks(1:ntrack)%flux,randi)
call randu(randomu)
birth%part(birth%cnt)%neut_type = neut_type
birth%part(birth%cnt)%energy = nbi%einj/real(neut_type)
birth%part(birth%cnt)%weight = flux_tot/inputs%n_birth
birth%part(birth%cnt)%ind = tracks(randi(1))%ind
birth%part(birth%cnt)%vi = vnbi
ri = tracks(randi(1))%pos + vnbi*(tracks(randi(1))%time*(randomu(1)-0.5))
birth%part(birth%cnt)%ri = ri
call get_fields(fields,pos=ri)
birth%part(birth%cnt)%pitch = dot_product(fields%b_norm,vnbi/norm2(vnbi))
call gyro_step(vnbi,fields,r_gyro)
birth%part(birth%cnt)%ri_gc = ri + r_gyro
birth%cnt = birth%cnt + 1
enddo
!$OMP END CRITICAL(ndmc_birth)
endif
enddo energy_fractions
enddo loop_over_markers
!$OMP END PARALLEL DO
#ifdef _MPI
!! Combine beam neutrals
call parallel_sum(neut%full)
call parallel_sum(neut%half)
call parallel_sum(neut%third)
call parallel_sum(nbi_outside)
if(inputs%calc_birth.ge.1) then
call parallel_sum(birth%dens)
endif
!! Combine spectra
if(inputs%calc_bes.ge.1) then
call parallel_sum(spec%full)
call parallel_sum(spec%half)
call parallel_sum(spec%third)
endif
#endif
if(nbi_outside.gt.0)then
if(inputs%verbose.ge.1) then
write(*,'(T4,a, f6.2,a)') 'Percent of markers outside the grid: ', &
100.*nbi_outside/(3.*inputs%n_nbi),'%'
endif
if(sum(neut%full).eq.0) stop 'Beam does not intersect the grid!'
endif
end subroutine ndmc