Calculates Direct Charge Exchange (DCX) neutral density and spectra
subroutine dcx
!+ Calculates Direct Charge Exchange (DCX) neutral density and spectra
integer :: ic,i,j,k,ncell,is
integer(Int64) :: idcx !! counter
real(Float64), dimension(3) :: ri !! start position
real(Float64), dimension(3) :: vihalo
integer,dimension(3) :: ind
integer,dimension(3) :: neut_types = [1,2,3]
!! 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) :: fi_correction, dcx_dens
!! Initialized Neutral Population
call init_neutral_population(neut%dcx)
papprox=0.d0
tot_denn=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(neut%full%dens(:,i,j,k)) + &
sum(neut%half%dens(:,i,j,k)) + &
sum(neut%third%dens(:,i,j,k))
papprox(i,j,k)= tot_denn*sum(plasma%deni)
enddo
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(inputs%n_dcx,papprox,nlaunch)
if(inputs%verbose.ge.1) then
write(*,'(T6,"# of markers: ",i10)') sum(nlaunch)
endif
!$OMP PARALLEL DO schedule(dynamic,1) private(i,j,k,ic,is,idcx,ind,vihalo, &
!$OMP& ri,tracks,ntrack,rates,denn,states,jj,photons,plasma,fi_correction)
loop_over_cells: do ic = istart, ncell, istep
call ind2sub(beam_grid%dims,cell_ind(ic),ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
loop_over_species: do is=1, n_thermal !This loop has to come first
!! Loop over the markers
loop_over_dcx: do idcx=1, nlaunch(i,j,k)
!! Calculate ri,vhalo and track
call mc_beam_grid(ind, ri)
call get_plasma(plasma, pos=ri)
call mc_halo(plasma, thermal_mass(is), vihalo)
call track(ri,vihalo,tracks,ntrack)
if(ntrack.eq.0) cycle loop_over_dcx
!! Calculate CX probability
call get_total_cx_rate(tracks(1)%ind, ri, vihalo, neut_types, rates)
if(sum(rates).le.0.) cycle loop_over_dcx
!! Solve collisional radiative model along track
call get_plasma(plasma,pos=tracks(1)%pos)
!! Weight CX rates by ion source density
if(beam_mass.eq.thermal_mass(is)) then
states = rates*(plasma%deni(is) + plasma%denf)
if(sum(states).eq.0) cycle loop_over_dcx
fi_correction = max(plasma%deni(is)/(plasma%deni(is)+plasma%denf),0.d0)
else
states = rates*plasma%deni(is)
if(sum(states).eq.0) cycle loop_over_dcx
fi_correction = 1.d0
endif
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_mass(is),vihalo,tracks(jj)%time,states,denn,photons)
call store_neutrals(tracks(jj)%ind,tracks(jj)%pos,vihalo,dcx_type,denn/nlaunch(i,j,k))
photons = fi_correction*photons !! Correct for including fast-ions in states
if((photons.gt.0.d0).and.(inputs%calc_dcx.ge.1)) then
call store_photons(tracks(jj)%pos,vihalo, thermal_lambda0(is), photons/nlaunch(i,j,k),&
&spec%dcx(:,:,:,is),spec%dcxstokes(:,:,:,:,is))
endif
if((photons.gt.0.d0).and.(inputs%calc_res.ge.1)) then
call store_photon_birth(tracks(1)%pos, photons/nlaunch(i,j,k), spatres%dcx)
endif
enddo loop_along_track
enddo loop_over_dcx
enddo loop_over_species
enddo loop_over_cells
!$OMP END PARALLEL DO
#ifdef _MPI
!! Combine densities
call parallel_merge_populations(neut%dcx)
if(inputs%calc_dcx.ge.1) then
call parallel_sum(spec%dcx)
endif
if(inputs%calc_res.ge.1) then
do jj=1,spec_chords%nchan
call parallel_merge_reservoirs(spatres%dcx(jj))
enddo
endif
#endif
dcx_dens = sum(neut%dcx%dens)
if((dcx_dens.eq.0).or.isnan(dcx_dens)) then
write(*,*) 'DCX: DCX density is zero or nan: ', dcx_dens
stop
endif
end subroutine dcx