subroutine nbi_bsource_conct(nfile,icd_f,fname_av,ier_c)
  !files with id (icd_f) will be concatenated and written to fname_av
  use iso_c_binding, only: fp => c_double
  use nbi_bscom
  use logmod
  use netcdf

  implicit none
  integer, intent(in) :: nfile !# of files to concatenate
  integer, intent(in) :: icd_f(nfile) ! file's ids (files have been opened in nbi_avg_bsource.f90)
  integer, intent(out) :: ier_c
  character*(*), intent(in) :: fname_av ! file name to where write data
  integer :: ii,i, ierr, icdf,ik,j
  character*30 :: cvar
  integer :: i_dim
  integer :: mibs_use
  character*32,dimension(:,:), allocatable :: cmibs
  character*32,dimension(20) :: cl ! labels used in files
  character*32 :: cstep
  logical,dimension(0:20) :: nlcl
  integer, dimension(:,:), allocatable :: nbsdep ! # of deposited events for a given species and file
  integer, dimension(:), allocatable :: nbsdep_s ! # of deposited events for a given species
  integer :: n_bsdep
  real(fp), dimension(:), allocatable :: tmp
  integer, dimension(:), allocatable :: itmp
  integer :: id_mibs_use,id_mc,i_in,i_lst
  integer, dimension(1) :: idims
  integer, dimension(:), allocatable :: id_nbsdep,id_bs_r,id_bs_z,id_bs_xksid,id_bs_einj,&
       &id_bs_wght,id_bs_zeta,id_bs_time,id_bs_ib,id_bs_rgc,id_bs_zgc
  integer :: ich_mc,id_ich_mc

  ier_c=0 
  ierr = nf90_create(path=trim(fname_av),cmode= or(nf90_clobber,nf90_64bit_offset),ncid=icdf)
  if(ierr.ne.0) then
     call errorlog(' nbi_bsource_conct: cannot create file'//trim(fname_av))
     do ii=1,nfile
        !close  opened files before exit
        ierr=nf90_close(icd_f(ii))
     enddo
     ierr = nf90_close( icdf)
     call bad_exit("nbi_avg_bsource :  cannot create file_av")
  endif
  
  mibs_use=0
  ! find max # of  species into all files
  do i=1,nfile
     cvar = 'mibs_use'
     i_dim=0
     call check_inq_dim(icd_f(i),trim(cvar),i_dim,ierr)
     if(ierr.ne.0.or.i_dim.eq.0) then
        !file is empty
        cycle
     else
       mibs_use =max(i_dim,mibs_use)
     endif
  enddo
  if(mibs_use.le.0) then
     ier_c=ier_c+1
     call errorlog(' nbi_bsource_conct: all files are empty')
     ierr = nf90_close( icdf)
     return
  endif


  !find species label

  if(allocated(cmibs)) deallocate(cmibs)
  allocate(cmibs(nfile,mibs_use))
  cmibs=' '
  
  do i=1,nfile
     !read labels from all files
     cvar = 'mclabel'
     i_dim=0
     call check_inq(icd_f(i),trim(cvar),i_dim,ierr) !get var. id
     call check( nf90_get_var(icd_f(i),i_dim,cmibs(i,:)),cvar)
  enddo

  !find all diff. labels

  cl=' ' ! not supposed to be ' '
  ii=0
  nlcl(0)=.true.
  do i=1,nfile
     do j=1,mibs_use
        ik=0
        do while(ik.le.ii)
           ik=ik+1
           nlcl(ik)=nlcl(ik-1).and.(cl(ik).ne.cmibs(i,j))
           if( nlcl(ik)) then
              if(ik.eq.ii+1) then
                 ii=ii+1
                 cl(ik)=cmibs(i,j)
                 exit
              endif
           else
              exit
           endif
        enddo
     enddo
  enddo
  if(allocated(cmibs)) deallocate(cmibs)

  if(ii.le.0) then
     call errorlog(' nbi_bsource_conct: inconsistency in species found')
     ier_c=ier_c+1
     return
  endif
  if(allocated(nbsdep)) deallocate(nbsdep)
  allocate(nbsdep(nfile,ii))
  nbsdep=0
  if(allocated(nbsdep_s)) deallocate(nbsdep_s)
  allocate(nbsdep_s(ii))
  nbsdep_s=0

  do j=1,ii
     do i=1,nfile
        cvar='nbsdep_'//trim(adjustl(cl(j)))
        i_dim=0
        call check_inq_dim(icd_f(i),trim(cvar),i_dim,ierr)
        if(ierr.ne.0.or.i_dim.eq.0) then
           !file is empty
           cycle
        else
           nbsdep(i,j)=i_dim
        endif
     enddo
  enddo
  nbsdep_s(1:ii)=sum(nbsdep(:,1:ii),dim=1)
  n_bsdep=maxval(nbsdep_s)
  if(n_bsdep.le.0) then
     call warnlog(' nbi_bsource_conct: zero events occured')
     ierr = nf90_close( icdf)
     ier_c = ier_c+1
     return
  endif

  if(allocated(tmp)) deallocate(tmp)
  allocate(tmp(n_bsdep))
  tmp=0.0_fp
  if(allocated(itmp)) deallocate(itmp)
  allocate(itmp(n_bsdep))
  itmp=0

  if(allocated(id_nbsdep)) deallocate(id_nbsdep)
  allocate(id_nbsdep(ii))

  if(allocated(id_bs_r)) deallocate(id_bs_r)
  allocate(id_bs_r(ii))

  if(allocated(id_bs_z)) deallocate(id_bs_z)
  allocate(id_bs_z(ii))

  if(allocated(id_bs_rgc)) deallocate(id_bs_rgc)
  allocate(id_bs_rgc(ii))

  if(allocated(id_bs_zgc)) deallocate(id_bs_zgc)
  allocate(id_bs_zgc(ii))

  if(allocated(id_bs_xksid)) deallocate(id_bs_xksid)
  allocate(id_bs_xksid(ii))

  if(allocated(id_bs_einj)) deallocate(id_bs_einj)
  allocate(id_bs_einj(ii))

  if(allocated(id_bs_wght)) deallocate(id_bs_wght)
  allocate(id_bs_wght(ii))

  if(allocated(id_bs_zeta)) deallocate(id_bs_zeta)
  allocate(id_bs_zeta(ii))

  if(allocated(id_bs_time)) deallocate(id_bs_time)
  allocate(id_bs_time(ii))

  if(allocated(id_bs_ib)) deallocate(id_bs_ib)
  allocate(id_bs_ib(ii))

  !all are ready, get data
  ! define mode
  call check(nf90_def_dim(icdf,'mibs_use',ii,id_mibs_use),'mibs_use')
  ich_mc=32
  call check(nf90_def_dim(icdf,'ich_mc',ich_mc,id_ich_mc),'ich_mc')
  cvar = 'mclabel'   
  call check(nf90_def_var(icdf,cvar,nf90_char,(/id_ich_mc,id_mibs_use/),id_mc),cvar)
  do i=1,ii
     cstep=trim(adjustl(cl(i)))
     cvar = 'nbsdep_'//trim(adjustl(cstep))   
     call check(nf90_def_dim(icdf,cvar,nbsdep_s(i),id_nbsdep(i)),cvar)
     idims = (/id_nbsdep(i)/)
     cvar = 'bs_r_'//trim(adjustl(cstep))   
     call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_bs_r(i)),cvar)
     cvar = 'bs_z_'//trim(adjustl(cstep))   
     call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_bs_z(i)),cvar)
     cvar = 'bs_rgc_'//trim(adjustl(cstep))   
     call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_bs_rgc(i)),cvar)
     cvar = 'bs_zgc_'//trim(adjustl(cstep))   
     call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_bs_zgc(i)),cvar)
     cvar = 'bs_xksid_'//trim(adjustl(cstep))   
     call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_bs_xksid(i)),cvar)
     cvar = 'bs_einj_'//trim(adjustl(cstep))   
     call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_bs_einj(i)),cvar)
     cvar = 'bs_wght_'//trim(adjustl(cstep))   
     call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_bs_wght(i)),cvar)
     cvar = 'bs_zeta_'//trim(adjustl(cstep))   
     call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_bs_zeta(i)),cvar)
     cvar = 'bs_time_'//trim(adjustl(cstep))   
     call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_bs_time(i)),cvar)
     cvar = 'bs_ib_'//trim(adjustl(cstep))   
     call check(nf90_def_var(icdf,cvar,nf90_int,idims,id_bs_ib(i)),cvar)
  enddo
  do i=1,ii
     cstep=trim(adjustl(cl(i)))
     cvar = 'bs_r_'//trim(adjustl(cstep))   
     call check(nf90_put_att(icdf,id_bs_r(i),cvar,"r at deposition, cm"),cvar)
     cvar = 'bs_z_'//trim(adjustl(cstep))   
     call check(nf90_put_att(icdf,id_bs_z(i),cvar,"z at deposition, cm"),cvar)
     cvar = 'bs_rgc_'//trim(adjustl(cstep))   
     call check(nf90_put_att(icdf,id_bs_rgc(i),cvar,"r at gc, cm"),cvar)
     cvar = 'bs_zgc_'//trim(adjustl(cstep))   
     call check(nf90_put_att(icdf,id_bs_zgc(i),cvar,"z at gc, cm"),cvar)
     cvar = 'bs_xksid_'//trim(adjustl(cstep))   
     call check(nf90_put_att(icdf,id_bs_xksid(i),cvar,"pitch angle, v_ii/v"),cvar)
     cvar = 'bs_einj_'//trim(adjustl(cstep))   
     call check(nf90_put_att(icdf,id_bs_einj(i),cvar,"energy, ev"),cvar)
     cvar = 'bs_wght_'//trim(adjustl(cstep))   
     call check(nf90_put_att(icdf,id_bs_wght(i),cvar,"weight deposited ptcl"),cvar)
     cvar = 'bs_zeta_'//trim(adjustl(cstep))   
     call check(nf90_put_att(icdf,id_bs_zeta(i),cvar,"toroidal angle at deposition, deg"),cvar)
     cvar = 'bs_time_'//trim(adjustl(cstep))   
     call check(nf90_put_att(icdf,id_bs_time(i),cvar,"time at deposition, sec"),cvar)
     cvar = 'bs_ib_'//trim(adjustl(cstep))   
     call check(nf90_put_att(icdf,id_bs_ib(i),cvar,"beam id"),cvar)
  enddo
  ierr = nf90_enddef(icdf) ! out of def mode and go to data mode
  
  cvar = 'mclabel'   
  call check( nf90_put_var(icdf,id_mc,cl(1:ii)),cvar)

  do j=1,ii ! # of species
     i_lst=0
     i_in=1
     cvar = 'bs_r_'//trim(adjustl(cl(j)))   
     do i=1,nfile !# of files
        i_in=i_lst+1
        i_lst=i_lst+nbsdep(i,j)
        i_dim=0
        call check_inq(icd_f(i),trim(cvar),i_dim,ierr) !get var. id
        call check( nf90_get_var(icd_f(i),i_dim,tmp(i_in:i_lst)),cvar)
     enddo
     call check( nf90_put_var(icdf,id_bs_r(j),tmp(1:nbsdep_s(j))),cvar)
     tmp=0.0_fp
     
     i_lst=0
     i_in=1
     cvar = 'bs_z_'//trim(adjustl(cl(j)))   
     do i=1,nfile !# of files
        i_in=i_lst+1
        i_lst=i_lst+nbsdep(i,j)
        i_dim=0
        call check_inq(icd_f(i),trim(cvar),i_dim,ierr) !get var. id
        call check( nf90_get_var(icd_f(i),i_dim,tmp(i_in:i_lst)),cvar)
     enddo
     call check( nf90_put_var(icdf,id_bs_z(j),tmp(1:nbsdep_s(j))),cvar)
     tmp=0.0_fp
     
     i_lst=0
     i_in=1
     cvar = 'bs_rgc_'//trim(adjustl(cl(j)))   
     do i=1,nfile !# of files
        i_in=i_lst+1
        i_lst=i_lst+nbsdep(i,j)
        i_dim=0
        call check_inq(icd_f(i),trim(cvar),i_dim,ierr) !get var. id
        call check( nf90_get_var(icd_f(i),i_dim,tmp(i_in:i_lst)),cvar)
     enddo
     call check( nf90_put_var(icdf,id_bs_rgc(j),tmp(1:nbsdep_s(j))),cvar)
     tmp=0.0_fp
     
     i_lst=0
     i_in=1
     cvar = 'bs_zgc_'//trim(adjustl(cl(j)))   
     do i=1,nfile !# of files
        i_in=i_lst+1
        i_lst=i_lst+nbsdep(i,j)
        i_dim=0
        call check_inq(icd_f(i),trim(cvar),i_dim,ierr) !get var. id
        call check( nf90_get_var(icd_f(i),i_dim,tmp(i_in:i_lst)),cvar)
     enddo
     call check( nf90_put_var(icdf,id_bs_zgc(j),tmp(1:nbsdep_s(j))),cvar)
     tmp=0.0_fp
     
     i_lst=0
     i_in=1
     cvar = 'bs_xksid_'//trim(adjustl(cl(j)))   
     do i=1,nfile !# of files
        i_in=i_lst+1
        i_lst=i_lst+nbsdep(i,j)
        i_dim=0
        call check_inq(icd_f(i),trim(cvar),i_dim,ierr) !get var. id
        call check( nf90_get_var(icd_f(i),i_dim,tmp(i_in:i_lst)),cvar)
     enddo
     call check( nf90_put_var(icdf,id_bs_xksid(j),tmp(1:nbsdep_s(j))),cvar)
     tmp=0.0_fp
     
     i_lst=0
     i_in=1
     cvar = 'bs_einj_'//trim(adjustl(cl(j)))   
     do i=1,nfile !# of files
        i_in=i_lst+1
        i_lst=i_lst+nbsdep(i,j)
        i_dim=0
        call check_inq(icd_f(i),trim(cvar),i_dim,ierr) !get var. id
        call check( nf90_get_var(icd_f(i),i_dim,tmp(i_in:i_lst)),cvar)
     enddo
     call check( nf90_put_var(icdf,id_bs_einj(j),tmp(1:nbsdep_s(j))),cvar)
     tmp=0.0_fp
     
     i_lst=0
     i_in=1
     cvar = 'bs_wght_'//trim(adjustl(cl(j)))   
     do i=1,nfile !# of files
        i_in=i_lst+1
        i_lst=i_lst+nbsdep(i,j)
        i_dim=0
        call check_inq(icd_f(i),trim(cvar),i_dim,ierr) !get var. id
        call check( nf90_get_var(icd_f(i),i_dim,tmp(i_in:i_lst)),cvar)
     enddo
     call check( nf90_put_var(icdf,id_bs_wght(j),tmp(1:nbsdep_s(j))),cvar)
     tmp=0.0_fp
     
     i_lst=0
     i_in=1
     cvar = 'bs_zeta_'//trim(adjustl(cl(j)))   
     do i=1,nfile !# of files
        i_in=i_lst+1
        i_lst=i_lst+nbsdep(i,j)
        i_dim=0
        call check_inq(icd_f(i),trim(cvar),i_dim,ierr) !get var. id
        call check( nf90_get_var(icd_f(i),i_dim,tmp(i_in:i_lst)),cvar)
     enddo
     call check( nf90_put_var(icdf,id_bs_zeta(j),tmp(1:nbsdep_s(j))),cvar)
     tmp=0.0_fp
     
     i_lst=0
     i_in=1
     cvar = 'bs_time_'//trim(adjustl(cl(j)))   
     do i=1,nfile !# of files
        i_in=i_lst+1
        i_lst=i_lst+nbsdep(i,j)
        i_dim=0
        call check_inq(icd_f(i),trim(cvar),i_dim,ierr) !get var. id
        call check( nf90_get_var(icd_f(i),i_dim,tmp(i_in:i_lst)),cvar)
     enddo
     call check( nf90_put_var(icdf,id_bs_time(j),tmp(1:nbsdep_s(j))),cvar)
     tmp=0.0_fp
     
     i_lst=0
     i_in=1
     cvar = 'bs_ib_'//trim(adjustl(cl(j)))   
     do i=1,nfile !# of files
        i_in=i_lst+1
        i_lst=i_lst+nbsdep(i,j)
        i_dim=0
        call check_inq(icd_f(i),trim(cvar),i_dim,ierr) !get var. id
        call check( nf90_get_var(icd_f(i),i_dim,itmp(i_in:i_lst)),cvar)
     enddo
     call check( nf90_put_var(icdf,id_bs_ib(j),itmp(1:nbsdep_s(j))),cvar)
     itmp=0.0_fp
  enddo
  if(allocated(tmp)) deallocate(tmp)
  if(allocated(itmp)) deallocate(itmp)
  if(allocated(id_nbsdep)) deallocate(id_nbsdep)
  if(allocated(id_bs_r)) deallocate(id_bs_r)
  if(allocated(id_bs_z)) deallocate(id_bs_z)
  if(allocated(id_bs_rgc)) deallocate(id_bs_rgc)
  if(allocated(id_bs_zgc)) deallocate(id_bs_zgc)
  if(allocated(id_bs_xksid)) deallocate(id_bs_xksid)
  if(allocated(id_bs_einj)) deallocate(id_bs_einj)
  if(allocated(id_bs_wght)) deallocate(id_bs_wght)
  if(allocated(id_bs_zeta)) deallocate(id_bs_zeta)
  if(allocated(id_bs_time)) deallocate(id_bs_time)
  if(allocated(id_bs_ib)) deallocate(id_bs_ib)
  
  ierr = nf90_close( icdf)
  if(ierr.ne.0) then
     call errorlog(' nbi_bsource_conct: cannot close file'//trim(fname_av))
     call bad_exit(" nbi_bsource_conct ")
  endif
  
contains 
  subroutine check(status,varnam)
    integer, intent ( in) :: status 
    character*(*),intent(in) :: varnam ! dimension name input 
 
    if(status /= nf90_noerr) then 
       call warnlog(' nbi_bsource_conct: '& 
            &// trim(nf90_strerror(status))//' for var.= '//trim(varnam) )
       return
    end if 
  end subroutine check 
  subroutine check_inq(i_cdf,varnam,id,ierrcdf)
    integer,intent(in) :: i_cdf
    integer,intent(out) :: id,ierrcdf
    character*(*),intent(in) :: varnam 
    ierrcdf = nf90_inq_varid(i_cdf,trim(varnam),id) 
    if(ierrcdf.ne.0) then 
       call warnlog(' nbi_bsource_conct: no such var='//trim(varnam))
       ierrcdf=1 
       return
    endif 
  end subroutine check_inq 
 
  subroutine check_inq_dim(i_cdf,varnam,idim,ierrcdf) 
    integer,intent(in) :: i_cdf
    integer,intent(out) :: ierrcdf 
    character*(*),intent(in) :: varnam ! dimension name input 
    integer, intent(out) :: idim        ! dimension size returned 
    integer ::id !var id 
    ierrcdf = nf90_inq_dimid(i_cdf,trim(varnam),id) 
    if(ierrcdf.ne.0) then 
       call warnlog(' nbi_bsource_conct: no such dim='//trim(varnam))
       ierrcdf=1 
       return
    endif 
    ierrcdf = nf90_inquire_dimension(i_cdf,id,len=idim) 
    if(ierrcdf.ne.0) then 
       call errorlog(' nbi_bsource_conct: error '//nf90_strerror(ierrcdf)//'for var.'//trim(varnam))
       call exit(ierrcdf) 
    endif 
  end subroutine check_inq_dim

end subroutine nbi_bsource_conct
