halo Subroutine

subroutine halo
    !+ Calculates halo neutral density and spectra
    integer :: i,j,k !indices of cells
    integer(Int64) :: ihalo,n_halo !! counter
    real(Float64), dimension(3) :: ri    !! start position
    real(Float64), dimension(3) :: vihalo!! velocity bulk plasma ion
    integer,dimension(3) :: ind    !! 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 :: 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
    !! Halo iteration
    integer :: hh !! counters
    real(Float64) :: dcx_dens, halo_iteration_dens, seed_dcx
    real(Float64) :: r_nlaunchijk, local_dens
    integer :: s1type  ! halo iteration
    integer :: s2type  ! halo iteration

    s1type = fida_type
    s2type = brems_type

    dcx_dens = halo_iter_dens(halo_type)
    halo_iter_dens(s1type) = 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
    inv_ng = 100.0/real(beam_grid%ngrid)
    neut%dens(:,s1type,:,:,:) = neut%dens(:,halo_type,:,:,:)
    n_halo = inputs%n_halo
    seed_dcx = 1.0
    iterations: do hh=1,200
        papprox=0.d0
        papprox_tot=0.d0
        tot_denn=0.d0
        halo_iter_dens(s2type) = 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(:,s1type,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(n_halo,papprox,papprox_tot,nlaunch)

        if(inputs%verbose.ge.1) then
            write(*,'(T6,"# of markers: ",i9," --- Seed/DCX: ",f5.3)') int(sum(nlaunch),Int64), seed_dcx
        endif
        ccnt=0.d0
        !$OMP PARALLEL DO schedule(dynamic,1) collapse(3) private(i,j,k,ihalo,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
                    local_dens = 0.0d0
                    !! Loop over the markers
                    loop_over_halos: do ihalo=1,int(nlaunch(i,j,k),Int64)
                        !! Calculate ri,vhalo and track
                        ind = [i, j, k]
                        r_nlaunchijk = 1.0d0/nlaunch(i,j,k)
                        call mc_halo(ind,vihalo,ri)
                        call track(ri,vihalo,tracks,ncell)
                        if(ncell.eq.0)cycle loop_over_halos

                        !! Calculate CX probability
                        call get_beam_cx_rate(tracks(1)%ind,ri,vihalo,thermal_ion,[s1type],rates)
                        if(sum(rates).le.0.)cycle loop_over_halos

                        !! 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

                        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, s2type)
                    enddo loop_over_halos

                    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
        !$OMP END PARALLEL DO

        if(halo_iter_dens(s2type)/halo_iter_dens(s1type).gt.1.0) then
            if(inputs%verbose.ge.0) then
                write(*,'(a)') "HALO: Halo generation density exceeded seed density. This shouldn't happen."
            endif
            exit iterations
        endif

        halo_iteration_dens = halo_iter_dens(s2type)
        halo_iter_dens(s1type) = halo_iter_dens(s2type)

        neut%dens(:,halo_type,:,:,:)= neut%dens(:,halo_type,:,:,:) &
                                           + neut%dens(:,s2type,:,:,:)
        neut%dens(:,s1type,:,:,:)= neut%dens(:,s2type,:,:,:)
        neut%dens(:,s2type,:,:,:)= 0.

        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
    !! set the neutral density in s1type(fida_type) and s2type (brems) to 0!
    neut%dens(:,s1type,:,:,:) = 0.d0
    neut%dens(:,s2type,:,:,:) = 0.d0

end subroutine halo