ndmc Subroutine

public subroutine ndmc()

Calculates neutral beam deposition and spectra

Arguments

None

Calls

proc~~ndmc~~CallsGraph proc~ndmc ndmc proc~mc_nbi mc_nbi proc~ndmc->proc~mc_nbi interface~randind randind proc~ndmc->interface~randind interface~store_neutrals store_neutrals proc~ndmc->interface~store_neutrals proc~randu randu proc~ndmc->proc~randu proc~track track proc~ndmc->proc~track proc~mc_nbi->proc~randu proc~in_plasma in_plasma proc~mc_nbi->proc~in_plasma proc~grid_intersect grid_intersect proc~mc_nbi->proc~grid_intersect proc~randn randn proc~mc_nbi->proc~randn proc~store_neutrals_track store_neutrals_track interface~store_neutrals->proc~store_neutrals_track proc~store_neutrals_cell store_neutrals_cell interface~store_neutrals->proc~store_neutrals_cell omp_get_thread_num omp_get_thread_num proc~randu->omp_get_thread_num proc~rng_uniform rng_uniform proc~randu->proc~rng_uniform proc~get_indices get_indices proc~track->proc~get_indices proc~track->proc~in_plasma interface~interpol_coeff interpol_coeff proc~in_plasma->interface~interpol_coeff proc~xyz_to_uvw xyz_to_uvw proc~in_plasma->proc~xyz_to_uvw proc~randn->omp_get_thread_num proc~rng_normal rng_normal proc~randn->proc~rng_normal proc~interpol2d_coeff interpol2D_coeff interface~interpol_coeff->proc~interpol2d_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~interpol1d_coeff_arr interpol1D_coeff_arr interface~interpol_coeff->proc~interpol1d_coeff_arr proc~rng_normal->proc~rng_uniform proc~interpol2d_coeff_arr->proc~interpol2d_coeff proc~interpol1d_coeff_arr->proc~interpol1d_coeff

Called by

proc~~ndmc~~CalledByGraph proc~ndmc ndmc program~fidasim fidasim program~fidasim->proc~ndmc

Contents

Source Code


Source Code

subroutine ndmc
    !+ Calculates neutral beam deposition and spectra
    integer :: neut_type !! full half third energy
    real(Float64)  :: nlaunch   !! nr. of markers
    real(Float64)  :: nneutrals !! # NBI particles
    real(Float64), dimension(3) :: vnbi !! velocities(full..)
    real(Float64), dimension(3) :: rnbi !! initial position

    integer(Int64) :: jj, ii, kk, cnt
    integer :: ncell
    type(ParticleTrack), dimension(beam_grid%ntrack) :: tracks
    integer(Int64), dimension(3) :: nl_birth
    type(LocalProfiles) :: plasma
    real(Float64), dimension(nlevs) :: states, dens
    real(Float64) :: photons, iflux, r_nlaunch
    integer(Int32), dimension(3) :: ind
    real(Float64), dimension(1) :: randomu
    integer, dimension(1) :: randi
    logical :: err

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

    !! # of injected neutrals = NBI power/energy_per_particle
    nneutrals=1.d6*nbi%pinj/ (1.d3*nbi%einj*e0 &
         *( nbi%current_fractions(1)      &
         +  nbi%current_fractions(2)/2.d0 &
         +  nbi%current_fractions(3)/3.d0 ) )

    nlaunch=real(inputs%n_nbi)
    r_nlaunch = 1.d0/nlaunch
    cnt = 0
    !$OMP PARALLEL DO schedule(guided) &
    !$OMP& private(vnbi,rnbi,tracks,ncell,plasma,nl_birth,randi, &
    !$OMP& states,dens,iflux,photons,neut_type,jj,ii,kk,ind,err)
    loop_over_markers: do ii=1,inputs%n_nbi
        if(inputs%calc_birth.ge.1) then
            !! Determine the type of birth particle (1, 2, or 3)
            nl_birth = 0
            do kk=1,inputs%n_birth
                call randind(nbi%current_fractions,randi)
                nl_birth(randi(1)) = nl_birth(randi(1)) + 1
            enddo
        endif
        energy_fractions: do neut_type=1,3
            !! (type = 1: full energy, =2: half energy, =3: third energy
            call mc_nbi(vnbi,neut_type,rnbi,err)
            if(err) cycle energy_fractions

            call track(rnbi,vnbi,tracks,ncell)
            if(ncell.eq.0) cycle energy_fractions

            !! Solve collisional radiative model along track
            states=0.d0
            states(1)=nneutrals*nbi%current_fractions(neut_type)/beam_grid%dv
            loop_along_track: do jj=1,ncell
                iflux = sum(states)
                ind = tracks(jj)%ind
                call get_plasma(plasma,pos=tracks(jj)%pos)

                call colrad(plasma,beam_ion,vnbi,tracks(jj)%time,states,dens,photons)
                tracks(jj)%dens = dens*r_nlaunch
                tracks(jj)%flux = (iflux - sum(states))*beam_grid%dv*r_nlaunch

                if(inputs%calc_birth.ge.1) then
                    call store_births(ind,neut_type,tracks(jj)%flux)
                endif

                if((photons.gt.0.d0).and.(inputs%calc_bes.ge.1)) then
                    call store_bes_photons(tracks(jj)%pos,vnbi,photons*r_nlaunch,neut_type)
                endif
            enddo loop_along_track
            call store_neutrals(tracks, ncell, neut_type)

            if(inputs%calc_birth.ge.1) then
                !! Sample according to deposited flux along neutral trajectory
                !$OMP CRITICAL(ndmc_birth)
                do kk=1,nl_birth(neut_type)
                    call randind(tracks(1:ncell)%flux,randi)
                    call randu(randomu)
                    birth%neut_type(birth%cnt) = neut_type
                    birth%ind(:,birth%cnt) = tracks(randi(1))%ind
                    birth%vi(:,birth%cnt) = vnbi
                    birth%ri(:,birth%cnt) = tracks(randi(1))%pos + &
                                            vnbi*(tracks(randi(1))%time*(randomu(1)-0.5))
                    birth%cnt = birth%cnt+1
                enddo
                !$OMP END CRITICAL(ndmc_birth)
            endif
        enddo energy_fractions
        if (inputs%verbose.ge.2)then
            cnt = cnt + 1
            WRITE(*,'(f7.2,"% completed",a,$)') 100*cnt/nlaunch,char(13)
        endif
    enddo loop_over_markers
    !$OMP END PARALLEL DO

    if(nbi_outside.gt.0)then
        if(inputs%verbose.ge.0) then
            write(*,'(T4,a, f6.2)') 'Percent of markers outside the grid: ', &
                                  100.*nbi_outside/(3.*inputs%n_nbi)
        endif
        if(sum(neut%dens).eq.0) stop 'Beam does not intersect the grid!'
    endif

end subroutine ndmc