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 :: ncell
type(ParticleTrack), dimension(beam_grid%ntrack) :: tracks
integer(Int64), dimension(3) :: nl_birth
type(LocalProfiles) :: plasma
real(Float64), dimension(nlevs) :: states, dens
real(Float64) :: photons, iflux, r_nlaunch
integer(Int32), dimension(3) :: ind
real(Float64), dimension(1) :: randomu
integer, dimension(1) :: randi
logical :: err
if(inputs%verbose.ge.1) then
write(*,'(T6,"# of markers: ",i9)') inputs%n_nbi
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)
r_nlaunch = 1.d0/nlaunch
cnt = 0
!$OMP PARALLEL DO schedule(guided) &
!$OMP& private(vnbi,rnbi,tracks,ncell,plasma,nl_birth,randi, &
!$OMP& states,dens,iflux,photons,neut_type,jj,ii,kk,ind,err)
loop_over_markers: do ii=1,inputs%n_nbi
if(inputs%calc_birth.ge.1) then
!! Determine the type of birth particle (1, 2, or 3)
nl_birth = 0
do kk=1,inputs%n_birth
call randind(nbi%current_fractions,randi)
nl_birth(randi(1)) = nl_birth(randi(1)) + 1
enddo
endif
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,ncell)
if(ncell.eq.0) cycle energy_fractions
!! Solve collisional radiative model along track
states=0.d0
states(1)=nneutrals*nbi%current_fractions(neut_type)/beam_grid%dv
loop_along_track: do jj=1,ncell
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)
tracks(jj)%dens = dens*r_nlaunch
tracks(jj)%flux = (iflux - sum(states))*beam_grid%dv*r_nlaunch
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*r_nlaunch,neut_type)
endif
enddo loop_along_track
call store_neutrals(tracks, ncell, neut_type)
if(inputs%calc_birth.ge.1) then
!! Sample according to deposited flux along neutral trajectory
!$OMP CRITICAL(ndmc_birth)
do kk=1,nl_birth(neut_type)
call randind(tracks(1:ncell)%flux,randi)
call randu(randomu)
birth%neut_type(birth%cnt) = neut_type
birth%ind(:,birth%cnt) = tracks(randi(1))%ind
birth%vi(:,birth%cnt) = vnbi
birth%ri(:,birth%cnt) = tracks(randi(1))%pos + &
vnbi*(tracks(randi(1))%time*(randomu(1)-0.5))
birth%cnt = birth%cnt+1
enddo
!$OMP END CRITICAL(ndmc_birth)
endif
enddo energy_fractions
if (inputs%verbose.ge.2)then
cnt = cnt + 1
WRITE(*,'(f7.2,"% completed",a,$)') 100*cnt/nlaunch,char(13)
endif
enddo loop_over_markers
!$OMP END PARALLEL DO
if(nbi_outside.gt.0)then
if(inputs%verbose.ge.0) then
write(*,'(T4,a, f6.2)') 'Percent of markers outside the grid: ', &
100.*nbi_outside/(3.*inputs%n_nbi)
endif
if(sum(neut%dens).eq.0) stop 'Beam does not intersect the grid!'
endif
end subroutine ndmc