halo Subroutine

public subroutine halo()

Calculates halo neutral density and spectra

Arguments

None

Calls

proc~~halo~~CallsGraph proc~halo halo proc~mc_halo mc_halo proc~halo->proc~mc_halo proc~get_nlaunch get_nlaunch proc~halo->proc~get_nlaunch proc~ind2sub ind2sub proc~halo->proc~ind2sub omp_get_max_threads omp_get_max_threads proc~halo->omp_get_max_threads proc~get_plasma get_plasma proc~halo->proc~get_plasma proc~track track proc~halo->proc~track interface~parallel_sum parallel_sum proc~halo->interface~parallel_sum omp_get_thread_num omp_get_thread_num proc~halo->omp_get_thread_num proc~mc_halo->proc~get_plasma interface~randn randn proc~mc_halo->interface~randn interface~randu randu proc~mc_halo->interface~randu proc~rng_init rng_init proc~get_nlaunch->proc~rng_init interface~randind_cdf randind_cdf proc~get_nlaunch->interface~randind_cdf proc~cumsum cumsum proc~get_nlaunch->proc~cumsum proc~uvw_to_xyz uvw_to_xyz proc~get_plasma->proc~uvw_to_xyz proc~in_plasma in_plasma proc~get_plasma->proc~in_plasma proc~get_position get_position proc~get_plasma->proc~get_position proc~xyz_to_uvw xyz_to_uvw proc~get_plasma->proc~xyz_to_uvw proc~get_fields get_fields proc~track->proc~get_fields proc~doppler_stark doppler_stark proc~track->proc~doppler_stark proc~track->proc~in_plasma proc~get_indices get_indices proc~track->proc~get_indices proc~parallel_sum_i1 parallel_sum_i1 interface~parallel_sum->proc~parallel_sum_i1 proc~parallel_sum_i0 parallel_sum_i0 interface~parallel_sum->proc~parallel_sum_i0 proc~parallel_sum_d4 parallel_sum_d4 interface~parallel_sum->proc~parallel_sum_d4 proc~parallel_sum_d5 parallel_sum_d5 interface~parallel_sum->proc~parallel_sum_d5 proc~parallel_sum_d2 parallel_sum_d2 interface~parallel_sum->proc~parallel_sum_d2 proc~parallel_sum_d3 parallel_sum_d3 interface~parallel_sum->proc~parallel_sum_d3 proc~parallel_sum_d0 parallel_sum_d0 interface~parallel_sum->proc~parallel_sum_d0 proc~parallel_sum_d1 parallel_sum_d1 interface~parallel_sum->proc~parallel_sum_d1 proc~parallel_sum_i2 parallel_sum_i2 interface~parallel_sum->proc~parallel_sum_i2 proc~get_fields->proc~uvw_to_xyz proc~get_fields->proc~in_plasma proc~get_fields->proc~xyz_to_uvw proc~calc_perp_vectors calc_perp_vectors proc~get_fields->proc~calc_perp_vectors proc~rng_seed rng_seed proc~rng_init->proc~rng_seed proc~my_rank my_rank proc~rng_init->proc~my_rank mpi_allreduce mpi_allreduce proc~parallel_sum_i1->mpi_allreduce proc~parallel_sum_i0->mpi_allreduce proc~parallel_sum_d4->mpi_allreduce proc~in_plasma->proc~xyz_to_uvw interface~interpol_coeff interpol_coeff proc~in_plasma->interface~interpol_coeff proc~cyl_to_uvw cyl_to_uvw proc~in_plasma->proc~cyl_to_uvw proc~cyl_to_xyz cyl_to_xyz proc~get_position->proc~cyl_to_xyz proc~parallel_sum_d5->mpi_allreduce proc~parallel_sum_d2->mpi_allreduce proc~parallel_sum_d3->mpi_allreduce proc~parallel_sum_d0->mpi_allreduce proc~parallel_sum_d1->mpi_allreduce proc~parallel_sum_i2->mpi_allreduce proc~cyl_interpol3d_coeff cyl_interpol3D_coeff interface~interpol_coeff->proc~cyl_interpol3d_coeff proc~interpol1d_coeff interpol1D_coeff interface~interpol_coeff->proc~interpol1d_coeff proc~interpol2d_coeff_arr interpol2D_coeff_arr interface~interpol_coeff->proc~interpol2d_coeff_arr proc~interpol2d_coeff interpol2D_coeff interface~interpol_coeff->proc~interpol2d_coeff proc~cyl_interpol3d_coeff_arr cyl_interpol3D_coeff_arr interface~interpol_coeff->proc~cyl_interpol3d_coeff_arr proc~interpol1d_coeff_arr interpol1D_coeff_arr interface~interpol_coeff->proc~interpol1d_coeff_arr proc~cyl_to_xyz->proc~uvw_to_xyz proc~cyl_to_xyz->proc~cyl_to_uvw proc~cyl_interpol3d_coeff->proc~interpol2d_coeff proc~interpol2d_coeff_arr->proc~interpol2d_coeff proc~cyl_interpol3d_coeff_arr->proc~cyl_interpol3d_coeff proc~cyl_interpol3d_coeff_arr->proc~interpol2d_coeff proc~interpol1d_coeff_arr->proc~interpol1d_coeff

Called by

proc~~halo~~CalledByGraph proc~halo halo program~fidasim fidasim program~fidasim->proc~halo

Contents

Source Code


Source Code

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