subroutine nubeam_step(ss,zcptime,ierr)

  ! nubeam component -- advance nubeam time step
  ! assumption:
  !   nubeam responsible for deposition and slowing down of 
  !   neutral beam and fusion product fast ion species

  use nbi_com
  use nbi_mpi_com
  use xplasma_obj_instance
  use plasma_state_mod
  use logmod

  implicit none

  !-----------------------------------
  ! passed:

  type (plasma_state)   :: ss

  real(fp), intent(inout) :: zcptime ! cpu time in setup routine (caller)
  !  time for inserting output into plasma state instance (ss) is added.

  integer, intent(out)  :: ierr    ! status code, 0=ok

  !-----------------------------------
  ! local:

  integer :: isb,ibeami,ifusi,ib,inb
  real(fp) :: ztim1,ztim2
  real(fp) :: zcptim1,zcptim2,zcptim_all
  !--------------------------------------------------------------------
  logical :: ifran_wall ! frantic post-processing => output plasma state:
  !            .true. to compute wall sources.
  logical :: ifran_halo ! frantic post-processing => output plasma state:
  !            .true. to compute halo sources.
  logical :: ifran_reco ! frantic post-processing => output plasma state:
  !            .true. to compute recombination sources.
  character*200 :: wrkpath


  call cpu_time(zcptim1)
  ! shared code, all processors:
  write(nonlin,*)' %nubeam_step ++++++++ nltrk_dep0',nltrk_dep0
  write(nonlin,*)' %nubeam_step ++++++++ nlbout (lost ptcl.)',nlbout
  write(nonlin,*)' %nubeam_step ++++++++ nbsteps ',nbsteps
  write(nonlin,*)' %nubeam_step ++++++++ ntest_max ',ntest_max
  call fort_flush(nonlin)
  call nubeam_step_child(ierr) 
  if(ierr.ne.0) return

  !-----------------------------------
  if (myidd.eq.0) then
     ! mpi cpu #0 only...
     call cpu_time(ztim1)

     write(nonlin,*) ' nubeam_load_state... '
     call fort_flush(nonlin)
     ibeami=0
     ifusi=0
     do isb=1,nsbeam
        if(nlfprod(isb)) then
           ifusi=ifusi+1
        else
           ibeami=ibeami+1
        endif
     enddo

     call nubeam_chk_state_grids(ss,ierr)

     if(ierr.eq.0) then
        ! copy current fractions which may not be reflected in the output state

        if(ibeami.gt.0) then
           inb=0
           do ib=1,nbeam
              if(ntrace(ib).eq.0) then
                 inb=inb+1
                 if(inb.gt.ss%nbeam) cycle
                 ss%frac_full(inb) = ffulla(ib)
                 ss%frac_half(inb) = fhalfa(ib)
              endif
           enddo
           ! sanity check
           if(inb.ne.ss%nbeam) then
              ierr=77
              call nbi_warning(' ??nubeam_step: (#beams).ne.ss%nbeam, output state dimension error!')
           endif
        endif
     endif

     if(ierr.eq.0) then
        call infolog('nubeam_step:nubeam_load_state')
        call nubeam_load_state(ss,(ibeami.gt.0),(ifusi.gt.0),ierr)
     endif


     call cpu_time(ztim2)
     zcptime = zcptime + (ztim2-ztim1)

     !-----------------------------------
     !  store ps i/o subroutines' cpu time now
     !  this makes it available via normal nubeam output scalars access

     cpbps = zcptime


     write(nonlin,*) ' nubeam_step: levmod_halo= ',levmod_halo
     
     if((levmod_halo.eq.1.or.levmod_halo.eq.2).and.ierr.eq.0) then
        !run  frantic after plasma state has been updated
        ! see nubeam_comp_exec l 1420
        ! ...
        ! define 
        write(nonlin,*) ' frantic... '
        ifran_wall=.false.
        ifran_halo=.true.
        ifran_reco=.false.
        wrkpath=trim(workpath)
        call infolog('nubeam_step:frantic start')
        call nubeam_frantic(ss,nonlin, &
             ifran_wall,ifran_halo,ifran_reco, &
             wrkpath,ierr)
        if(ierr.ne.0) then
           call errorlog(' nubeam_step: frantic failed')
        else
  ! load frantic stuff into xplasma
           
           call infolog('nubeam_step:frantic done')
           call xplasma_author_set(s,'nubeam',ierr)
           if(ierr.ne.0) go to 999
           call nbo_load_frantic
           call xplasma_author_clear(s,'nubeam',ierr)
        endif
        
     endif
     
  endif

  !-----------------------------------

999 continue
  call cpu_time(zcptim2)
  zcptim_all = zcptim2 - zcptim1

end subroutine nubeam_step

!======================================================

subroutine nubeam_step_child(ierr)

  ! this routine contains all the mpi parts of nubeam
  ! but has a simpler argument list than nubeam_step: just a status code return
  !   (myidd.gt.0) threads can call this directly w/o going through nubeam_step

  use nbi_com
  use nbi_mpi_com
  use logmod
  implicit none

  !-----------------------------------
  ! passed:

  integer, intent(out)  :: ierr    ! status code, 0=ok

  !-----------------------------------
  ! local:

  integer :: iorbtot,j
  real(fp) :: zcptim1,zcptim2,zcptim_loc
  character*180 :: suffix_ptcl ! nubeam's ptcl state suffix
  integer myid_debug, numproc_debug
  character*180 logfile_level, log_file ! log file level and name
  integer istat,ilen,nbsteps_prv,isb
  character*250 :: f_name
  !--------------------------------------------------------------------

  call enterlog("nubeam_step_child") ! all returns through exitlog()

  call cpu_time(zcptim1)
#ifndef __MPI
      if(debug_mpi.and.nltest_output) then
         myid_debug=myidd
         numproc_debug=numprocd
      endif
#endif

      call nbi_mpi_com_init

#ifndef __MPI
      if(debug_mpi.and.nltest_output) then
         myidd=myid_debug       ! do not redefine  myidd and nuprocd for 
         numprocd=numproc_debug ! debugging of mpi version by serial code 
      endif
#endif

  write(nonlin,*) ' ==> entering nubeam_step_child, mpi myidd = ',myidd
  call fort_flush(nonlin)

#ifdef __MPI
  call infolog('call nubeam_xplasma_bcast: numprocd=',numprocd)
  call nubeam_xplasma_bcast   

  if (.not.nbi_init_flag) then
     call infolog('call nubeam_nbinputs_bcast_init: numprocd=',numprocd)
     call nubeam_nbinputs_bcast_init
  else
     call infolog('call nubeam_nbinputs_bcast_step: numprocd',numprocd)
     call nubeam_nbinputs_bcast_step
  endif
#endif

  if(nbi_init_flag) then  
     
     if(need_read_ptcl) then !need to read ptcl_state data 
        ! fld state has been initialized  by nbi_fld_state
        ! in resume for nubeam in transp or in
        ! nubeam_copm_exec for standalone nubeam.
        ! ptcl.state file need to be read
        suffix_ptcl='_nbi_ptcl_state_cpu' ! suffix for ptcl list file
        ! need to be consistent with
        ! suffix in nubeam_comp_exe/nubeam_comp_exec.f90
        ! and outcor/wrstrt.for
        ! and trcore/resume.for
        ! and trmpi/trmpi_listener.f90
        ! and nubeam/nubeam_step.f90

        call infolog('call nbi_ptcl_read:numprocd',numprocd)
        if(abs(ntest_max).eq.1) then ! run bi-maxwellian distribution for fast ion deposition
! for this test all ptcl will be deposited from scratch
           nbsteps=0
           puseyp=0.0_fp
           nlbflr=.false.
           nlfbmflr=.false.
           orbrzv_zzerr_con=1.0e-6_fp !set high tolerence for rz integrator, important
        endif
        call nbi_ptcl_read(suffix_ptcl)
        need_read_ptcl = .false. ! myidd do not read these file during next step
     endif
  endif

        

  !--------------------------------------------------------------------------
  !  if(myidd.eq.0) 
  write(nonlin,*) ' nbi_interp_profiles... '
  call fort_flush(nonlin)
  call infolog('call nbi_interp_profiles')
  call nbi_interp_profiles(ierr)
  if(ierr.ne.0) goto 900

  !-----------------------------------
!  if(myidd.eq.0) 
  write(nonlin,*) ' nbstart... '
  call fort_flush(nonlin)
  call infolog('call nbstart')
!  call nubeam_workpath_reset(workpath)
  call nbstart(ierr)
  if(ierr.ne.0) goto 900

  !-----------------------------------
!  if(myidd.eq.0) 
  write(nonlin,*) ' depall... '
  if(abs(ntest_max).eq.1) then
     write(nonlin,*) ' bi_maxwellian test... '
     write(nonlin,*) ' bi_max_pr =',bi_max_pr
     write(nonlin,*) ' bi_max_pl =',bi_max_pl
  endif
  call fort_flush(nonlin)
  call infolog('call depall')
  call depall(ierr)
  write(nonlin,*) ' depall exited ',ierr
  call fort_flush(nonlin)
!c save nubeam state, for debug
  if(abs(ntest_max).eq.1) then
! for this test save ptcl list as deposited.
     ilen=len(trim(workpath))
     nbsteps_prv=nbsteps ! save
     nbsteps=1 !to allow write ptcl list to the disk
     ebavps=0.0_fp ! commented  to save average beam energy used to 
                      ! simulate this bi-maxwellian distribution
     ebapl2s=0.0_fp
     ebapp2s=0.0_fp
     do isb=1,nsbeam
        lsbeam=isb
        if(wghtsave(isb).gt.0.0_fp) then
#ifdef __MPI
           call nbi_ebav_reduce
#endif
           wghtay(:,isb)=wghtay(:,isb)*wghtsave(isb)
           do j=lcentr,ledge
              if(bdenss(j,isb).gt.0.0_fp) then
                 ebapl2s(j,isb)=ebapl2sw(j,isb)/bdenss(j,isb) !ev
                 ebapp2s(j,isb)=ebapp2sw(j,isb)/bdenss(j,isb) !ev
                 ebavps(j,isb)=ebapl2s(j,isb)+ebapp2s(j,isb) !ev
                 bdenss(j,isb)=bdenss(j,isb)*wghtsave(isb)/delt ! beam density #, consistant with outptb.for
              endif
           enddo
        endif
#ifdef __DEBUG
        call nubeam_maxwell_plot(.true.)
#endif
     enddo     

!#ifdef __DEBUG

     if (myidd.eq.0) then
        f_name = trim(workpath)//trim(runid)//'_max.cdf'
        call nubeam2toric_wr(f_name)
     endif
     call nbi_states(workpath(1:ilen)//trim(runid(1:lrunid))//'_max_nbi_fld_state.cdf',&
          '_max_nbi_ptcl_state_cpu','s',0,ierr) ! see nubeam_test

!#endif
     nbsteps=nbsteps_prv !back to previous one
     do isb=1,nsbeam
        if(wghtsave(isb).gt.0.0_fp) then
           wghtay(:,isb)=wghtay(:,isb)/wghtsave(isb)
        endif
     enddo
     bdenss=0.0_fp
     ebavps=0.0_fp
     ebapl2sw=0.0_fp
     ebapp2sw=0.0_fp
     ebapl2s=0.0_fp
     ebapp2s=0.0_fp

  endif
!m  ierr=1


!mg
!mg  goto 900
  if(ierr.ne.0) goto 900

  !-----------------------------------
!  if(myidd.eq.0) 
  write(nonlin,*) ' orball... '
  call fort_flush(nonlin)
  call infolog('call orball')
  call orball(iorbtot,ierr)
#ifdef __DEBUG
  if(abs(ntest_max).eq.1) then
     do isb=1,nsbeam
        lsbeam=isb
        call nubeam_maxwell_plot(.false.)
     enddo
  endif
#endif
  if ((myidd.eq.0).and.nltest_output) then
     if(abs(ntest_max).eq.1) then
        f_name = trim(workpath)//trim(runid)//'_max_after_orb.cdf'
     elseif(ntest_max.eq.0.or.abs(ntest_max).eq.10) then
        f_name = trim(workpath)//trim(runid)//'_nonmax_after_orb.cdf'
     endif
     if(ntest_max.eq.-1.or.abs(ntest_max).eq.10) then
        call nubeam2toric_wr(f_name)
     endif
  endif


  write(nonlin,*) ' orball exited ',ierr
  call fort_flush(nonlin)
  if(ierr.ne.0) goto 900

  !-----------------------------------
  if (myidd.eq.0) then
     write(nonlin,*) ' nbfinish... '
     call fort_flush(nonlin)
     call infolog('call nbfinish')
     call nbfinish(ierr)
     if(ierr.ne.0) goto 900
  endif

  write(nonlin,*) ' ==> exiting nubeam_step_child, mpi myidd = ',myidd
  call fort_flush(nonlin)

  call cpu_time(zcptim2)
  zcptim_loc = zcptim2 - zcptim1
  
  if(ierr.eq.0) then
     nbsteps = nbsteps + 1
  endif

  call infolog('nubeam_step.f90: done')
  
900 continue
  call exitlog("nubeam_step_child") 

end subroutine nubeam_step_child
