dcx Subroutine

subroutine dcx
    !+ Calculates Direct Charge Exchange (DCX) neutral density and spectra
    integer :: i,j,k !indices of cells
    integer(Int64) :: idcx !! counter
    real(Float64), dimension(3) :: ri    !! start position
    real(Float64), dimension(3) :: vihalo
    integer,dimension(3) :: ind    !! actual cell
    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 :: ncell
    type(ParticleTrack), dimension(beam_grid%ntrack) :: tracks  !! Particle tracks
    integer :: jj       !! counter along track
    real(Float64):: tot_denn, photons  !! photon flux
    real(Float64), dimension(beam_grid%nx,beam_grid%ny,beam_grid%nz) :: papprox, nlaunch !! approx. density
    real(Float64) :: papprox_tot, ccnt, inv_ng, r_nlaunchijk

    halo_iter_dens(halo_type) = 0.d0
    papprox=0.d0
    papprox_tot=0.d0
    tot_denn=0.d0
    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)
                tot_denn = sum(neut%dens(:,nbif_type,i,j,k)) + &
                           sum(neut%dens(:,nbih_type,i,j,k)) + &
                           sum(neut%dens(:,nbit_type,i,j,k))
                papprox(i,j,k)= tot_denn*(plasma%denp-plasma%denf)
                if(plasma%in_plasma) papprox_tot=papprox_tot+papprox(i,j,k)
            enddo
        enddo
    enddo

    call get_nlaunch(inputs%n_dcx,papprox,papprox_tot,nlaunch)

    if(inputs%verbose.ge.1) then
       write(*,'(T6,"# of markers: ",i9)') int(sum(nlaunch),Int64)
    endif

    ccnt=0.d0
    inv_ng = 100.0/real(beam_grid%ngrid)
    !$OMP PARALLEL DO collapse(3) schedule(dynamic,1) private(k,j,i,idcx,ind,vihalo, &
    !$OMP& ri,tracks,ncell,rates,denn,states,jj,photons,plasma,r_nlaunchijk)
    loop_along_z: do k = 1, beam_grid%nz
        loop_along_y: do j = 1, beam_grid%ny
            loop_along_x: do i = 1, beam_grid%nx
                !! Loop over the markers
                loop_over_dcx: do idcx=1,int(nlaunch(i,j,k),Int64)
                    !! Calculate ri,vhalo and track
                    ind = [i, j, k]
                    r_nlaunchijk = 1.d0/nlaunch(i,j,k)
                    call mc_halo(ind,vihalo,ri)
                    call track(ri,vihalo,tracks,ncell)
                    if(ncell.eq.0) cycle loop_over_dcx

                    !! Calculate CX probability
                    call get_beam_cx_rate(tracks(1)%ind,ri,vihalo,thermal_ion,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
                    states = rates*(plasma%denp - plasma%denf)

                    loop_along_track: do jj=1,ncell
                        call get_plasma(plasma,pos=tracks(jj)%pos)

                        call colrad(plasma,thermal_ion,vihalo,tracks(jj)%time,states,denn,photons)
                        tracks(jj)%dens = denn*r_nlaunchijk

                        if((photons.gt.0.d0).and.(inputs%calc_bes.ge.1)) then
                          call store_bes_photons(tracks(jj)%pos,vihalo,photons*r_nlaunchijk,halo_type)
                        endif

                    enddo loop_along_track
                    call store_neutrals(tracks, ncell, halo_type)
                enddo loop_over_dcx
                ccnt=ccnt+1
                if (inputs%verbose.ge.2)then
                    WRITE(*,'(f7.2,"% completed",a,$)') ccnt*inv_ng,char(13)
                endif
            enddo loop_along_x
        enddo loop_along_y
    enddo loop_along_z

end subroutine dcx