  subroutine nbi_bsource_write
    use nbi_bscom
    use nbi_com,only : workpath,runid,mibs2use,mibs_use,nsfast,mibs
    use nbi_com,only : nbs_step,nsbeam,nlfprod
    use nbi_mpi_com
    use netcdf
    use logmod
    implicit none
    character*150 :: suffix  ! birth source filenames
    character*150 :: fname
    character*30 :: cvar
    character*12 cstep
    integer :: icdf,i,ierr
    integer :: id_mibs_use,id_mc
    integer, dimension(1) :: idims
    integer,dimension(mibs_use) :: 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
    character*32 zlbls(mibs),zlbl  ! fast ion labels
    integer ilbls(mibs),ilbl
    integer :: inum,i_mc,ich_mc,id_ich_mc
    character*32  cmibs(mibs_use) !object's label

    zlbls=' '
    call nbfi_lbl(zlbls,mibs,inum)
    if(inum.ne.nsfast) then
       call errorlog(' nbi_bs_write: expected: inum = nsfast!')
       call errorlog(' nbi_bs_write: inum = ',inum)
       call errorlog(' nbi_bs_write: nsfast = ',nsfast)
       call bad_exit("nbi_bs_write : inconsistency in number of fast species!")
    endif
    do i=1,inum
       ilbls(i)=len_trim(zlbls(i))
    enddo
    i_mc=0
    do i=1,nsbeam
       if(nlfprod(i)) cycle
       if(mibs2use(i).gt.0) then
          zlbl=zlbls(i)
          ilbl=ilbls(i)
          if(index(zlbl,'_MCBEAM ').le.0) then
             cycle
          else
             i_mc=i_mc+1
             cmibs(i_mc)=zlbl(1:ilbl)
          endif
       endif
    enddo
    if(i_mc.eq.0) then
       call errorlog(' nbi_bs_write: cannot find mcbeam species')
       return
    endif
    write(cstep,'(i5)') nbs_step+1
    suffix='_birth_cpu'
    fname=' '
    call default_filename_ptcl(workpath,runid,fname,suffix,myidd,numprocd)
    fname=trim(fname)//'_'//trim(adjustl(cstep))
    ierr = nf90_create(path=trim(fname),cmode= or(nf90_clobber,nf90_64bit_offset),ncid=icdf)
    if(ierr.ne.0) then
       call errorlog(' nbi_bs_write: cannot create file'//trim(fname))
       call bad_exit("nbi_bs_write :  cannot create file")
    endif
    call check(nf90_def_dim(icdf,'mibs_use',mibs_use,id_mibs_use),'mibs_use')
    ich_mc=32 !dimension for character array
    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,i_mc
       cstep=trim(cmibs(i))
       cvar = 'nbsdep_'//trim(adjustl(cstep))   
       call check(nf90_def_dim(icdf,cvar,nbs_dep(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,i_mc
       cstep=trim(cmibs(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)
       cstep=trim(cmibs(i))
       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,cmibs(1:i_mc)),cvar)
    do i=1,i_mc
       cstep=trim(cmibs(i))
       cvar = 'bs_r_'//trim(adjustl(cstep))   
       call check( nf90_put_var(icdf,id_bs_r(i),bs_r(1:nbs_dep(i),i)),cvar)
       cvar = 'bs_z_'//trim(adjustl(cstep))   
       call check( nf90_put_var(icdf,id_bs_z(i),bs_z(1:nbs_dep(i),i)),cvar)
       cvar = 'bs_rgc_'//trim(adjustl(cstep))   
       call check( nf90_put_var(icdf,id_bs_rgc(i),bs_rgc(1:nbs_dep(i),i)),cvar)
       cvar = 'bs_zgc_'//trim(adjustl(cstep))   
       call check( nf90_put_var(icdf,id_bs_zgc(i),bs_zgc(1:nbs_dep(i),i)),cvar)
       cvar = 'bs_xksid_'//trim(adjustl(cstep))   
       call check( nf90_put_var(icdf,id_bs_xksid(i),bs_xksid(1:nbs_dep(i),i)),cvar)
       cvar = 'bs_einj_'//trim(adjustl(cstep))   
       call check( nf90_put_var(icdf,id_bs_einj(i),bs_einj(1:nbs_dep(i),i)),cvar)
       cvar = 'bs_wght_'//trim(adjustl(cstep))   
       call check( nf90_put_var(icdf,id_bs_wght(i),bs_wght(1:nbs_dep(i),i)),cvar)
       cvar = 'bs_zeta_'//trim(adjustl(cstep))   
       call check( nf90_put_var(icdf,id_bs_zeta(i),bs_zeta(1:nbs_dep(i),i)),cvar)
       cvar = 'bs_time_'//trim(adjustl(cstep))   
       call check( nf90_put_var(icdf,id_bs_time(i),bs_time(1:nbs_dep(i),i)),cvar)
       cvar = 'bs_ib_'//trim(adjustl(cstep))   
       call check( nf90_put_var(icdf,id_bs_ib(i),bs_ib(1:nbs_dep(i),i)),cvar)
    enddo
    ierr = nf90_close( icdf)
    if(ierr.ne.0) then
       call errorlog(' nbi_bs_write: netcdf close failure '//trim(fname))
       call bad_exit("nbi_bs_write :  cannot create file")
    endif
    call nbi_bs_dalloc
    
    
  contains 
    subroutine check(status,varnam)
      integer, intent ( in) :: status 
      character*(*),intent(in) :: varnam ! dimension name input 
      
      if(status /= nf90_noerr) then 
         call errorlog(' nbi_bs_write: error:'& 
              &// trim(nf90_strerror(status))//' for var.= '//trim(varnam))
         call exit(status) 
      end if
    end subroutine check
    
  end subroutine nbi_bsource_write
