!this version has mpi_reduce and mpi_gatherv for routines
!used by orball
!*********************************************************************
  subroutine nubeam_buf_reduce 
!
      use nbi_com
      use nbi_mpi_com
      use nbi_pllbuf
      use logmod
      implicit none

      integer :: mpierr,jk
      character(len=ncfnm)   :: fname

      integer :: inum,i,j,k,i_wrk,ii, i_cx_do
      real(fp) :: r_wrk,rwrk(1)
      real(fp), dimension(:,:,:), allocatable :: zfbm

      call nbi_pput_orb
#ifdef __MPI
      mpierr=0
      if(myidd.eq.0) then
         if(.not.allocated(r8buf1)) then
            allocate(r8buf1(n_data))
            r8buf1=0.0_fp
         endif
      endif
      call infolog('nubeam_buf_reduce: nubeam_reduce_chunks_r8 start')
      call mpi_reduce(r8buf,r8buf1,n_data,MPI_REAL8,MPI_SUM,0,nbi_comm,mpierr)
!      if(myidd.eq.0) then
!         call nubeam_reduce_chunks_r8(r8buf,r8buf1,n_data,'MPI_SUM',20)
!      else
!         rwrk=0.0_fp
! myid > 0  does not need memory for r8buf1
!         call nubeam_reduce1_chunks_r8(r8buf,rwrk,n_data,'MPI_SUM',20)
!      endif
      call infolog('nubeam_buf_reduce: nubeam_reduce_chunks_r8 done')

!  write out orball buffer

      if (debug_mpi.and.(myidd.eq.0).and.(numprocd.gt.1)) then

       write(fname,'(a,a13,i1,a4)') trim(workpath)//trim(runid),'-orb-buf-out-',lsbeam,'.tst'
       call nbi_buf_netcdf('W',trim(adjustl(fname)),n_data,r8buf1)
      endif

      if (myidd.eq.0) then
         call nbi_pget_orb !from r8buf1
      endif
      nmax_data=0
      n_data=0
      if(allocated(r8buf)) deallocate(r8buf)
      if(allocated(r8buf1)) deallocate(r8buf1)

! reduce # of splittings
      mpierr=0
      i_wrk=0
      call  mpi_reduce(j_spl,i_wrk,1,MPI_INTEGER,  &
           MPI_SUM,0,nbi_comm,mpierr)
      if(myidd.eq.0) then
         j_spl=i_wrk
      endif


!------------------------------------------------------------------
!!! also reduce fbm and fbm_ptcl -- directly, not through buffer

      if(allocated(zfbm)) deallocate(zfbm)
      allocate(zfbm(mimxbe,mimxba,mimxbzf))
      zfbm=0.0_fp
      inum = mimxbe*mimxba*mimxbzf

      mpierr=0
      call infolog('nubeam_buf_reduce: nubeam_reduce_chunks_r8 start fbm')
      call mpi_reduce(fbm(1,1,1,lsbeam),zfbm,inum,MPI_REAL8,MPI_SUM,0,nbi_comm,mpierr)

!      if(myidd.eq.0) then
!         call nubeam_reduce_chunks_r8(fbm(1,1,1,lsbeam),zfbm,inum,'MPI_SUM',21)
!      else
!         rwrk=0.0_fp
!         call nubeam_reduce1_chunks_r8(fbm(1,1,1,lsbeam),rwrk,inum,'MPI_SUM',21)
!      endif
      call infolog('nubeam_buf_reduce: nubeam_reduce_chunks_r8 done fbm')

      if(myidd.eq.0) then
         fbm(1:mimxbe,1:mimxba,1:mimxbzf,lsbeam)=zfbm(1:mimxbe,1:mimxba,1:mimxbzf) !this is faster
!         do k=1,mimxbzf
!            do j=1,mimxba
!               do i=1,mimxbe
!                  fbm(i,j,k,lsbeam)=zfbm(i,j,k)
!               enddo
!            enddo
!         enddo
      endif
      if(nlfbmflr) then
         mpierr=0
         inum = mimxbe*mimxba*mimxbzf
         call infolog('nubeam_buf_reduce: nubeam_allreduce_chunks_r8 start fbm_ptcl')
!mg fbm_ptcl need to be allreduced to a 3d halo model
!         if(.not.allocated(zfbm)) allocate(zfbm(mimxbe,mimxba,mimxbzf))
         do ii=1,mgyro
            zfbm=0.0_fp
            call nubeam_allreduce_chunks_r8(fbm_ptcl(1,1,1,lsbeam,ii),zfbm,inum,'MPI_SUM',22)
            fbm_ptcl(1:mimxbe,1:mimxba,1:mimxbzf,lsbeam,ii)=zfbm(1:mimxbe,1:mimxba,1:mimxbzf) !this is faster
!            do k=1,mimxbzf
!               do j=1,mimxba
!                  do i=1,mimxbe
!                     fbm_ptcl(i,j,k,lsbeam,ii)=zfbm(i,j,k)
!                  enddo
!               enddo
!            enddo
         enddo
         call infolog('nubeam_buf_reduce: nubeam_allreduce_chunks_r8 done fbm_ptcl')
      endif !nlfbmflr
      
!------------------------------------------------------------------
!!! also reduce fbm_cx -- directly, not through buffer 
      if(allocated(fbm_cx)) then
         if(allocated(zfbm))deallocate(zfbm)
         allocate(zfbm(mimxbe,mimxba,mimxbz)) !all cpu, includes sol
         inum = mimxbe*mimxba*mimxbz
         zfbm=0.0_fp
         mpierr=0
         call mpi_allreduce(fbm_cx(1,1,1,lsbeam),zfbm,inum,MPI_REAL8,MPI_SUM,nbi_comm,mpierr)
         call infolog('nubeam_buf_reduce: nubeam_reduce_chunks_r8 start fbm_cx')
!         call nubeam_allreduce_chunks_r8(fbm_cx(1,1,1,lsbeam),zfbm,inum,'MPI_SUM',23)
!         call infolog('nubeam_buf_reduce: nubeam_reduce_chunks_r8 done fbm_cx')
         fbm_cx(1:mimxbe,1:mimxba,1:mimxbz,lsbeam)=zfbm(1:mimxbe,1:mimxba,1:mimxbz) 
!         do k=1,mimxbz
!            do j=1,mimxba
!               do i=1,mimxbe
!                  fbm_cx(i,j,k,lsbeam)=zfbm(i,j,k)
!               enddo
!            enddo
!         enddo
         
      endif

      if(allocated(zfbm)) deallocate(zfbm)

! reduce  # cx tracks from orbiting  from all cpu
      mpierr=0
      r_wrk=0 
      i_wrk=0
      if(xncx0_orb(lsbeam).gt.0) then
         i_cx_do=1
      else
         i_cx_do=0
      endif
      call  mpi_allreduce(i_cx_do,i_wrk,1,MPI_INTEGER,MPI_SUM, nbi_comm,mpierr)
      if(i_wrk.lt.numprocd) then  !skip improving cx stat, not all cpu has  cx tracks
         xncx0_orb(lsbeam)=0.0_fp
      endif
      call mpi_allreduce(xncx0_orb(lsbeam),r_wrk,1,MPI_REAL8,MPI_SUM,nbi_comm,mpierr)
      xncx0_orb(lsbeam)=r_wrk !now total # cx tracks from orbiting on all cpu
#else
#endif


      call infolog('nubeam_buf_reduce: nubeam_buf_reduce done')
      

      end subroutine nubeam_buf_reduce 
!---------------------------------------------------------------------
!*********************************************************************
  subroutine nubeam_cxbuf_reduce 
!
      use nbi_com
      use nbi_mpi_com
      use nbi_pllbuf
      use logmod
      implicit none

      integer :: mpierr,jk
      character(len=ncfnm)   :: fname

      integer :: inum,i,j,k,i_wrk
      real(fp), dimension(:,:,:), allocatable :: zfbm
      real(fp) :: rwrk(1)

      call nbi_pput_ntr


#ifdef __MPI
      mpierr=0
      if(myidd.eq.0) then
         if(.not.allocated(r8buf1)) then
            allocate(r8buf1(n_data))
            r8buf1=0.0_fp
         endif
      endif
      call infolog('nubeam_cxbuf_reduce: nubeam_reduce_chunks_r8 start')
      if(myidd.eq.0) then
         call nubeam_reduce_chunks_r8(r8buf,r8buf1,n_data,'MPI_SUM',20)
      else
         rwrk=0.0_fp
! myid > 0  does not need memory for r8buf1
         call nubeam_reduce1_chunks_r8(r8buf,rwrk,n_data,'MPI_SUM',20)
      endif
      call infolog('nubeam_cxbuf_reduce: nubeam_reduce_chunks_r8 done')

!  write out orball buffer / append ntr% data

      if (debug_mpi.and.(myidd.eq.0).and.(numprocd.gt.1)) then
         
         write(fname,'(a,a15,i1,a4)') trim(workpath)//trim(runid),'-orb-cxbuf-out-',lsbeam,'.tst'
         call nbi_buf_netcdf('W',trim(adjustl(fname)),n_data,r8buf1)
      endif


      if (myidd.eq.0) then

       call nbi_pget_ntr

      endif

      nmax_data=0
      n_data=0
      if(allocated(r8buf)) deallocate(r8buf)
      if(allocated(r8buf1)) deallocate(r8buf1)
      call infolog('nubeam_cxbuf_reduce: nubeam_buf_reduce done')
#endif      
      

      end subroutine nubeam_cxbuf_reduce 
!---------------------------------------------------------------------
!*********************************************************************
      subroutine nbi_ptcl_wght_reduce 
!
      use nbi_com
      use nbi_mpi_com
!
      implicit none
      real(fp),dimension(:), allocatable :: wrk_ptcl
      integer mpierr
!
#ifdef __MPI
      if(allocated(wrk_ptcl)) deallocate(wrk_ptcl)
      allocate(wrk_ptcl(mj))
      wrk_ptcl=0._fp
      mpierr=0
      call mpi_reduce(nmc(1:mj,lsbeam),wrk_ptcl,mj,MPI_REAL8,  &
           MPI_SUM,0,nbi_comm,mpierr)
      if(myidd.eq.0) then
            nmc(:,lsbeam)=wrk_ptcl(:)
      endif
      wrk_ptcl=0._fp
      mpierr=0
      call mpi_reduce(wnmc(1:mj,lsbeam),wrk_ptcl,mj,MPI_REAL8,  &
           MPI_SUM,0,nbi_comm,mpierr)
      if(myidd.eq.0) then
            wnmc(:,lsbeam)=wrk_ptcl(:)
      endif
      wrk_ptcl=0._fp
      mpierr=0
      call mpi_reduce(wghtot,wrk_ptcl(1:1),1,MPI_REAL8,  &
           MPI_SUM,0,nbi_comm,mpierr)
      if(myidd.eq.0) then
           wghtot =wrk_ptcl(1)
      endif
      deallocate(wrk_ptcl)
#endif
      
    end subroutine nbi_ptcl_wght_reduce
    
    subroutine nbi_ebav_reduce 
!
      use nbi_com
      use nbi_mpi_com
!
      implicit none
      real(fp),dimension(:), allocatable :: wrk_ptcl
      integer mpierr
!
#ifdef __MPI
      if(allocated(wrk_ptcl)) deallocate(wrk_ptcl)
      allocate(wrk_ptcl(mj))
      wrk_ptcl=0._fp
      mpierr=0
      call mpi_reduce(ebapl2sw(1:mj,lsbeam),wrk_ptcl,mj,MPI_REAL8,  &
           MPI_SUM,0,nbi_comm,mpierr)
      if(myidd.eq.0) then
            ebapl2sw(:,lsbeam)=wrk_ptcl(:)
      endif
      wrk_ptcl=0._fp
      mpierr=0
      call mpi_reduce(ebapp2sw(1:mj,lsbeam),wrk_ptcl,mj,MPI_REAL8,  &
           MPI_SUM,0,nbi_comm,mpierr)
      if(myidd.eq.0) then
            ebapp2sw(:,lsbeam)=wrk_ptcl(:)
      endif
      wrk_ptcl=0._fp
      mpierr=0
      call mpi_reduce(bdenss(1:mj,lsbeam),wrk_ptcl,mj,MPI_REAL8,  &
           MPI_SUM,0,nbi_comm,mpierr)
      if(myidd.eq.0) then
            bdenss(:,lsbeam)=wrk_ptcl(:)
      endif
      deallocate(wrk_ptcl)
#endif
      
    end subroutine nbi_ebav_reduce
     
    subroutine nbi_aep_reduce
      use nbi_com
      use nbi_aep_com,only : np_e,np_p,np_m, pbaep_ms,nkaep_use
      use nbi_aep_com,only :  pbaep_pol,pph_ms,aep_count,n_norm_pol
      use nbi_mpi_com
      use logmod
!
      implicit none
      real(fp) :: wrk_ptcl,w_naeptot
      integer mpierr,i,inum
      real(fp), dimension(:,:,:), allocatable :: zpbaep
!
#ifdef __MPI
      if(allocated(pbaep_pol)) then ! it is not allocated for RBQ model
         wrk_ptcl=0._fp
         mpierr=0
         w_naeptot=real(naeptot(lsbeam))
         call mpi_reduce(w_naeptot,wrk_ptcl,1,MPI_REAL8,  &
              MPI_SUM,0,nbi_comm,mpierr)
         if(myidd.eq.0) then
            naeptot(lsbeam) =wrk_ptcl
         endif
      endif
      if(nltest_output) then
         inum = np_e*np_p*np_m
         if(inum.gt.0) then
            if(allocated(zpbaep)) deallocate(zpbaep)
            allocate(zpbaep(np_e,np_p,np_m))
            mpierr=0
            call infolog('nbi_aep_reduce: nbi_aep_reduce start pbaep_ms')
            do i=1,nkaep_use 
               zpbaep=0.0_fp
               call mpi_reduce(pbaep_ms(1:np_e,1:np_p,1:np_m,i:i,1:1),zpbaep,inum,MPI_REAL8,&
                    MPI_SUM,0,nbi_comm,mpierr)
               if(myidd.eq.0) then
                  pbaep_ms(1:np_e,1:np_p,1:np_m,i,1)=zpbaep(1:np_e,1:np_p,1:np_m)
               endif
               zpbaep=0.0_fp
               call mpi_reduce(aep_count(1:np_e,1:np_p,1:np_m,i:i,1:1),zpbaep,inum,MPI_REAL8,&
                    MPI_SUM,0,nbi_comm,mpierr)
               if(myidd.eq.0) then
                  aep_count(1:np_e,1:np_p,1:np_m,i,1)=zpbaep(1:np_e,1:np_p,1:np_m)
               endif
               zpbaep=0.0_fp
               call mpi_reduce(pph_ms(1:np_e,1:np_p,1:np_m,i:i,1:1),zpbaep,inum,MPI_REAL8,&
                    MPI_SUM,0,nbi_comm,mpierr)
               if(myidd.eq.0) then
                  pph_ms(1:np_e,1:np_p,1:np_m,i,1)=zpbaep(1:np_e,1:np_p,1:np_m)
               endif
            enddo
         endif
         if(allocated(pbaep_pol)) then ! it is not allocated for RBQ model
            inum=nkaep_use*n_norm_pol
            if(inum.gt.0) then
               if(allocated(zpbaep)) deallocate(zpbaep)
               allocate(zpbaep(n_norm_pol,nkaep_use,1))
               zpbaep=0.0_fp
               call mpi_reduce(pbaep_pol(1:n_norm_pol,1:nkaep_use,1:1),&
                    zpbaep(1:n_norm_pol,1:nkaep_use,1:1),&
                    inum,MPI_REAL8,MPI_SUM,0,nbi_comm,mpierr)
               if(myidd.eq.0) then
                  pbaep_pol(1:n_norm_pol,1:nkaep_use,1:1)=zpbaep(1:n_norm_pol,1:nkaep_use,1:1)
               endif
            endif
         endif
      endif
      if(allocated(zpbaep)) deallocate(zpbaep)
#endif
    end subroutine nbi_aep_reduce

!*********************************************************************
      subroutine nbi_ptcl_wght_bcast
!
      use nbi_com
      use nbi_mpi_com
      use logmod
!
      implicit none
      integer :: mpierr,iost
#ifdef __MPI
      mpierr=0
      call mpi_bcast(wbav_kill(1:mj,lsbeam),mj,MPI_REAL8,0,nbi_comm,mpierr)
      if (mpierr.ne.0) then
         call errorlog(' ??nbi_ptcl_wght_bcast: mpi_bcast error - wbav_kill')
         
         call errset_mpi(myidd,5)
      endif
      mpierr=0
      call mpi_bcast(wbav(1:mj,lsbeam),mj,MPI_REAL8,0,nbi_comm,mpierr)
      if (mpierr.ne.0) then
         call errorlog(' ??nbi_ptcl_wght_bcast: mpi_bcast error - wbav')
         
         call errset_mpi(myidd,5)
      endif
#endif
    end subroutine nbi_ptcl_wght_bcast
!---------------------------------------------------------------------
    subroutine beam_in_box_reduce(il_3dhalo_den,ierr)
      use nbi_com
      use nbi_mpi_com
      use logmod
      implicit none
      logical, intent(in) :: il_3dhalo_den ! if 3d halo model is set
                           ! and all beams have been set in 3d-boxes
                           ! then il_3dhalo_den will be set .true.
                           ! and halo density, vol av halo e and 
                           ! vol av halo omega from 3d halo model 
                           ! (not frantic one, set in parms.for)
                           ! will be used
      integer,intent(out) :: ierr
      real(fp),dimension(:), allocatable :: wrk      
      integer,dimension(:), allocatable :: iwrk
      integer :: mpierr,isize,i,j,k
      real(fp) :: rwrk(1)
      real(fp), dimension(:,:,:), allocatable :: zfbm
#ifdef __MPI
      ierr=0
      isize=nbbox*nxbox*nybox*nlbox
      if(levmod_halo.eq.2.and.nsplt.gt.0) then
         isize=isize*3
      endif
      if(isize.ne.size(boxn0).or.isize.le.0) then
         ierr=1
         call errorlog(' ??beam_in_box_reduce: size boxn0 does not match on myidd=', myidd)
         call errset_mpi(myidd,ierr)
         return
      endif
      if(allocated(wrk)) deallocate(wrk)
      allocate(wrk(isize))
      wrk=0._fp
      mpierr=0
      call infolog('beam_in_box_reduce: nubeam_allreduce_chunks_r8 start')
      call nubeam_allreduce_chunks_r8(boxn0,wrk,isize,'MPI_SUM',35)
      call infolog('beam_in_box_reduce: nubeam_allreduce_chunks_r8 done')
!mg      call mpi_allreduce(boxn0,wrk,isize,MPI_REAL8,MPI_SUM,nbi_comm,mpierr)
!mg      if (mpierr.ne.0) then
!mg         ierr=2
!mg         call errorlog(' ??beam_in_box_reduce: mpi_reduce error - boxn0')
!mg         call errset_mpi(myidd,ierr)
!mg         if(allocated(wrk)) deallocate(wrk)
!mg         return
!mg      endif

      boxn0=wrk !now total statistic 'ndepbox_use' is there on all cpu

      if(allocated(wrk)) deallocate(wrk) 

      if(myidd.eq.0) then
         if(allocated(zfbm)) deallocate(zfbm)
         allocate(zfbm(mj,mig,mibs))
         zfbm=0.0_fp
      endif
      isize = mj*mig*mibs

      mpierr=0
      call infolog('beam_in_box_reduce: sn0s3dh')
      if(myidd.eq.0) then
         call nubeam_reduce_chunks_r8(sn0s3dh(1,1,1),zfbm,isize,'MPI_SUM',41)
      else
         rwrk=0.0_fp
         call nubeam_reduce1_chunks_r8(sn0s3dh(1,1,1),rwrk,isize,'MPI_SUM',41)
      endif
      call infolog('beam_in_box_reduce: done sn0s3dh')

      if(myidd.eq.0) then
         do k=1,mibs
            do j=1,mig
               do i=1,mj
                  sn0s3dh(i,j,k)=zfbm(i,j,k)
               enddo
            enddo
         enddo
      endif
      if(il_3dhalo_den) then
         if(allocated(zfbm)) deallocate(zfbm)
         allocate(zfbm(mig,mj,mis))
         zfbm=0.0_fp
         isize = mj*mig*mis
         
         mpierr=0
         call infolog('beam_in_box_reduce: nbi_den0mn')
         call nubeam_allreduce_chunks_r8(nbi_den0mn(1,1,1),zfbm,isize,'MPI_SUM',41)
         call infolog('beam_in_box_reduce: done nbi_den0mn')
         
         nbi_den0mn=zfbm !#
         
         isize = mj*mig
         
         mpierr=0
         call infolog('beam_in_box_reduce: vol')
         zfbm=0.0_fp
         call nubeam_allreduce_chunks_r8(voln0(1,1),zfbm(:,:,1),isize,'MPI_SUM',42)
         voln0(:,:)=zfbm(:,:,1)
         zfbm=0.0_fp
         call nubeam_allreduce_chunks_r8(evol0(1,1),zfbm(:,:,1),isize,'MPI_SUM',42)
         evol0(:,:)=zfbm(:,:,1)
         zfbm=0.0_fp
         call nubeam_allreduce_chunks_r8(ovol0(1,1),zfbm(:,:,1),isize,'MPI_SUM',42)
         ovol0(:,:)=zfbm(:,:,1)
         call infolog('beam_in_box_reduce: done vol')
         
         if(allocated(zfbm)) deallocate(zfbm)
         allocate(zfbm(mig,mimxbz,1))
         
         isize = mimxbz*mig
         mpierr=0
         call infolog('beam_in_box_reduce: vol2')
         zfbm=0.0_fp
         call nubeam_allreduce_chunks_r8(voln02(1,1),zfbm(:,:,1),isize,'MPI_SUM',43)
         voln02(:,:)=zfbm(:,:,1)
         zfbm=0.0_fp
         call nubeam_allreduce_chunks_r8(evol02(1,1),zfbm(:,:,1),isize,'MPI_SUM',43)
         evol02(:,:)=zfbm(:,:,1)
         zfbm=0.0_fp
         call nubeam_allreduce_chunks_r8(ovol02(1,1),zfbm(:,:,1),isize,'MPI_SUM',43)
         ovol02(:,:)=zfbm(:,:,1)
         call infolog('beam_in_box_reduce: done vol2')
      endif

      if(allocated(zfbm)) deallocate(zfbm)

      if(allocated(iwrk)) deallocate(iwrk)      
      allocate(iwrk(nbbox))
      iwrk=0._fp
      mpierr=0
      call mpi_allreduce(nsample_box,iwrk,nbbox,MPI_INTEGER,MPI_SUM,nbi_comm,mpierr)
      if (mpierr.ne.0) then
         ierr=2
         call errorlog(' ??beam_in_box_reduce: mpi_reduce error - nsample_box')
         call errset_mpi(myidd,ierr)
         if(allocated(iwrk)) deallocate(iwrk)
         return
      endif

      nsample_box=iwrk
      if(allocated(iwrk)) deallocate(iwrk)      
      
#else
#endif
      
    end subroutine beam_in_box_reduce
!---------------------------------------------------------------------
     subroutine orball_lost_gather
!
      use nbi_com
      use nbi_mpi_com
      use nbi_pllbuf
      use nbi_lost_ptcl, only: mlost_all    !arrays containing the number of elements 
                                            !that are received from each process

      use nbi_lost_ptcl, only: num_lst_arrays !number of arrays to gather lost ptcl
      use logmod

      implicit none
      integer :: mpierr,i,ii
      integer, dimension(:), allocatable :: idispls  ! displacement relative to r8buf1
                                                     ! at which to place the incoming data from 
                                                     ! process  i 
      integer, dimension(:), allocatable :: nn_data  !total dimension of receiving buffer
#ifdef __MPI
! gather #of lost ptcls. from all cpu to all cpu
      mpierr=0
      if(mlost_data.eq.0) then
         n_data=0
      else
! put data of lost ptcls. to buffer
         call nbi_pput_lst
      endif
      if (.not.allocated(nn_data))  allocate(nn_data(numprocd))
      nn_data = 0
      
      call mpi_allgather(n_data,1,MPI_INTEGER,nn_data,1,MPI_INTEGER,nbi_comm,mpierr)

      if (sum(nn_data).eq.0) then
         if (allocated(nn_data))  deallocate(nn_data)
         call infolog('orball_lost_gather: no data to be gathered')
         return                                        !no data to be gathered
      endif

! allocate buffer to receive data from all cpu
      call infolog('orball_lost_gather start')
      call nubeam_gatherv_chunks(nn_data)
      call infolog('orball_lost_gather: lost ptcl. gather done')
      
      if(myidd.eq.0) then
         if(allocated(mlost_all))deallocate(mlost_all)
         allocate(mlost_all(numprocd))
         mlost_all = nn_data/num_lst_arrays  !num_lst_arrays not supposed to be 0, check nbspec.dat
                                             !for a presence of lst% arrays
         call nbi_pget_lst !cpu0 - allocate and get lost data to mlost_data from all cpu
                           !now mlost_data=sum(mlost_all) on cpu0 only
      endif
      nmax_data=0
      n_data=0
      if(myidd.eq.0) then
         if(allocated(r8buf1))deallocate(r8buf1)
         if(allocated(mlost_all))deallocate(mlost_all)
      endif
      if(allocated(r8buf))deallocate(r8buf)
      if(allocated(nn_data))deallocate(nn_data)
#endif
      end subroutine orball_lost_gather
!---------------------------------------------------------------------
     subroutine depall_buf_reduce 
!
      use nbi_com
      use nbi_mpi_com
      use nbi_pllbuf
      use logmod

      implicit none

      integer :: mpierr,jk
      character(len=ncfnm)   :: fname 

      call nbi_pput_dep

#ifdef __MPI
      mpierr=0
! that combines the values stored by each process into a single value 
! available to all processes
      call infolog('depall_buf_reduce: nubeam_allreduce_chunks_r8 start')
      call mpi_allreduce(r8buf,r8buf1,n_data,MPI_REAL8,MPI_SUM,nbi_comm,mpierr)
!      call nubeam_allreduce_chunks_r8(r8buf,r8buf1,n_data,'MPI_SUM',25)
      call infolog('depall_buf_reduce: nubeam_allreduce_chunks_r8 done')

      if (debug_mpi.and.(myidd.eq.0).and.nltest_output) then
! for mpi nubeam debugging write to the file 
! summarized over all prosessors data from the deposition.

         write(fname,'(a,a13,i1,a4)') trim(workpath)//trim(runid),'-dep-buf-out-',lsbeam,'.tst'
         call nbi_buf_netcdf('W',trim(adjustl(fname)),n_data,r8buf1)
      endif

#else
#endif


      call nbi_pget_dep

      nmax_data=0
      n_data=0
      deallocate(r8buf,r8buf1)

      end subroutine depall_buf_reduce 
!---------------------------------------------------------------------
!---------------------------------------------------------------------
     subroutine depall_dpt_buf_gather
!
      use nbi_com
      use nbi_mpi_com
      use nbi_pllbuf
      use logmod
      implicit none
      integer, dimension(:), allocatable :: nn_data  !total dimension of receiving buffer
      integer :: mpierr,i
      integer :: mtrk_0(mibs_use),mwrk(mibs_use)
      integer :: ijmem,fsystem_mem

      mwrk=0
#ifdef __MPI
      mtrk_0=0 !total# of zones
      mtrk_0(1:mibs_use)=sum(xntrack(1:mdep0,1:mibs_use),dim=1)+0.5_fp!cpu's # of track for each species 
      call  mpi_reduce(mtrk_0,mwrk,mibs_use,MPI_INTEGER,MPI_SUM,0,nbi_comm,mpierr)
      call nbi_pput_dpt
      !now all track data in black box and deallocated
      mtrk_0=mwrk
      mtrk0=sum(mtrk_0(1:mibs_use)) ! only cpu0 knows total number,  child's setting mtrk0=0

! that combines the values stored by each process into a single value 
! available to all processes

      if (.not.allocated(nn_data))  allocate(nn_data(numprocd))
      nn_data(:) = 0
      call infolog('depall_dpt_buf_gather: n_data = ',n_data)
      call mpi_allgather(n_data,1,MPI_INTEGER,nn_data,1,MPI_INTEGER,nbi_comm,mpierr)

      if(myidd.eq.0) then
         do i=1,numprocd
           call infolog('depall_dpt_buf_gather: nn_data = ',nn_data(i)) 
         enddo
      endif
      if (sum(nn_data).eq.0) then
         if (allocated(nn_data))  deallocate(nn_data)
         call infolog('depall_dpt_buf_gather: no data to be gathered')
         return                                        !no data to be gathered
      endif

! allocate buffer to receive data from all cpu

      call warnlog('depall_dpt_buf_gather: start')
      call nubeam_gatherv_chunks(nn_data)
      call warnlog('depall_dpt_buf_gather: track data gather done')

      call nbi_dalloc_track ! to deallocate track arrays for  all processes
      call nbi_dalloc_track_zone ! to deallocate track arrays for  all processes
      call infolog('depall_dpt_buf_gather: dealloc_track')
         
      if(myidd.eq.0) then
         call warnlog('depall_dpt_buf_gather: start nbi_pget_dpt')
         call nbi_pget_dpt(mtrk_0) ! allocate track arrays for full list of track (ndep0), only on cpu0
                           ! and read data from buffer
         call warnlog('depall_dpt_buf_gather: data written to netcdf file')
         
      endif


      nmax_data=0
      n_data=0
      if(allocated(r8buf1))deallocate(r8buf1)
      if(allocated(r8buf))deallocate(r8buf)
      if(allocated(nn_data))deallocate(nn_data)

      ijmem = fsystem_mem() ! print out memmory usage

#else
#endif

      end subroutine depall_dpt_buf_gather
!---------------------------------------------------------------------
!---------------------------------------------------------------------
      subroutine deptube_buf_reduce 
!
      use nbi_com
      use nbi_mpi_com
      use tube_mod

      implicit none

      integer :: ier,mpierr,inumdata,inum1

      call deptube_buf_put(ier) 

      inumdata=size(deptubebuf)
      inum1=size(deptubebuf1)

#ifdef __MPI
      mpierr=0
      call infolog('deptube_buf_reduce: nubeam_allreduce_chunks_r8 start')
      call nubeam_allreduce_chunks_r8(deptubebuf,deptubebuf1,inumdata,'MPI_SUM',26)
      call infolog('deptube_buf_reduce: nubeam_allreduce_chunks_r8 done')

#else
#endif

      ier=0
      call deptube_buf_get(ier)

      end subroutine deptube_buf_reduce 
!---------------------------------------------------------------------
      subroutine depall_mpi_split(innp1,idep0,nusey) 
!
!     particles from innp1 to idep0 are spread among processors
!     myidd, numprocd from nbi_mpi_com
!     nusey - #ptcl on each processor
! mg 03/2009

        use nbi_com
        use nbi_mpi_com
        use dep_mpi_split
        
        implicit none
        integer, intent(in) :: innp1,idep0
        integer :: idt,idk
        integer, intent(out) :: nusey
       
        nusey=0
        if(idep0.lt.innp1) return
        
        do idt=innp1+myidd,idep0,numprocd
           nusey=nusey+1
        enddo
        if (allocated(ninddumn)) deallocate(ninddumn)
        allocate(ninddumn(nusey))
        
        ninddumn=0.0
        
        idk=0
        do idt=innp1+myidd,idep0,numprocd
           idk=idk+1
           ninddumn(idk)=idt 
        enddo
        
      end subroutine depall_mpi_split
!---------------------------------------------------------------------
      subroutine depall_buf_read 
!
      use nbi_com
      use nbi_mpi_com
      use nbi_pllbuf
      use logmod
      

      implicit none

      integer :: lk,nn_data,read_stat,iost
      character(len=ncfnm)   :: fname 
      integer, parameter ::eor=-2

      nn_data=0
      write(fname,'(a,a13,i1,a4)')  trim(workpath)//trim(runid),'-dep-buf-out-',lsbeam,'.tst'
      call nbi_buf_netcdf('d',trim(adjustl(fname)),nn_data,r8buf1) !nn_data will be read, r8buf1 nor allocated yet

      if(nn_data.le.0) then
         call errorlog(' ?depall_buf_read: error reading deposition buf')
         return
      endif
      if (allocated(r8buf1)) deallocate(r8buf1)
      allocate(r8buf1(nn_data))

      call nbi_buf_netcdf('R',trim(adjustl(fname)),nn_data,r8buf1) !r8buf1 will be read

      call nbi_pget_dep

      if (allocated(r8buf1)) deallocate(r8buf1)

      end subroutine depall_buf_read 
!---------------------------------------------------------------------
!----------------------------------------------------------------------
      subroutine nbi_ptcl_write(suffix_in)
! 
!  suffix_in - filename suffix,
!  write file - workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>.cdf or,
!               if last 3 letters of <suffix_in> are tmp
!               then write to the file
!               /workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>_cdf.tmp for transp.
!               workpath, runid - from nbi_com
!               myid, numprocd  - from nbi_mpi_com
!  mg 03/2009
!
        use ezcdf
        use nbi_com
        use nbi_mpi_com
        use logmod
        
        implicit none
        
        character*(*), intent(in) :: suffix_in  ! ptcl state filenames suffix
        character*150 :: suffix  ! ptcl state filenames
        integer :: jcdf,ierr,lsuff,i_wrk,mpierr,iertmp
        character*150 :: fname, fname_list
        character*300 :: zcmd
        character*9 cmyid,cnumprocs
        integer :: ilun_tmp,iost,iwarn
        logical :: ex
!--------------------------
        if(nbsteps.eq.0) then
           call infolog('nbi_ptcl_write: nbsteps=0 there is no particle state to write')
           
           return                ! no steps have been completed, so,
                                 ! there is no particle state to read, yet...
        endif
        if(.not.nltest_output) return
        cmyid=' '
        write(cmyid,'(i5)') myidd
        cnumprocs=' '
        write(cnumprocs,'(i5)') numprocd
        
        lsuff=len_trim(workpath)
        if(lsuff.eq.0) then
           workpath='./'
        else
           suffix=workpath(lsuff:lsuff)
           if(trim(suffix).ne.'/') then
              workpath=trim(adjustl(workpath))//'/'
           endif
        endif

!        call find_io_unit(ilun_tmp)


!        fname=trim(adjustl(workpath))//'file-list-'//trim(adjustl(cmyid))// &
!           '-'//trim(adjustl(cnumprocs))//'.dat'
!        open(unit=ilun_tmp,file=trim(fname),status='unknown')
        
!        close(unit=ilun_tmp)

 
       
        suffix=suffix_in
        fname=' '
        call default_filename_ptcl(workpath,runid,fname,suffix,myidd,numprocd)
        call infolog('nbi_ptcl_write: ptcl. write to '//trim(fname))
!mg,06/2019  before writing ptcl state file delete old one if exists
        inquire(file=trim(fname),exist=ex)
        if(ex) then
           call fdelete(trim(fname),iwarn)
        endif

!mg,06/2019  no more         call frename_ok(trim(fname),trim(fname)//'~',iwarn) !rename to *cdf~

        if(len(trim(fname)).gt.0) then
           ierr=0
           call ezcdf_open(jcdf,fname,'W',ierr)
           if(ierr.eq.0) then
              call nbi_state_write_ptcl(jcdf)
              call ezcdf_close(jcdf)
              call infolog('nbi_ptcl_write: write ok to file = '//trim(fname))
              
              
           else
              call errorlog('nbi_ptcl_write: error on write to file = '//trim(fname))
           endif
        else
           ierr=1
           call errorlog(' nbi_ptcl_write: nbi_states(....,"s",...) has blank ptcl state filename')
           call errorlog(' nbi_ptcl_write: filename must be provided')
        endif
           
! take an action if error flag /= 0           
#ifdef __MPI
        ! reduce  error flags
        mpierr=0
        i_wrk=0
        call  mpi_reduce(ierr,i_wrk,1,MPI_INTEGER,  &
             MPI_SUM,0,nbi_comm,mpierr)
        if(mpierr.ne.0) then
           call errorlog(' nbi_ptcl_write: mpi_reduce failed. myid= ',myid)
           call errset_mpi(myid,111)
        endif
        if(myidd.eq.0) then
           if(i_wrk.ne.0)  then
              ierr = max(1,i_wrk)
              write(0,*) ' nbi_ptcl_write:  error -- not all writes were successful'        
              write(0,*) ' nbi_ptcl_write:  myid= 0 calls mpi_abort'        
              call mpi_abort(nbi_comm,ierr,iertmp)
              call exit(1)
           endif
        endif
#else
        if(ierr.ne.0) then
           call errset_mpi(myid,111)
        endif
#endif

      end subroutine nbi_ptcl_write
!---------------------------------------------------------------------
!---------------------------------------------------------------------
      subroutine nbi_ptcl_delete(suffix_in)
!
!  suffix_in - filename suffix
!  delete  file - workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>.cdf or,
!               if last 3 letters of <suffix_in> are tmp
!               then read file
!               /workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>_cdf.tmp for transp,
!               or try /workpath/<runid>//<suffix_in> for nubeam_test debugging.
!               workpath, runid - from nbi_com
!               myid, numprocd  - from nbi_mpi_com
!  mg 03/2009
!
        use ezcdf
        use nbi_com
        use nbi_mpi_com
        use logmod
        
        implicit none
        
        character*(*), intent(in) :: suffix_in  ! state filenames
        integer :: jcdf,ierr,lsuff
        character*150 fname,suffix
        character*9 cmyid,cnumprocs
        character*200 zcmd
        integer :: iwarn

        if(nbsteps.eq.0) then
           call infolog('nbi_ptcl_read: nbsteps=0 there is no particle state to delete')
           
           return                ! no steps have been completed, so,
                                 ! there is no particle state to read, yet...
        endif
        cmyid=' '
        write(cmyid,'(i5)') myidd
        cnumprocs=' '
        write(cnumprocs,'(i5)') numprocd
        lsuff=len_trim(workpath)
        if(lsuff.eq.0) then
           workpath='./'
        else
           suffix=workpath(lsuff:lsuff)
           if(trim(suffix).ne.'/') then
              workpath=trim(adjustl(workpath))//'/'
           endif
        endif
        suffix=suffix_in
        fname=' '
        call default_filename_ptcl(workpath,runid,fname,suffix,myidd,numprocd)

        if(len(trim(fname)).gt.0) then
           call infolog('nbi_ptcl_delete: '//trim(fname))
           zcmd = 'rm -f '//trim(fname)//' '//trim(fname)//'_save'
	   call fdelete(trim(fname),iwarn)
           if (iwarn/=0) then
              call errorlog('nbi_ptcl_delete: failed to delete' &
		//trim(fname))
           endif	
	   call fdelete(trim(fname)//'_save',iwarn)
           if (iwarn/=0) then
              call errorlog('nbi_ptcl_delete: failed to delete' &
		//trim(fname)//'_save')
           endif	
        endif
          
      end subroutine nbi_ptcl_delete
!---------------------------------------------------------------------
!---------------------------------------------------------------------
      subroutine nbi_ptcl_backup(suffix_in)
!
!  suffix_in - filename suffix
!  backup  file - workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>.cdf or,
!               if last 3 letters of <suffix_in> are tmp
!               then read file
!               /workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>_cdf.tmp for transp,
!               or try /workpath/<runid>//<suffix_in> for nubeam_test debugging.
!               workpath, runid - from nbi_com
!               myid, numprocd  - from nbi_mpi_com
!               to workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>.cdf_save or
!               /workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>_cdf.tmp_save for transp
!  mg 03/2009
!
        use ezcdf
        use nbi_com
        use nbi_mpi_com
        use logmod
      
        implicit none
        
        character*(*), intent(in) :: suffix_in  ! state filenames
        integer :: jcdf,ierr,lsuff
        character*150 fname,suffix
        character*9 cmyid,cnumprocs
        character*200 zcmd
        integer :: iwarn

        cmyid=' '
        write(cmyid,'(i5)') myidd
        cnumprocs=' '
        write(cnumprocs,'(i5)') numprocd
        lsuff=len_trim(workpath)
        if(lsuff.eq.0) then
           workpath='./'
        else
           suffix=workpath(lsuff:lsuff)
           if(trim(suffix).ne.'/') then
              workpath=trim(adjustl(workpath))//'/'
           endif
        endif
        suffix=suffix_in
        fname=' '
        call default_filename_ptcl(workpath,runid,fname,suffix,myidd,numprocd)

        if(len(trim(fname)).gt.0) then
           call infolog('nbi_ptcl_backup: '//trim(fname)//'_save')
           
           zcmd = 'cp '//trim(fname)//' '//trim(fname)//'_save'
	   call fcopy(trim(fname),trim(fname)//'_save',iwarn)
           if(iwarn.ne.0) then
              ierr=1
              call errorlog('nbi_ptcl_backup(...,"b",...): backup failure')
              call errorlog('    fname was: '//trim(fname))
           endif
        endif
          
      end subroutine nbi_ptcl_backup
!---------------------------------------------------------------------
!---------------------------------------------------------------------
      subroutine nbi_ptcl_retrieve(suffix_in)
!
!  suffix_in - filename suffix
!  retrieve  file - workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>.cdf_save or,
!               if last 3 letters of <suffix_in> are tmp
!               then read file
!               /workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>_cdf.tmp_save for transp,
!               or try /workpath/<runid>//<suffix_in> for nubeam_test debugging.
!               workpath, runid - from nbi_com
!               myid, numprocd  - from nbi_mpi_com
!  mg 03/2009
!
        use ezcdf
        use nbi_com
        use nbi_mpi_com
        use logmod
        
        implicit none
        
        character*(*), intent(in) :: suffix_in  ! state filenames
        integer :: jcdf,ierr,lsuff
        character*150 fname,suffix
        character*9 cmyid,cnumprocs
        character*200 zcmd
        integer :: iwarn

        cmyid=' '
        write(cmyid,'(i5)') myidd
        cnumprocs=' '
        write(cnumprocs,'(i5)') numprocd
        lsuff=len_trim(workpath)
        if(lsuff.eq.0) then
           workpath='./'
        else
           suffix=workpath(lsuff:lsuff)
           if(trim(suffix).ne.'/') then
              workpath=trim(adjustl(workpath))//'/'
           endif
        endif
        suffix=suffix_in
        fname=' '
        call default_filename_ptcl(workpath,runid,fname,suffix,myidd,numprocd)

        if(len(trim(fname)).gt.0) then
           call infolog('nbi_ptcl_retrieve: '//trim(fname)//'_save')
           
           zcmd = 'cp '//trim(fname)//'_save  '//trim(fname)
           call fcopy( trim(fname)//'_save', trim(fname), iwarn)
           if(iwarn.ne.0) then
              ierr=1
              call errorlog('nbi_ptcl_retrieve(...,"c",...): retrieve failure')
              call errorlog('    fname was: '//trim(fname))
           endif
        endif
          
      end subroutine nbi_ptcl_retrieve
!---------------------------------------------------------------------
!---------------------------------------------------------------------
      subroutine nbi_ptcl_read(suffix_in)
!
!  suffix_in - filename suffix
!  read  file - workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>.cdf or,
!               if last 3 letters of <suffix_in> are tmp
!               then read file
!               /workpath/<runid>_<suffix_in>_cpu<myid>_<numprocd>_cdf.tmp for transp,
!               or try /workpath/<runid>//<suffix_in> for nubeam_test debugging.
!               workpath, runid - from nbi_com
!               myid, numprocd  - from nbi_mpi_com
!  mg 03/2009
!
        use ezcdf
        use nbi_com
        use nbi_mpi_com
        use logmod
  
        implicit none
        
        character*(*), intent(in) :: suffix_in  ! state filenames
        integer :: jcdf,ierr,ilen,lsuff,iost
        character*150 fname,suffix
        character*9 cmyid,cnumprocs

        if(nbsteps.eq.0) then
           nuse=0
           ntotptp=0
           call infolog('nbi_ptcl_read: nbsteps=0 there is no particle state to read')
           
           return                ! no steps have been completed, so,
                                 ! there is no particle state to read, yet...
        endif
        cmyid=' '
        write(cmyid,'(i5)') myidd
        cnumprocs=' '
        write(cnumprocs,'(i5)') numprocd
        lsuff=len_trim(workpath)
        if(lsuff.eq.0) then
           workpath='./'
        else
           suffix=workpath(lsuff:lsuff)
           if(trim(suffix).ne.'/') then
              workpath=trim(adjustl(workpath))//'/'
           endif
        endif
        suffix=suffix_in
        fname=' '
        call default_filename_ptcl(workpath,runid,fname,suffix,myidd,numprocd)

        if(len(trim(fname)).gt.0) then
           call infolog('nbi_ptcl_read: '//trim(fname))
           
           call ezcdf_open(jcdf,fname,'r',ierr)
           
           if(ierr.ne.0) then
              ilen=len(trim(fname))
              fname(ilen+1:ilen+6) = '_save'
              call infolog('nbi_ptcl_read: retry, filename is: '//trim(fname))
              
              call ezcdf_open(jcdf,fname,'r',ierr)
              if(ierr.eq.0) then
                 call infolog('nbi_ptcl_read: file read ok.')
              else
                 call errorlog('nbi_ptcl_read: error reading file: '// &
                      trim(fname))
                 
#ifndef __MPI
! for nubeam_test  to debug mpi version as a serial version . try file = suffix
                 if(debug_mpi.and.nltest_output) then
                    fname=trim(workpath)//trim(runid)//trim(suffix)
                    call infolog('nbi_ptcl_read retry, filename is: '//trim(fname))
                    
                    call ezcdf_open(jcdf,fname,'r',ierr)
                    if(ierr.eq.0) then
                       call infolog('nbi_ptcl_read: file read ok.')
                       
                    else 
                       call errorlog(' ??nbi_ptcl_read: error reading file: '// &
                            trim(fname))
                       call errset_mpi(myidd,110)
                    endif
                 endif
#else
                 call errset_mpi(myidd,110)
#endif
              endif
           endif
           
           if (ierr.eq.0) then 
              
              call nbi_state_read_ptcl(jcdf)
              call ezcdf_close(jcdf)
              call infolog( 'nbi_ptcl_read: file read ok.')
           else
              call errorlog( 'nbi_ptcl_read: file does not appear to exist '//trim(fname(1:ilen)))
           endif
        else
           ierr=1
              call errorlog(' ??nbi_ptcl_read: a filename must be provided')
           
           call errset_mpi(myidd,110)
        endif
          
      end subroutine nbi_ptcl_read
!---------------------------------------------------------------------
!---------------------------------------------------------------------
      subroutine nubeam_buf_read 
!
        use nbi_com
        use nbi_mpi_com
        use nbi_pllbuf
        use logmod
        
        implicit none
        
        integer :: lk,nn_data
        character(len=ncfnm)   :: fname 
#ifdef __MPI
        
        write(fname,'(a,a13,i1,a4)') trim(workpath)//trim(runid),'-orb-buf-out-',lsbeam,'.tst'
        call nbi_buf_netcdf('d',trim(adjustl(fname)),nn_data,r8buf1) !nn_data will be read, r8buf1 nor allocated yet
        
        if(nn_data.le.0) then
           call errorlog(' ?nubeam_buf_read: error reading deposition buf')
           return
        endif
        if (allocated(r8buf1)) deallocate(r8buf1)
        if (allocated(r8buf1)) deallocate(r8buf1)
        allocate(r8buf1(nn_data))
        
        call nbi_buf_netcdf('R',trim(adjustl(fname)),nn_data,r8buf1) !r8buf1 will be read
        
        call nbi_pget_orb
#else
#endif
        if (allocated(r8buf1)) deallocate(r8buf1)
        
        
      end subroutine nubeam_buf_read
!---------------------------------------------------------------------
!---------------------------------------------------------------------
      subroutine nubeam_cxbuf_read 
        !
        use nbi_com
        use nbi_mpi_com
        use nbi_pllbuf
        use logmod
        
        implicit none
        
        integer :: lk,nn_data
        character(len=ncfnm)   :: fname 
#ifdef __MPI
        
        write(fname,'(a,a15,i1,a4)') trim(workpath)//trim(runid),'-orb-cxbuf-out-',lsbeam,'.tst'
        call nbi_buf_netcdf('d',trim(adjustl(fname)),nn_data,r8buf1) !nn_data will be read, r8buf1 nor allocated yet
        
        if(nn_data.le.0) then
           call errorlog(' ?nubeam_cxbuf_read: error reading deposition buf')
           return
        endif
        if (allocated(r8buf1)) deallocate(r8buf1)
        allocate(r8buf1(nn_data))
        
        call nbi_buf_netcdf('R',trim(adjustl(fname)),nn_data,r8buf1) !r8buf1 will be read
        
        call nbi_pget_ntr
#else
#endif
        if (allocated(r8buf1)) deallocate(r8buf1)
        
        
      end subroutine nubeam_cxbuf_read
      !---------------------------------------------------------------------
!---------------------------------------------------------------------
subroutine nubeam_xplasma_bcast

  ! broadcast xplasma data, used by nbi_interp_profiles to load nbi_com
  ! at this point, also broadcast the value of nbi_init_flag, converted
  ! to integer

  ! note, this routine, like most in nbi_mpisubs, is only called if 
  ! #ifdef __MPI holds...

  use nbi_com
  use nbi_mpi_com
  use nbi_bcastbuf
  use xplasma_obj_instance
  use logmod

  implicit none

  character*200 :: zfilename1

    !  pointers for pack/unpack test...
  integer, dimension(:), pointer :: ipack
  real(fp), dimension(:), pointer :: r8pack
  integer, dimension(3) :: intgr_buf
  integer :: ilenint,ilenreal

  integer mpierr, ier,iost
  integer, dimension(4) :: lintbuf

  !---------------------------------------------
  call infolog('nubeam_xplasma_bcast: start')
  
  !---------------------------------------------
  nullify(ipack,r8pack)
  intgr_buf=0
  lintbuf=0
  if (myidd.eq.0) then
      call xplasma_pack(s, ipack, r8pack, ier)
      if (ier.ne.0) then
         call errorlog(' ??nubeam_xplasma_bcast: ier.ne.0: xplasma_pack difficulty')
         call xplasma_error(s, ier, 0)
         call errset_mpi(myid,9)
      endif
      intgr_buf(1)=ipack(5)
      intgr_buf(2)=ipack(6)
      if(nbi_init_flag) then
         intgr_buf(3)=1
      else
        intgr_buf(3)=0
     endif
   endif

#ifdef __MPI
   call mpi_bcast(intgr_buf,3,MPI_INTEGER,0,nbi_comm,mpierr)
  if (mpierr.ne.0) then
     call errorlog(' ??nubeam_xplasma_bcast: mpi_bcast error - intgr_buf')
     call errset_mpi(myidd,10)
  endif
#endif

  ilenint= intgr_buf(1)
  ilenreal= intgr_buf(2)

  if (myidd.gt.0)  then
     nbi_init_flag = (intgr_buf(3).eq.1)
     allocate(ipack(ilenint),r8pack(ilenreal))
  endif

#ifdef __MPI
  mpierr=0
  call infolog('nubeam_xplasma_bcast: mpi_bcast  - ipack size: ',ilenint)
  call nubeam_bcast_chunks_int(ipack(1:ilenint),ilenint,11)
  call infolog('nubeam_xplasma_bcast: mpi_bcast  - ipack, done')

  mpierr=0
  call infolog('nubeam_xplasma_bcast: mpi_bcast  - r8pack size: ',ilenreal)
  call nubeam_bcast_chunks_r8(r8pack(1:ilenreal),ilenreal,12)
  call infolog('nubeam_xplasma_bcast: mpi_bcast  - r8pack, done')

#endif

  if (myidd.gt.0) then
   call eqm_clear
   call xplasma_unpack(s, ipack, r8pack, ier)
   if (ier.ne.0) then
     call errorlog(' ??nubeam_xplasma_bcast: ier.ne.0 xplasma_unpack difficulty')
     call xplasma_error(s, ier, 0)
     call errset_mpi(myidd,13)
   endif
  endif

  deallocate(ipack,r8pack)

end subroutine nubeam_xplasma_bcast
!---------------------------------------------
!--------------------------------------------------
subroutine nubeam_nbinputs_bcast_init

  use nbi_com
  use nbi_mpi_com
  use nbi_bcastbuf
  use logmod

  implicit none

  integer mpierr,iost
  integer, dimension(4) :: lintbuf

  !---------------------------------------------
  call warnlog(' %nubeam_nbinputs_bcast_init start, myidd=',myidd)
  !---------------------------------------------

! first put the dimensions in a buffer (only root does this)

  if (myidd.eq.0) then
     call nbi_initdim_bufput
  endif

! broadcast the buffer size and allocate the buffer on child processes

#ifdef __MPI
  mpierr=0
  call mpi_bcast(ni_data,1,MPI_INTEGER,0,nbi_comm,mpierr)
  if (mpierr.ne.0) then
     call errorlog('nubeam_nbinputs_bcast_init: ni_data broadcast error')
!     call mpi_finalize(mpierr)
!     call exit(1)
     call bad_exit
  endif
#endif

  if (myidd.gt.0) then
     if (allocated(bufint))deallocate(bufint)
     allocate(bufint(ni_data)) 
  endif

! now broadcast the buffer containing demensions

#ifdef __MPI
  mpierr=0
  call infolog('nubeam_nbinputs_bcast_init: integer broadcast size: ',ni_data)
  call nubeam_bcast_chunks_int(bufint,ni_data,13)
  call infolog('nubeam_nbinputs_bcast_init: integer broadcast done')
#endif

! if child process, retrieve the contents of the buffer
! containing the dimensions

  if (myidd.gt.0) then
     call nbi_initdim_bufget
  endif

  ni_data=0
  nmax_idata=0
  if (allocated(bufint)) deallocate(bufint) 

! now put other data into respective buffers (only root does this)

  if (myidd.eq.0) then
     call nbi_bufput_data
     lintbuf(1)=ni_data
     lintbuf(2)=nr_data
     lintbuf(3)=nc_data
     lintbuf(4)=nl_data
  endif

! broadcast the buffer containing buffer dimensions to all
! participating processes

#ifdef __MPI
  mpierr=0
  call mpi_bcast(lintbuf,4,MPI_INTEGER,0,nbi_comm,mpierr)
  if (mpierr.ne.0) then
     call errorlog('nubeam_nbinputs_bcast_init: lintbuf broadcast error')
!     call mpi_finalize(mpierr)
!     call exit(1)
     call bad_exit
  endif
#endif

! now allocate the buffers on child processes to receive the data

  if (myidd.gt.0) then
     ni_data=lintbuf(1)
     nr_data=lintbuf(2)
     nc_data=lintbuf(3)
     nl_data=lintbuf(4)
     if (allocated(bufint)) deallocate(bufint)
     if (allocated(bufr8))  deallocate(bufr8)
     if (allocated(bufchr)) deallocate(bufchr)
     if (allocated(buflgc)) deallocate(buflgc)
     allocate(bufint(ni_data),bufr8(nr_data),bufchr(csize*nc_data),buflgc(nl_data))
  endif

! mpi broadcast the rest of the data to all the participating processes.

#ifdef __MPI
  mpierr=0
  call infolog('nubeam_nbinputs_bcast_init: integer broadcast 2 size: ',ni_data)
  call nubeam_bcast_chunks_int(bufint,ni_data,14)
  call infolog('nubeam_nbinputs_bcast_init: integer broadcast 2 done ')

  mpierr=0
  call infolog('nubeam_nbinputs_bcast_init: real broadcast size: ',nr_data)
  call nubeam_bcast_chunks_r8(bufr8,nr_data,15)
  call infolog('nubeam_nbinputs_bcast_init: real broadcast done')

  mpierr=0
  call mpi_bcast(bufchr,csize*nc_data,mpi_character,0,nbi_comm,mpierr)
  if (mpierr.ne.0) then
     call errorlog('nubeam_nbinputs_bcast_init: character broadcast error')
!     call mpi_finalize(mpierr)
!     call exit(1)
     call bad_exit
  endif

  mpierr=0
  call mpi_bcast(buflgc,nl_data,mpi_logical,0,nbi_comm,mpierr)
  if (mpierr.ne.0) then
     call errorlog('nubeam_nbinputs_bcast_init: logical broadcast error')
!     call mpi_finalize(mpierr)
!     call exit(1)
     call bad_exit
  endif
#endif

  if (myidd.gt.0) then 
     call nbi_bufget_data
  endif

  ni_data=0
  nmax_idata=0
  nr_data=0
  nmax_rdata=0
  nc_data=0
  nmax_cdata=0
  nl_data=0
  nmax_ldata=0

  if (allocated(bufint)) deallocate(bufint)
  if (allocated(bufr8))  deallocate(bufr8)
  if (allocated(bufchr)) deallocate(bufchr)
  if (allocated(buflgc)) deallocate(buflgc)

end subroutine nubeam_nbinputs_bcast_init
!------------------------------------------------
!---------------------------------------------------
subroutine nubeam_nbinputs_bcast_step

  use nbi_com
  use nbi_mpi_com
  use nbi_bcastbuf
  use logmod

  implicit none

  integer mpierr
  integer, dimension(4) :: lintbuf

  !---------------------------------------------
  call infolog('nubeam_nbinputs_bcast_step: start')
  !---------------------------------------------

  ! first put the dimensions in a buffer (only root does this)

  if (myidd.eq.0) then
     call nbi_dim_fld_bufput 
  endif

  ! broadcast the buffer size and allocate the buffer on child processes

#ifdef __MPI
  mpierr=0
  call mpi_bcast(ni_data,1,MPI_INTEGER,0,nbi_comm,mpierr)
  if (mpierr.ne.0) then
     call errorlog('nubeam_nbinputs_bcast_step: ni_data broadcast error')
!     call mpi_finalize(mpierr)
!     call exit(1)
     call bad_exit
  endif
#endif
       
  if (myidd.gt.0) then
     if (allocated(bufint)) deallocate(bufint)
     allocate(bufint(ni_data)) 
  endif

  ! now broadcast the buffer containing dimensions

#ifdef __MPI
  mpierr=0
  call infolog('nubeam_nbinputs_bcast_step: integer broadcast 1 size: ',ni_data)
  call nubeam_bcast_chunks_int(bufint,ni_data,16)
  call infolog('nubeam_nbinputs_bcast_step: integer broadcast 1 done')
#endif

  ! if child process, retrieve the contents of the buffer
  ! containing the dimensions

  if (myidd.gt.0) then
     call nbi_dim_fld_bufget 
  endif

  ni_data=0
  nmax_idata=0
  if (allocated(bufint)) deallocate(bufint) 

  ! now put other data into respective buffers (only root does this)

  if (myidd.eq.0) then
     call nbi_fld_bufput
     lintbuf(1)=ni_data
     lintbuf(2)=nr_data
     lintbuf(3)=nc_data
     lintbuf(4)=nl_data
  endif

  ! broadcast the buffer containing buffer dimensions to all
  ! participating processes

#ifdef __MPI
  mpierr=0
  call mpi_bcast(lintbuf,4,MPI_INTEGER,0,nbi_comm,mpierr)
  if (mpierr.ne.0) then
     call errorlog('nubeam_nbinputs_bcast_step: lintbuf broadcast  error')
!     call mpi_finalize(mpierr)
!     call exit(1)
     call bad_exit
  endif
#endif

  ! now allocate the buffers on child processes to receive the data

  if (myidd.gt.0) then
     ni_data=lintbuf(1)
     nr_data=lintbuf(2)
     nc_data=lintbuf(3)
     nl_data=lintbuf(4)
     if (allocated(bufint)) deallocate(bufint)
     if (allocated(bufr8))  deallocate(bufr8)
     if (allocated(bufchr)) deallocate(bufchr)
     if (allocated(buflgc)) deallocate(buflgc)
     allocate(bufint(ni_data),bufr8(nr_data),bufchr(csize*nc_data),buflgc(nl_data))
  endif

  ! mpi broadcast the rest of the data to all the participating processes.

#ifdef __MPI
  mpierr=0
  call infolog('nubeam_nbinputs_bcast_step: integer broadcast 2: ',ni_data)
  call nubeam_bcast_chunks_int(bufint,ni_data,17)
  call infolog('nubeam_nbinputs_bcast_step: integer broadcast 2 done')

  mpierr=0
  call infolog('nubeam_nbinputs_bcast_step: real broadcast: ',nr_data)
  call nubeam_bcast_chunks_r8(bufr8,nr_data,18)
  call infolog('nubeam_nbinputs_bcast_step: real broadcast done')

  mpierr=0
  call mpi_bcast(bufchr,csize*nc_data,mpi_character,0,nbi_comm,mpierr)
  if (mpierr.ne.0) then
     call errorlog('nubeam_nbinputs_bcast_step: character broadcast error')
!     call mpi_finalize(mpierr)
!     call exit(1)
     call bad_exit       
  endif

  mpierr=0
  call mpi_bcast(buflgc,nl_data,mpi_logical,0,nbi_comm,mpierr)
  if (mpierr.ne.0) then
     call errorlog('nubeam_nbinputs_bcast_step: logical broadcast error')
!     call mpi_finalize(mpierr)
!     call exit(1)
     call bad_exit       
  endif
#endif

  if (myidd.gt.0) then 
     call nbi_fld_bufget
  endif

  ni_data=0
  nmax_idata=0
  nr_data=0
  nmax_rdata=0
  nc_data=0
  nmax_cdata=0
  nl_data=0
  nmax_ldata=0

  if (allocated(bufint)) deallocate(bufint)
  if (allocated(bufr8))  deallocate(bufr8)
  if (allocated(bufchr)) deallocate(bufchr)
  if (allocated(buflgc)) deallocate(buflgc)

  !---------------------------------------------
  call infolog('nubeam_nbinputs_bcast_step: end')
  !---------------------------------------------

end subroutine nubeam_nbinputs_bcast_step 
!---------------------------------------------
!---------------------------------------------
!---------------------------------------------
subroutine nubeam_gatherv_chunks(nn_data)
  use nbi_com
  use nbi_mpi_com
  use nbi_pllbuf
  use logmod
  implicit none
  integer mpierr,i,j,l1,l2,l1_first
  integer, dimension(:,:), allocatable :: idispls,l_chunk
  integer, dimension(numprocd) :: nn_data !all cpu have this data as of mpi_allgathered
  integer :: l_displ,n_chunk_max,iertmp
  integer :: i_max_blk_size
  integer, dimension(:), allocatable :: n_chunk,n_remain
  integer :: nchunk,nremain
  real(fp) :: scpu,fcpu
  real(fp) :: t_gatherv=0.0_fp
!all myidd, find # of chunks to be transfer to root
  call cpu_time(scpu)

  if(allocated(n_chunk)) deallocate(n_chunk)
  if(allocated(n_remain)) deallocate(n_remain)
  allocate(n_chunk(numprocd),n_remain(numprocd),stat=iertmp)
	if (iertmp/=0) then 
	   call errorlog('nubeam_gatherv_chunks: allocat n_chunk stat=',iertmp)
           call errset_mpi(myidd,iertmp)
  	endif
  n_remain=0
  n_chunk=0

  do i=1,numprocd
!
! find block size, # of chunks for each processes
! 
  call nbi_chunks_find_fp(nn_data(i),i_max_blk_size,nchunk,nremain)

     n_chunk(i)=nchunk
     n_remain(i)=nremain

  enddo

  n_chunk_max=maxval(n_chunk)
  if(n_chunk_max.eq.0) then
     if(allocated(n_chunk)) deallocate(n_chunk)
     if(allocated(n_remain)) deallocate(n_remain)
     return
  endif

  if(allocated(l_chunk)) deallocate(l_chunk)
  allocate(l_chunk(n_chunk_max,numprocd),stat=iertmp)
  if (iertmp/=0) then 
     call errorlog('nubeam_gatherv_chunks: allocat l_chunk stat=',iertmp)
     call errset_mpi(myidd,iertmp)
  endif
  l_chunk=0
  do i=1,numprocd
     if(n_chunk(i).eq.0) cycle
     l_chunk(1:n_chunk(i),i)=i_max_blk_size
     if(n_remain(i).ne.0) then
        l_chunk(n_chunk(i),i)=abs(n_remain(i))
     endif
  enddo
  if(myidd.eq.0) then     
     if(allocated(r8buf1)) deallocate(r8buf1)
     allocate(r8buf1(sum(nn_data)),stat=iertmp)
     if (iertmp/=0) then 
        call errorlog('nubeam_gatherv_chunks: allocat r8buf1 stat=',iertmp)
        call errset_mpi(myidd,iertmp)
     endif
  endif
  if(allocated(idispls)) deallocate(idispls)
  allocate(idispls(n_chunk_max,numprocd),stat=iertmp)
  if (iertmp/=0) then 
     call errorlog('nubeam_gatherv_chunks: allocat idispls stat=',iertmp)
     call errset_mpi(myidd,iertmp)
  endif          
  idispls = 0 
  
  l_displ=0 !displacement
  do i=1,numprocd
     do j=1,n_chunk_max
        idispls(j,i)=l_displ
        l_displ=l_displ+l_chunk(j,i)
     enddo
     if((n_chunk(i).lt.n_chunk_max)) then
        if(n_chunk(i).ne.0) then
           idispls(n_chunk(i)+1:n_chunk_max,i)=idispls(n_chunk(i)+1:n_chunk_max,i)-1
        else
           idispls(:,i)=0
        endif
     endif
  enddo
     
  mpierr=0
  l1_first=1
  l2=0
  if(n_chunk(myidd+1).eq.0) then
     if(allocated(r8buf)) deallocate(r8buf)
     allocate(r8buf(1:1),stat=iertmp)
     if (iertmp/=0) then 
        call errorlog('nubeam_gatherv_chunks: allocat r8buf stat=',iertmp)
        call errset_mpi(myidd,iertmp)
     endif
     r8buf=0.0_fp
  endif
#ifdef __MPI
  call infolog('gatherv n_chunk_max=',n_chunk_max)
  do i=1,n_chunk_max
     call infolog('gatherv i chunk=',i)
     l1=l2+l1_first
     l2=l2+l_chunk(i,myidd+1)
     if(i.le.n_chunk(myidd+1)) then
        call mpi_gatherv(r8buf(l1:l2),l_chunk(i,myidd+1),MPI_REAL8,r8buf1,nn_data,idispls(i,1:numprocd),&
             & MPI_REAL8,0,nbi_comm,mpierr)
     else
        if(n_chunk(myidd+1).eq.0) then
           call mpi_gatherv(r8buf,0,MPI_REAL8,r8buf1,nn_data,idispls(i,1:numprocd),&
                & MPI_REAL8,0,nbi_comm,mpierr)
        else
           call mpi_gatherv(r8buf(l2:l2),1,MPI_REAL8,r8buf1,nn_data,idispls(i,1:numprocd),&
                & MPI_REAL8,0,nbi_comm,mpierr)
        endif
     endif
  enddo
#endif
  if(allocated(n_chunk)) deallocate(n_chunk)
  if(allocated(n_remain)) deallocate(n_remain)
  if(allocated(l_chunk)) deallocate(l_chunk)
  if(allocated(idispls)) deallocate(idispls)

  call cpu_time(fcpu)
  t_gatherv=t_gatherv+(fcpu-scpu)
  call infolog('nubeam_gatherv_chunks: gatherv cpu time = ',t_gatherv)
end subroutine nubeam_gatherv_chunks
!---------------------------------------------
!---------------------------------------------
!---------------------------------------------
subroutine nubeam_bcast_chunks_int(name_blk,len_blk,istat)
  use nbi_mpi_com
  use nbi_com, only : blk_mpi_bcast_int, blk_mpi_bcast_r8,fp
  use logmod
  implicit none
  integer,intent(in) :: len_blk
  integer,dimension(1:len_blk) :: name_blk
  integer, intent(in) :: istat  ! error code (user defined), .gt.0 expected...
! local
  integer :: i_max_blk_size
  integer :: n_chunk,n_remain
  integer,dimension(:), allocatable :: l_chunk
  integer :: mpierr,i,l1,l2,l1_first,iertmp
  real(fp) :: scpu,fcpu
  real(fp) :: t_bcast_int=0.0_fp

  call cpu_time(scpu)
!
! find block size, # of chunks 
! 
  call nbi_chunks_find_int(len_blk,i_max_blk_size,n_chunk,n_remain)
  if(n_chunk.eq.0) then
     return
  endif

  call infolog('nubeam_bcast_chunks_int, n_chunk=',n_chunk)

  if(allocated(l_chunk)) deallocate(l_chunk)
  allocate(l_chunk(n_chunk))
  l_chunk=i_max_blk_size

  if(n_remain.ne.0) then
     l_chunk(n_chunk)=abs(n_remain)
  endif

#ifdef __MPI

  l1_first=1
  l2=0
  do i=1,n_chunk
     l1=l2+l1_first
     l2=l2+l_chunk(i)
     mpierr=0

     call infolog('nubeam_bcast_chunks_int, l_chunk(i)=',l_chunk(i))

    call mpi_bcast(name_blk(l1:l2),l_chunk(i),MPI_INTEGER,0,nbi_comm,mpierr)

     if (mpierr.ne.0) then
        call errorlog('nubeam_bcast_chunks_int: broadcast error')
        call errset_mpi(myidd,istat)
     endif
     
  enddo
  
#endif
  
  if(allocated(l_chunk)) deallocate(l_chunk)
  call cpu_time(fcpu)
  t_bcast_int=t_bcast_int+(fcpu-scpu)
  call infolog('nubeam_bcast_chunks_int: bcast int cpu time = ',t_bcast_int)


end subroutine nubeam_bcast_chunks_int
!---------------------------------------------
!---------------------------------------------
!---------------------------------------------
subroutine nubeam_bcast_chunks_r8(name_blk,len_blk,istat)
  use nbi_mpi_com
  use nbi_com, only : blk_mpi_bcast_int, blk_mpi_bcast_r8,fp
  use logmod
  implicit none
  integer,intent(in) :: len_blk
  real(fp),dimension(1:len_blk) :: name_blk
  integer, intent(in) :: istat  ! error code (user defined), .gt.0 expected...
! local
  integer :: i_max_blk_size
  integer :: n_chunk,n_remain
  integer,dimension(:), allocatable :: l_chunk
  integer :: mpierr,i,l1,l2,l1_first,iertmp
  real(fp) :: scpu,fcpu
  real(fp) :: t_bcast_fp=0.0_fp
  call cpu_time(scpu)
!
! find block size, # of chunks 
! 
  call nbi_chunks_find_fp(len_blk,i_max_blk_size,n_chunk,n_remain)

  if(n_chunk.eq.0) then
     return
  endif

  call infolog('nubeam_bcast_chunks_r8, n_chunk=',n_chunk)

  if(allocated(l_chunk)) deallocate(l_chunk)
  allocate(l_chunk(n_chunk))
  l_chunk=i_max_blk_size

  if(n_remain.ne.0) then
     l_chunk(n_chunk)=abs(n_remain)
  endif


#ifdef __MPI

  l1_first=1
  l2=0

  do i=1,n_chunk
     l1=l2+l1_first
     l2=l2+l_chunk(i)
     mpierr=0

     call infolog('nubeam_bcast_chunks_r8, l_chunk(i)=',l_chunk(i))

     call mpi_bcast(name_blk(l1:l2),l_chunk(i),MPI_REAL8,0,nbi_comm,mpierr)
     
     if (mpierr.ne.0) then
        call errorlog('nubeam_bcast_chunks_r8: broadcast error')
        call errset_mpi(myidd,istat)
     endif
     
  enddo

#endif

  if(allocated(l_chunk)) deallocate(l_chunk)
  call cpu_time(fcpu)
  t_bcast_fp=t_bcast_fp+(fcpu-scpu)
  call infolog('nubeam_bcast_chunks_r8: bcast r8 cpu time = ',t_bcast_fp)
 

end subroutine nubeam_bcast_chunks_r8
!---------------------------------------------
!---------------------------------------------
!---------------------------------------------
subroutine nubeam_reduce_chunks_r8(name_blk,name_blk_recv,len_blk,mpi_op,istat)
!only for myid = 0
  !.. nubeam_reduce_chunks_r8          <description>
  !
  !  modifications
  !
  !  1.00  06jan2010                   jim.conboy@ccfe.ac.uk
  !                                    fix mpi_op test ; 
  !                                    cant call uupper on mpi_op ( its input )
  !                                    cant trim the arg itself ( fails if debug compilation )
  !----_^---------------------------------===
  use nbi_mpi_com
  use nbi_com, only : blk_mpi_bcast_int, blk_mpi_bcast_r8,fp
  use logmod
  implicit none
  integer,intent(in) :: len_blk
  real(fp),intent(in) :: name_blk(len_blk)
  real(fp),intent(out) :: name_blk_recv(len_blk)
  character(len=*),intent(in) :: mpi_op !reduce operation (handle) 
  integer, intent(in) :: istat  ! error code (user defined), .gt.0 expected...
! local
  integer                           :: i_max_blk_size
  integer                           :: n_chunk,n_remain
  integer,dimension(:), allocatable :: l_chunk
  integer                           :: mpierr,i,l1,l2,l1_first,iertmp
  real(fp)                            :: scpu,fcpu
  real(fp)                            :: t_reduce_fp=0.0_fp
  character(len=32)                 :: mpi_op_test
! -------------------------------------------------
  call cpu_time(scpu)
!
! check eligibility
!
!-call uupper(trim(mpi_op))  jec 22dec2009  fails under debug compilation
  mpi_op_test = mpi_op //' '
  call uupper( mpi_op_test )
  if( mpi_op_test(:8) .ne. 'MPI_SUM ' )    then
     call errorlog('nubeam_reduce_chunks_r8: reduce operation is eligible only for  MPI_SUM')
     call errorlog('nubeam_reduce_chunks_r8: reduce operation = '//trim(mpi_op))
     print *,'mpi_op_test :- "',mpi_op_test(:8),'"'
     call exit(1)
  endif
!
! find block size, # of chunks 
! 
  call nbi_chunks_find_fp(len_blk,i_max_blk_size,n_chunk,n_remain)

  if(n_chunk.eq.0) then
     return
  endif

  call infolog('nubeam_reduce_chunks_r8, n_chunk=',n_chunk)

  if(allocated(l_chunk)) deallocate(l_chunk)
  allocate(l_chunk(n_chunk))
  l_chunk=i_max_blk_size

  if(n_remain.ne.0) then
     l_chunk(n_chunk)=abs(n_remain)
  endif


#ifdef __MPI


  l1_first=1
  l2=0

  do i=1,n_chunk
     l1=l2+l1_first
     l2=l2+l_chunk(i)
     mpierr=0

     call infolog('nubeam_reduce_chunks_r8, l_chunk(i)=',l_chunk(i))

     call mpi_reduce(name_blk(l1:l2),name_blk_recv(l1:l2),l_chunk(i),MPI_REAL8,MPI_SUM,0,nbi_comm,mpierr)
     
     if (mpierr.ne.0) then
        call errorlog('nubeam_reduce_chunks_r8: broadcast error, chunk# = ',1)
        call errset_mpi(myidd,istat)
     endif
     
  enddo

#endif

  if(allocated(l_chunk)) deallocate(l_chunk)
  call cpu_time(fcpu)
  t_reduce_fp=t_reduce_fp+(fcpu-scpu)
  call infolog('nubeam_reduce_chunks_r8: reduce r8 cpu time = ',t_reduce_fp)

end subroutine nubeam_reduce_chunks_r8
!---------------------------------------------
!---------------------------------------------
subroutine nubeam_reduce1_chunks_r8(name_blk,name_blk_recv,len_blk,mpi_op,istat)
!only for myid > 0
  !.. nubeam_reduce_chunks_r8          <description>
  !
  !  modifications
  !
  !  1.00  06jan2010                   jim.conboy@ccfe.ac.uk
  !                                    fix mpi_op test ; 
  !                                    cant call uupper on mpi_op ( its input )
  !                                    cant trim the arg itself ( fails if debug compilation )
  !----_^---------------------------------===
  use nbi_mpi_com
  use nbi_com, only : blk_mpi_bcast_int, blk_mpi_bcast_r8,fp
  use logmod
  implicit none
  integer,intent(in) :: len_blk
  real(fp),intent(in) :: name_blk(len_blk)
  real(fp),intent(out) :: name_blk_recv(1)
  character(len=*),intent(in) :: mpi_op !reduce operation (handle) 
  integer, intent(in) :: istat  ! error code (user defined), .gt.0 expected...
! local
  integer                           :: i_max_blk_size
  integer                           :: n_chunk,n_remain
  integer,dimension(:), allocatable :: l_chunk
  integer                           :: mpierr,i,l1,l2,l1_first,iertmp
  real(fp)                            :: scpu,fcpu
  real(fp)                            :: t_reduce_fp=0.0_fp
  character(len=32)                 :: mpi_op_test
! -------------------------------------------------
  call cpu_time(scpu)
!
! check eligibility
!
!-call uupper(trim(mpi_op))  jec 22dec2009  fails under debug compilation
  mpi_op_test = mpi_op //' '
  call uupper( mpi_op_test )
  if( mpi_op_test(:8) .ne. 'MPI_SUM ' )    then
     call errorlog('nubeam_reduce_chunks_r8: reduce operation is eligible only for  MPI_SUM')
     call errorlog('nubeam_reduce_chunks_r8: reduce operation = '//trim(mpi_op))
     print *,'mpi_op_test :- "',mpi_op_test(:8),'"'
     call exit(1)
  endif
!
! find block size, # of chunks 
! 
  call nbi_chunks_find_fp(len_blk,i_max_blk_size,n_chunk,n_remain)

  if(n_chunk.eq.0) then
     return
  endif

  call infolog('nubeam_reduce_chunks_r8, n_chunk=',n_chunk)

  if(allocated(l_chunk)) deallocate(l_chunk)
  allocate(l_chunk(n_chunk))
  l_chunk=i_max_blk_size

  if(n_remain.ne.0) then
     l_chunk(n_chunk)=abs(n_remain)
  endif


#ifdef __MPI


  l1_first=1
  l2=0

  do i=1,n_chunk
     l1=l2+l1_first
     l2=l2+l_chunk(i)
     mpierr=0

     call infolog('nubeam_reduce_chunks_r8, l_chunk(i)=',l_chunk(i))

     call mpi_reduce(name_blk(l1:l2),name_blk_recv,l_chunk(i),MPI_REAL8,MPI_SUM,0,nbi_comm,mpierr)
     
     if (mpierr.ne.0) then
        call errorlog('nubeam_reduce_chunks_r8: broadcast error, chunk# = ',1)
        call errset_mpi(myidd,istat)
     endif
     
  enddo

#endif

  if(allocated(l_chunk)) deallocate(l_chunk)
  call cpu_time(fcpu)
  t_reduce_fp=t_reduce_fp+(fcpu-scpu)
  call infolog('nubeam_reduce_chunks_r8: reduce r8 cpu time = ',t_reduce_fp)

end subroutine nubeam_reduce1_chunks_r8
!---------------------------------------------
!---------------------------------------------
!---------------------------------------------
subroutine nubeam_allreduce_chunks_r8(name_blk,name_blk_recv,len_blk,mpi_op,istat)
  !.. nubeam_allreduce_chunks_r8          <description>
  !
  !  modifications
  !
  !  1.00  06jan2010                   jim.conboy@ccfe.ac.uk
  !                                    fix mpi_op test ; 
  !                                    cant call uupper on mpi_op ( its input )
  !                                    cant trim the arg itself ( fails if debug compilation )
  !----_^---------------------------------===
  use nbi_mpi_com
  use nbi_com, only : blk_mpi_bcast_int, blk_mpi_bcast_r8,fp
  use logmod
  implicit none
  integer,intent(in) :: len_blk
  real(fp),intent(in) :: name_blk(len_blk)
  real(fp),intent(out) :: name_blk_recv(len_blk)
  character(len=*),intent(in) :: mpi_op !reduce operation (handle) 
  integer, intent(in) :: istat  ! error code (user defined), .gt.0 expected...
! local
  integer :: n_chunk,n_remain
  integer :: i_max_blk_size
  integer,dimension(:), allocatable :: l_chunk
  integer :: mpierr,i,l1,l2,l1_first,iertmp
  real(fp) :: scpu,fcpu
  real(fp) :: t_allreduce_fp=0.0_fp
  character(len=32)                 :: mpi_op_test
! -------------------------------------------------
  call cpu_time(scpu)
!
! check eligibility
!
!-call uupper(trim(mpi_op))  jec 22dec2009  fails under debug compilation
  mpi_op_test = mpi_op //' '
  call uupper( mpi_op_test )
  if( mpi_op_test(:8) .ne. 'MPI_SUM ' )    then
     call errorlog('nubeam_allreduce_chunks_r8: reduce operation is eligible only for  MPI_SUM')
     call errorlog('nubeam_allreduce_chunks_r8: reduce operation = '//trim(mpi_op))
     call exit(1)
  endif
!
! find block size, # of chunks 
! 
  call nbi_chunks_find_fp(len_blk,i_max_blk_size,n_chunk,n_remain)

  if(n_chunk.eq.0) then
     return
  endif

  call infolog('nubeam_allreduce_chunks_r8, n_chunk=',n_chunk)
!
!find chunk length
!
  if(allocated(l_chunk)) deallocate(l_chunk)
  allocate(l_chunk(n_chunk))
  l_chunk=i_max_blk_size

  if(n_remain.ne.0) then
     l_chunk(n_chunk)=abs(n_remain)
  endif


#ifdef __MPI


  l1_first=1
  l2=0

  do i=1,n_chunk

     l1=l2+l1_first
     l2=l2+l_chunk(i)
     mpierr=0

     call infolog('nubeam_allreduce_chunks_r8, l_chunk(i)=',l_chunk(i))

     call mpi_allreduce(name_blk(l1:l2),name_blk_recv(l1:l2),l_chunk(i),MPI_REAL8,MPI_SUM,nbi_comm,mpierr)
     
     if (mpierr.ne.0) then
        call errorlog('nubeam_allreduce_chunks_r8: broadcast error, chunk# = ',1)
        call errset_mpi(myidd,istat)
     endif
     
  enddo

#endif

  if(allocated(l_chunk)) deallocate(l_chunk)
  call cpu_time(fcpu)
  t_allreduce_fp=t_allreduce_fp+(fcpu-scpu)
  call infolog('nubeam_allreduce_chunks_r8: allreduce r8 cpu time = ',t_allreduce_fp)

end subroutine nubeam_allreduce_chunks_r8
!---------------------------------------------
!---------------------------------------------
!---------------------------------------------
subroutine nbi_chunks_find_int(len_blk,i_max_blk_size,n_chunk,n_remain)
  use nbi_mpi_com
  use nbi_com, only :  blk_mpi_bcast_int
  use logmod
  implicit none
  integer,intent(in) :: len_blk
  integer,intent(out) :: n_chunk,n_remain,i_max_blk_size
  integer :: i_max_blk_size_int = 0
  integer :: atoi,iertmp,mpierr
  character*200 chbuf
!
! get environmental variable max_mpi_int_block, if not set use default value
!
  n_chunk=0
  n_remain=0
  i_max_blk_size=0

  if(len_blk.eq.0) then
     call infolog('nbi_chunks_find_int: len_blk = ',len_blk)
     call infolog('nbi_chunks_find_int: no data available')
     return
  endif

  if(i_max_blk_size_int.eq.0) then
        
     call mpi_sget_env('max_mpi_int_block',chbuf,iertmp)
     call nbi_ckerr('nbi_chunks_find_int: mpi_sget_env("max_mpi_int_block")',iertmp)
     
! if env. variable max_mpi_int_block was set it will be set same on all nodes
! should not be a problem further, for example when master has no such variable 
! and will try to broadcast to slaves, which have this one set and did not expect 
! to get any broadcasted value
!
     if(trim(chbuf).eq.' ') then

#ifdef __MPI
        mpierr=0
        call mpi_bcast(blk_mpi_bcast_int,1,MPI_INTEGER,0,nbi_comm,mpierr)
        if (mpierr.ne.0) then
           call errorlog('nbi_chunks_find_int: blk_mpi_bcast_int broadcast error')
!           call mpi_finalize(mpierr)
!           call exit(1)
           call bad_exit       
        endif
#endif
        call warnlog(' nbi_chunks_find_int:'// &
             ' max_mpi_int_block environment variable: undefined.'// &
             ' use default  = ',blk_mpi_bcast_int)

        i_max_blk_size_int=blk_mpi_bcast_int/4 !default converted from  bytes to the number of elements in array 

     else

        i_max_blk_size_int=atoi(chbuf,1)/4 !convert bytes to the number of elements in array

        if(i_max_blk_size_int.le.0) then !use default value

           i_max_blk_size_int=blk_mpi_bcast_int/8 !default, convert bytes to the number of elements in array

           call warnlog(' nbi_chunks_find_int:'// &
                ' max_mpi_int_block environment variable: badly defined.'// &
                ' use default = ',blk_mpi_bcast_int)
        endif

     endif
     
  endif

  i_max_blk_size=i_max_blk_size_int

  if (i_max_blk_size.le.0) then
     call errorlog(' nbi_chunks_find_int: bad definition for the  block size, i_max_blk_size =  ',i_max_blk_size)
     call exit(1)
  endif

!
!find number of chunks to be broadcast 
!
  n_chunk=len_blk/i_max_blk_size
  n_remain=mod(len_blk,i_max_blk_size)
  if(n_remain.ne.0) then
     n_chunk=n_chunk+1
  endif

end subroutine nbi_chunks_find_int
!---------------------------------------------
!---------------------------------------------
!---------------------------------------------
subroutine nbi_chunks_find_fp(len_blk,i_max_blk_size,n_chunk,n_remain)
  use nbi_mpi_com
  use nbi_com, only :  blk_mpi_bcast_r8
  use logmod
  implicit none
  integer,intent(in) :: len_blk
  integer,intent(out) :: n_chunk,n_remain,i_max_blk_size
  integer :: i_max_blk_size_fp = 0
  integer :: atoi,iertmp,mpierr
  character*200 chbuf
!
! get environmental variable max_mpi_fp_block, if not set use default valie
!
  n_chunk=0
  n_remain=0
  i_max_blk_size=0
  
  if(len_blk.eq.0) then
     call infolog('nbi_chunks_find_fp: len_blk = ',len_blk)
     call infolog('nbi_chunks_find_fp: no data available')
     return
  endif

  if(i_max_blk_size_fp.eq.0) then
        
     call mpi_sget_env('max_mpi_fp_block',chbuf,iertmp)
    
     call nbi_ckerr('nbi_chunks_find_fp: mpi_sget_env("max_mpi_fp_block")',iertmp)
     
! if env. variable max_mpi_fp_block was set it will be set same on all nodes
! should not be a problem further, for example when master has no such variable 
! and will try to broadcast to slaves, which have this one set and did not expect 
! to get any broadcasted value
!
     if(trim(chbuf).eq.' ') then

#ifdef __MPI
        mpierr=0
        call mpi_bcast(blk_mpi_bcast_r8,1,MPI_INTEGER,0,nbi_comm,mpierr)
        if (mpierr.ne.0) then
           call errorlog('nbi_chunks_find_fp: blk_mpi_bcast_r8 broadcast error')
!           call mpi_finalize(mpierr)
!           call exit(1)
           call bad_exit       
        endif
#endif
        call warnlog(' nbi_chunks_find_fp:'// &
             ' max_mpi_fp_block environment variable: undefined.'// &
             ' use default  = ',blk_mpi_bcast_r8)

        i_max_blk_size_fp=blk_mpi_bcast_r8/8 !default, convert bytes to the number of elements in array

     else

        i_max_blk_size_fp=atoi(chbuf,1)/8 !convert bytes to the number of elements in array

        if(i_max_blk_size_fp.le.0) then !use default value

           i_max_blk_size_fp=blk_mpi_bcast_r8/8 !default, convert bytes to the number of elements in array

           call warnlog(' nbi_chunks_find_fp:'// &
                ' max_mpi_fp_block environment variable: badly defined.'// &
                ' use default = ',blk_mpi_bcast_r8)
        endif
        
     endif
     
  endif

  i_max_blk_size=i_max_blk_size_fp

  if (i_max_blk_size.le.0) then
     call errorlog(' nbi_chunks_find_fp: bad definition for the  block size, i_max_blk_size =  ',i_max_blk_size)
     call exit(1)
  endif
!
!find number of chunks to be broadcast 
!
  n_chunk=len_blk/i_max_blk_size
  n_remain=mod(len_blk,i_max_blk_size)
  if(n_remain.ne.0) then
     n_chunk=n_chunk+1
  endif

  
end subroutine nbi_chunks_find_fp
!---------------------------------------------
!---------------------------------------------
!---------------------------------------------
subroutine nbi_ckerr(str,iertmp)
  use logmod
  implicit none
  character*(*), intent(in) :: str
  integer :: iertmp
  
  if(iertmp.ne.0) then
     call errorlog( 'error code detected on call: '//trim(str))
  endif
end subroutine nbi_ckerr
!---------------------------------------------
!---------------------------------------------
!---------------------------------------------
subroutine nbi_buf_netcdf(cact,zfile_cdf,ndata,rbuf)
  use iso_c_binding, only: fp => c_double
  use netcdf
  character*(*) :: cact
  character*(*) :: zfile_cdf
  integer :: ndata
  real(fp), dimension(*) :: rbuf
  integer, save :: icdf
  integer :: id_ndata,id_rbuf

  call uupper(trim(adjustl(cact)))
 
  if(cact.eq.'W') then

     call check(nf90_create(path=trim(adjustl(zfile_cdf)),cmode=nf90_clobber,ncid=icdf),'create file')
     call check(nf90_def_dim(icdf,'ndata',ndata,id_ndata),'ndata')
     call check(nf90_def_var(icdf,'r_buf',nf90_double,(/id_ndata/),id_rbuf),'def r_buf')
     call check(nf90_enddef(icdf),'end_def') ! out of def mode and go to data mode
     call check(nf90_put_var(icdf,id_rbuf,rbuf(1:ndata)),'put r_buf')
     call check(nf90_close(icdf),'close file')

  elseif(cact.eq.'D') then
     
     call check(nf90_open(trim(adjustl(zfile_cdf)),nf90_nowrite,icdf),'create file')
     call check(nf90_inq_dimid(icdf,'ndata',id_ndata),'id ndata')
     call check(nf90_inquire_dimension(icdf,id_ndata,len=ndata) ,'get ndata')

  elseif(cact.eq.'R') then
     
     call check(nf90_inq_varid(icdf,'r_buf',id_rbuf),'id r_buf')
     call check(nf90_get_var(icdf,id_rbuf,rbuf(1:ndata)),'get rbuf')
     call check(nf90_close(icdf),'close file')

  endif
  

contains

  subroutine check(status,varnam)
    integer, intent ( in) :: status 
    character*(*),intent(in) :: varnam ! dimension name input 
    
    if(status /= nf90_noerr) then 
       write(nonlin,*) ' ?nbi_buf_netcdf error:'//trim(nf90_strerror(status))//' for var.= '//trim(varnam) 
       call exit(status) 
    end if
  end subroutine check
  

end subroutine nbi_buf_netcdf
