store_neutrals Subroutine

subroutine store_neutrals(ind, neut_type, dens, vn, store_iter)
    !Store neutrals in [[libfida:neut]] at indices `ind`
    integer(Int32), dimension(3), intent(in) :: ind
        !+ [[libfida:beam_grid]] indices
    integer, intent(in)                      :: neut_type
        !+ Neutral type
    real(Float64), dimension(:), intent(in)  :: dens
        !+ Neutral density [neutrals/cm^3]
    real(Float64), dimension(3), intent(in)  :: vn
        !+ Neutral particle velocity [cm/s]
    logical, intent(in), optional            :: store_iter
        !+ Store DCX/Halo iteration density in [[libfida:halo_iter_dens]]
    logical :: iter

    integer :: i, j, k
    if(present(store_iter)) then
        iter = store_iter
    else
        iter = .False.
    endif

    i = ind(1) ; j = ind(2) ; k = ind(3)

    !$OMP CRITICAL(store_neutrals_1)
    if(iter) halo_iter_dens(neut_type) = halo_iter_dens(neut_type) + sum(dens)
    select case (neut_type)
        case (nbif_type)
            neut%full(:,i,j,k) = neut%full(:,i,j,k) + dens
        case (nbih_type)
            neut%half(:,i,j,k) = neut%half(:,i,j,k) + dens
        case (nbit_type)
            neut%third(:,i,j,k) = neut%third(:,i,j,k) + dens
        case (dcx_type)
            neut%dcx(:,i,j,k) = neut%dcx(:,i,j,k) + dens
        case (halo_type)
            neut%halo(:,i,j,k) = neut%halo(:,i,j,k) + dens
        case default
            if(inputs%verbose.ge.0) then
                write(*,'("STORE_NEUTRALS: Unknown neutral type: ",i2)') neut_type
            endif
            stop
    end select
    !$OMP END CRITICAL(store_neutrals_1)

end subroutine store_neutrals