Calculates halo neutral density and spectra
subroutine halo
!+ Calculates halo neutral density and spectra
integer :: ic,i,j,k,ncell
integer(Int64) :: ihalo !! counter
real(Float64), dimension(3) :: ri !! start position
real(Float64), dimension(3) :: vihalo!! velocity bulk plasma ion
integer,dimension(3) :: ind,tind !! actual cell
!! Determination of the CX probability
type(LocalProfiles) :: plasma
real(Float64), dimension(nlevs) :: denn !! neutral dens (n=1-4)
real(Float64), dimension(nlevs) :: rates !! CX rates
!! Collisiional radiative model along track
real(Float64), dimension(nlevs) :: states ! Density of n-states
integer :: ntrack
type(ParticleTrack), dimension(beam_grid%ntrack) :: tracks !! Particle Tracks
integer :: jj !! counter along track
real(Float64) :: tot_denn, photons !! photon flux
integer, dimension(beam_grid%ngrid) :: cell_ind
real(Float64), dimension(beam_grid%nx,beam_grid%ny,beam_grid%nz) :: papprox
integer(Int32), dimension(beam_grid%nx,beam_grid%ny,beam_grid%nz) :: nlaunch
real(Float64), dimension(nlevs,beam_grid%nx,beam_grid%ny,beam_grid%nz) :: dens_prev
real(Float64), dimension(:,:,:,:,:),allocatable :: dens_cur
real(Float64) :: local_iter_dens
integer :: n_slices, cur_slice, is
#ifdef _OMP
integer, external :: omp_get_max_threads, omp_get_thread_num
#endif
!! Halo iteration
integer(Int64) :: hh, n_halo !! counters
real(Float64) :: dcx_dens, halo_iteration_dens,seed_dcx
integer :: prev_type ! previous iteration
integer :: cur_type ! current iteration
real(Float64) :: fi_correction
prev_type = fida_type
cur_type = brems_type
cur_slice = 1
n_slices = 1
#ifdef _OMP
n_slices = omp_get_max_threads()
#endif
allocate(dens_cur(nlevs,beam_grid%nx,beam_grid%ny,beam_grid%nz,n_slices))
dens_prev = 0.d0
dens_cur = 0.d0
dcx_dens = halo_iter_dens(dcx_type)
halo_iter_dens(prev_type) = dcx_dens
if(dcx_dens.eq.0) then
if(inputs%verbose.ge.0) then
write(*,'(a)') 'HALO: Density of DCX-neutrals is zero'
endif
stop
endif
dens_prev = neut%dcx
n_halo = inputs%n_halo
seed_dcx = 1.0
iterations: do hh=1,200
papprox=0.d0
tot_denn=0.d0
halo_iter_dens(cur_type) = 0.d0
do ic=1,beam_grid%ngrid
call ind2sub(beam_grid%dims,ic,ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
call get_plasma(plasma,ind=ind)
if(.not.plasma%in_plasma) cycle
tot_denn = sum(dens_prev(:,i,j,k))
papprox(i,j,k)= tot_denn*(plasma%denp-plasma%denf)
enddo
cell_ind = 0
ncell = 0
do ic=1,beam_grid%ngrid
call ind2sub(beam_grid%dims,ic,ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
if(papprox(i,j,k).gt.0.0) then
ncell = ncell + 1
cell_ind(ncell) = ic
endif
enddo
call get_nlaunch(n_halo, papprox, nlaunch)
if(inputs%verbose.ge.1) then
write(*,'(T6,"# of markers: ",i10," --- Seed/DCX: ",f5.3)') sum(nlaunch), seed_dcx
endif
local_iter_dens = halo_iter_dens(cur_type)
!$OMP PARALLEL DO schedule(dynamic,1) private(i,j,k,ic,ihalo,ind,vihalo, &
!$OMP& ri,tracks,ntrack,rates,denn,states,jj,photons,plasma,tind, cur_slice,fi_correction) &
!$OMP& reduction(+: local_iter_dens )
loop_over_cells: do ic=istart,ncell,istep
#ifdef _OMP
cur_slice = omp_get_thread_num()+1
#endif
call ind2sub(beam_grid%dims,cell_ind(ic),ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
!! Loop over the markers
loop_over_halos: do ihalo=1, nlaunch(i,j,k)
!! Calculate ri,vhalo and track
call mc_halo(ind,vihalo,ri)
call track(ri,vihalo,tracks,ntrack)
if(ntrack.eq.0)cycle loop_over_halos
!! Get plasma parameters at particle location
call get_plasma(plasma, pos=ri)
!! Calculate CX probability
tind = tracks(1)%ind
call bt_cx_rates(plasma, dens_prev(:,tind(1),tind(2),tind(3)), vihalo, thermal_ion, rates)
if(sum(rates).le.0.)cycle loop_over_halos
!! Get plasma parameters at mean point in cell
call get_plasma(plasma,pos=tracks(1)%pos)
!! Weight CX rates by ion source density
states = rates*plasma%denp
fi_correction = max((plasma%denp - plasma%denf)/plasma%denp,0.d0)
loop_along_track: do jj=1,ntrack
call get_plasma(plasma,pos=tracks(jj)%pos)
if(.not.plasma%in_plasma) exit loop_along_track
call colrad(plasma,thermal_ion,vihalo,tracks(jj)%time,states,denn,photons)
!! Store Neutrals
tind = tracks(jj)%ind
dens_cur(:,tind(1),tind(2),tind(3),cur_slice) = &
dens_cur(:,tind(1),tind(2),tind(3),cur_slice) + denn/nlaunch(i,j,k)
local_iter_dens = &
local_iter_dens + sum(denn)/nlaunch(i,j,k)
if((photons.gt.0.d0).and.(inputs%calc_halo.ge.1)) then
photons = fi_correction*photons !! Correct for including fast-ions in states
call store_bes_photons(tracks(jj)%pos,vihalo,photons/nlaunch(i,j,k),halo_type)
endif
enddo loop_along_track
enddo loop_over_halos
enddo loop_over_cells
!$OMP END PARALLEL DO
halo_iter_dens(cur_type) = local_iter_dens
#ifdef _OMP
!$OMP PARALLEL DO private(i,j,is)
do i=1,beam_grid%nz
do j=1,beam_grid%ny
do is = 2, n_slices
dens_cur(:,:,j,i,1) = dens_cur(:,:,j,i,1) + dens_cur(:,:,j,i,is)
enddo
enddo
enddo
#endif
! at this point, dens_cur(*,1) contains all the info
#ifdef _MPI
!! Combine densities
call parallel_sum(dens_cur(:,:,:,:,1))
call parallel_sum(halo_iter_dens(cur_type))
#endif
if(halo_iter_dens(cur_type)/halo_iter_dens(prev_type).gt.1.0) then
write(*,'(a)') "HALO: Halo generation density exceeded seed density. This shouldn't happen."
exit iterations
endif
halo_iteration_dens = halo_iter_dens(cur_type)
halo_iter_dens(prev_type) = halo_iter_dens(cur_type)
!$OMP PARALLEL DO private(i)
do i=1,beam_grid%nz
neut%halo(:,:,:,i) = neut%halo(:,:,:,i) + dens_cur(:,:,:,i,1)
dens_prev(:,:,:,i) = dens_cur(:,:,:,i,1)
dens_cur(:,:,:,i,:) = 0.d0
enddo
seed_dcx = halo_iteration_dens/dcx_dens
n_halo=int(inputs%n_halo*seed_dcx,Int64)
if(seed_dcx.lt.0.01) exit iterations
enddo iterations
#ifdef _MPI
!! Combine Spectra
if(inputs%calc_halo.ge.1) then
call parallel_sum(spec%halo)
endif
#endif
deallocate(dens_cur)
end subroutine halo