!----------------------------------------------------------------------
!  nbstart -- set up for another step on the beam calculation
!
subroutine nbstart(ierr)

  use nbi_com
  use nbi_idhdr
  use nbi_random
  use nbi_mpi_com
  use mpi_env_mod, only: mpi_env_comm,  loc_myid, num_procs
  use logmod

  implicit none

  integer, intent(out) :: ierr    ! completion code (0=ok)
  integer :: insbeam
  real(fp) :: time1,time2
  integer :: mpierr
  !
  !-------------------
  !
  !  dmc jan 1992 -- changes to support heavy neutral beams
  !
  !
  !  dmc july 1990 -- this code was formerly duplicated in routines
  !   nbmain and fpmain -- now consolidated here
  !
  !  this routine calls further routines to zero out monte carlo sums,
  !  and compute tables for use in the beam deposition and slowing down
  !  calculation, based on current plasma parameters
  !
  !--------------------
  !%      external philb,rygrad
  !--------------------
  !
  real(fp) :: ztest,ztim1,ztim2
  integer :: isb,ib,iwarn,ien
  integer :: ilen
  real(fp) :: ztime1,ztime2,ztime3,ztime4,ztime5,ztime6
  real(fp) :: ztime7,ztime8,ztime9,ztime10,ztime11
  real(fp) :: zdtim11,zdtim12,zdtim21,zdtim22,zdtim31,zdtim32
  !
  integer :: iuns(2)
  !
  integer :: initpl
  data initpl/0/
  integer :: reckon_minb      ! integer function to calculate minb
  character(len=5) :: cib,cie(3)
  !cpu time elapsed in nbatom
  real(fp) :: ctime_nbatom_in,ctime_nbatom_out,dctime_nbatom, &
       dctime_nbatom_all
  !wall time elapsed in nbatom
  integer :: itime_nbatom_in,itime_nbatom_out,idtime_nbatom, &
       idtime_nbatom_all


#ifndef __MPI
  real(fp) :: zwrk(mj)
  integer :: myid_debug, numproc_debug
  !     this will be a path for nubeam_test, when debugging mpi nubeam for one ptcl.sublist
  if(debug_mpi.and.nltest_output) then
     myid_debug=myidd
     numproc_debug=numprocd
  end if
#endif

  call nbi_mpi_com_init

  call infolog('nbstart: nbi_comm = ',nbi_comm)

#ifndef __MPI
  if(debug_mpi.and.nltest_output) then
     myidd=myid_debug       ! reset, do not redefine  myidd and nuprocd for
     numprocd=numproc_debug ! debugging of mpi version by serial code
  end if
#endif
  !--------------------
  !
  ierr=0
  !
  ilen=len(trim(workpath))
  if(ilen.eq.0) then
     workpath='./'
  else
     if(workpath(ilen:ilen).ne.'/') then
        ilen=ilen+1
        workpath(ilen:ilen)='/'
     end if
  end if
  !
  insbeam=nsbeam
#ifdef __MPI
  ttime11=mpi_wtime()
#endif

  if(only_io) then
     call errorlog('nbstart: improper call, only_io = .true.')
     ierr=666
     return
  end if
  !
  call cpu_time(ztim1)
  !
  !  find available i/o unit numbers
  !
  call find_io_unit_list(2,iuns)
  if((lunnbx.le.0).or.(lunnbx.gt.99)) lunnbx = iuns(1)
  if((lunres.le.0).or.(lunres.gt.99)) lunres = iuns(2)
  !
  !-------------------------------------------------------------
  !  check 2d mc grid & call bmzini...
  !
  call nbstart_init(ierr)
  if(ierr.ne.0) return

  !  *** particle list size minb ***
  !

  if(nbsteps.eq.0) then
     minb = reckon_minb(numprocd) !find # of particles (if mpi -  per each processor)
     call nbi_alloc_minb    !allocate particle list arrays
  end if

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


#ifdef __MPI
  ztime1=mpi_wtime()
#endif
  !
  ilunrng=nonlin


#ifdef __MPI
  ztime2=mpi_wtime()
#endif
  !


  if(dt_adj0.le.0.0_fp) dt_adj0=1.0_fp
  do isb=1,nsbeam
     if(dt_adj0s(isb).le.0.0_fp) dt_adj0s(isb)=1.0_fp
  end do
  !
  ! save state -- for standalone code, in case of hard crash
  ! (for debug restart in standaline beam code we will need:
  !     --xplasma files
  !     --state of nubeam calculation at beginning of step
  !
  ! the latter is taken care of here; but first, translate xplasma id
  !  integr codes into name strings
  !
  if((.not.nbi_idhdr_skip).and.(nltest_output)) then

     if (myidd.eq.0) then

        nbi_curhdr=0
        call nbi_id_save    ! save xplasma names
        !
        ! save xplasma state
        call eq_save(workpath(1:ilen)// &
             runid(1:lrunid)//'_debug_xplasma_nubeam.cdf',ierr)
        if(ierr.ne.0) then
           call errorlog('nbstart:failed to write:'//runid(1:lrunid) &
                //'_xplasma_nubeam.cdf')
           call errorlog('nbstart: error ignored')

           ierr=0
        end if
        !
     end if
     ! save nubeam state

     call nbi_states(workpath(1:ilen)// &
          runid(1:lrunid)//'_debug_nbi_fld_state.cdf', &
          '_debug_nbi_ptcl_state_cpu', &
          's',0,ierr) ! see nubeam_test

     if(ierr.ne.0) then

        if(myidd.eq.0) then
           call errorlog( &
                'nbstart:  failed to write:  '//runid(1:lrunid)// &
                '_nbi_fld(ptcl)_state.cdf')

           call errorlog('nbstart: error ignored')

           ierr=0
        else
           call errorlog( &
                'nbstart: failed to write:debug_nbi_ptcl_state_cpu')

           call errorlog('nbstart: error ignored')

           ierr=0
        end if
     end if

  end if

#ifdef __MPI
  ztime3=mpi_wtime()
#endif
  !
  !----------------------------------
  ! warnings: unusual switch settings...
  iwarn=0

  !      if(nlbout) then
  !         call warn(
  !     >        'nlbout=t, <runid>nb.out file output: disabled!!')
  !         nlbout=.false.
  !      end if

  if(nlfdep) then
     call warn( &
          'nlfdep=t, deposition file output: disabled!!')
     nlfdep=.false.
  end if

  if(nlbfpp) then
     call warn( &
          'nlbfpp=t, deposition data output using nltrk_dep0=t.')
     !         nltrk_dep0=.true.
  end if

  if(goocon.lt.(9.9_fp)) call warn('goocon < 10')
  if(nchdvp.ne.3) call warn('nchdvp .ne. 3')
  call t1warn(xdepmod,'xdepmod .ne. 1.0; deposition modified!')
  if(.not.nlbeamcx) then
     nlbbcx=.false.
     call warn('nlbeamcx=f, beam cx losses disabled!')
  end if
  if(.not.nlbbcx) then
     call warn('nlbeamcx.or.nlbbcx = f, no beam-beam deposition.')
  end if

  call t1warn(xcfanbi,'xcfanbi .ne. 1.0; electron drag changed!')
  call t1warn(xdfanbi,'xdfinbi .ne. 1.0; ion drag & edif changed!')
  call t1warn(xefanbi,'xefanbi .ne. 1.0; pa scattering changed!')

  call t1warn(xcfafus,'xcfafus .ne. 1.0; electron drag changed!')
  call t1warn(xdfafus,'xdfifus .ne. 1.0; ion drag & edif changed!')
  call t1warn(xefafus,'xefafus .ne. 1.0; pa scattering changed!')

  if(.not.nlbcde) call warn('nlbcde=f, energy diffusion disabled!')
  if(.not.nlbcoh) call warn('nlbcoh=f, vloop accel disabled!')
  if(.not.nlbcpa) call warn('nlbcpa=f, pa scattering disabled!')
  if(nlorbo) call warn('nlorbo=t, orbit only mode in effect!')

  if(.not.nlbflr) call warn('nlbflr=f, flr corrections turned off!')
  if(nrip.ne.0) call warn('nrip .gt. 0, ripple loss model active!')
  if(nlfbon) call warn('nlfbon=t, fishbone model is active!')
  if(nbbox.gt.0) call warn('nbbox > 0, beam in box calculation.')

  if(.not.nlcprb) call warn('nlcprb=f, db/dt correction disabled!')
  if(nmcurb.eq.0) call warn('nmcurb=0, beam current drive off!')

  if(iwarn.gt.0) call warnlog(' ===============================')

  !----------------------------------
  !  some error checks -- to stop run immediately

  if(goocon.lt.5.0_fp) then
     call errorlog('nbstart: detected goocon = ',goocon)

     call errorlog( &
          'nbstart: a minimum value of 5.0 is required for this input.')

     ierr=5
     return
  end if
  !
  !----------------------------------
  !
  ngoocon_vpvbin=max(0,ngoocon_vpvbin)
  ngoocon_ebin=max(1,ngoocon_ebin)
  ngoocon_rbin=max(1,ngoocon_rbin)

  if((goocon.gt.(10.001_fp)).or.(ngoocon_vpvbin.gt.0).or. &
       (ngoocon_ebin.gt.1).or.(ngoocon_rbin.gt.1)) then
     if(dt_acc.gt.1.0e-4_fp) then
        call warnlog( &
             'nbstart: non-default goocon controls detected.')

        call warnlog( &
             'nbstart:  therefore dt_acc <= 1.0e-4 is enforced.')
        dt_acc = 1.0e-4_fp

     end if
  end if
  !
  if(nmdifb.eq.4) then
     if(dt_acc.gt.1.0e-4_fp) then
        call warnlog( &
             'nbstart: pitch dependent radial trnasport db(e,x,t)')
        call warnlog( &
             'nbstart:  therefore dt_acc <= 1.0e-4 is enforced.')
        dt_acc = 1.0e-4_fp
     end if
  end if
  !
  !----------------------------------
  !  set some flags now
  !
  !   if initial goose=1 bounce is needed when resuming orbiting of a partially
  !   slowed down ion, set the switch here

  nlgoose_reset = .false.
  if(orbrzv_rf_option.eq.1) nlgoose_reset = .true. ! rf kick
  if(nlfbon) nlgoose_reset = .true. ! fishbone model needs this
  if(nmdifb.eq.4) then
     nlgoose_reset = .true. ! db(e,x,t) model also needs this
     ngoocon_ebin = max(1,ngoocon_ebin)
     ngoocon_rbin = max(1,ngoocon_rbin)
     if(dif2d_dpa.le.0.0_fp) dif2d_dpa = 1.0_fp/goocon
     dif2d_dpa=max(1.0e-3_fp,dif2d_dpa)  ! nmdifb.eq.4 pitch variance ctrl
  end if

  fbtrap_depth = max(0.0_fp,min(1.0_fp,fbtrap_depth))
  nfbon_vpvopt = max(0,min(1,nfbon_vpvopt))
  nfbon_species= max(1,min(3,nfbon_species))

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

  call nbi_nbsums                   !  update max beam voltage & power sums
  ! (above call sets xninj, xninja(1:nbeam)
  call nbi_check_beam_species ! check if beam is on to get track data (from bmzini.for)
  ! nltrk_dep0 could be changed there
  ! check for deposition neutral track data gathering...

  ! standard settings for track data gathering:

  ndep0_use = max(ndep0_use,ndep0)     ! gather this many...
  ! (this controls deposition statistics, even if no track data gathering
  ! is done)

  frac_cur_ndep0 = 0.0  ! no data gathering if left at 0.0
  ndep_beam_use = 0     ! 0 means, all beams
  ndep_ien_use  = 0     ! 0 means, all energy fractions

  if (nltrk_dep0) then
     !  track data gathering requested...

     if(ndep_set_beam.gt.0) then
        call check_beam_current !check if chosen beam and energy fraction
        !has enough ptcls. current (more than 1%)
        !if not -track data will be collected for
        !all beam sources, if the  ptcls. current
        !more than 1% but less 10% then the number
        !collected is reduced.
     else
        frac_cur_ndep0 = 1.0    !collecting tracks on all beams
     end if
  end if


  !----------------------------------
  !  define lsbeama
  do ib=1,nbeam
     do isb=1,nsbeam
        if(.not.nlfprod(isb)) then
           if((abs(abeama(ib)-abeams(isb)).le.0.05_fp).and. &
                (abs(xzbeama(ib)-xzbeams(isb)).le.0.05_fp)) then
              lsbeama(ib)=isb
              exit
           end if
        end if
     end do
  end do

  !----------------------------------
  ! set aep opertaor if available
  nphipzc=0
  if(maxval(abs(aep_ascale)).ne.0) then !put aep data to black box to be passed to nubeam
! if aep_ascale=-1.0 then RBQ will be applied
     call nbi_aep_get(nkaep,nsbeam,nonlin,ierr)  !aep data from blkbx
   end if


  !
  !----------------------------------
  !  clear output grid ids
  !
  id_rhozb=0
  id_rhozc=0
  !
  !----------------------------------
  !  miscellaneous setup...
  !
  if((nclass.gt.0).and.(fppcon.lt.16.0_fp)) then
     fppcon=16.0_fp
  end if
  !
  !  sawtooth operator...
  !
  call pre_sawnbi_setup

#ifdef __MPI
  ztime4=mpi_wtime()
#endif
  !
  !----------------------------------
  ! c% lines for debug plotting
  !
  !
  !%      nbplist=1
  !%      if(initpl.eq.0) then
  !%        call uinita(0)
  !%        initpl=1
  !%      end if
  !
  !  timestep
  if (delt.ne.0.0_fp) then
     if(tbm1.ne.tbm2p) then
        call warnlog( 'nbstart: warning -- ')
        call warnlog('   stop time of previous timestep = ',tbm2p)
        call warnlog('   start time of current timestep = ',tbm1)

     end if
     tfracp = (tbm2-tbm1)/delt
     bdenssp = bdenss
  else
     ! must be the first timestep (saved delt from prev. step is zero).
     tfracp = 1.0_fp
     bdenssp = 0.0_fp
     bdenss  = 0.0_fp
     frac_survive = 1.0_fp
  end if

  delt=(tbm2-tbm1)
  tbm2p=tbm2

  if(delt.le.0.0_fp) then
     call errorlog( &
          'nbstart: a finite positive time step is required.')
     call errorlog('  tbm1: ',tbm1)
     call errorlog('  tbm2: ',tbm2)
     call errorlog('  difference: ',tbm2-tbm1)
     call errorlog('  it is ok if you run nubeam_test')

     ierr=111
     !mg i commented as it is ok to run nubeam_test         return
  end if
  !
  !  there used to be a loop here to adjust goosay elements based on the
  !  tfracp value -- this has been subsumed into the parallelized orball
  !  particle loop.  dmc nov. 2005
  !
  if (tfluxp.eq.0.0_fp) then
     eldot=0.0_fp
  else
     eldot=(tflux-tfluxp)/(tflux+tfluxp)
  end if
  !
  !  init output smoothing routines
  !
  call r8bsmoo_init(mj,lcentr,ledge,abs(dxbsmoo),xi_r8)

  !  initialize gaussian sampler (no longer necessary)
  !xxx      call sgauss
  !
  !  re-evaluate beam zone volumes (actually, fetch xplasma values)
  !
  call bmzvol(ierr)
  if(ierr.ne.0) return
  !
  !  open special output file(s) (if any)
  !
  if(nlbout) then
     mlost_data = 0 ! zeroed number of lost particles
  end if
  !

  nbs=0                  ! clear species setting
  lsbeam=0
  abeam=0.0_fp
  xzbeami=0.0_fp
  !
  !---------------/end dmc 1-dec-1998 change
  !
  !  set 2d thermal neutral density, etc.
  !
  call parms(ierr)
  if(ierr.ne.0) return
  !
  ! extrapolate beam-beam relative energy average if needed
  call bb_extrp(ierr)
  if(ierr.ne.0) return
  !
  !  compute fusion product sources if needed
  !  [dmc bugfix -- call rdbbsg before fusion: 15 mar 2011]
  !-------------------------------------------
  !  init beam-beam & beam-target neutron calculations
  !
  call rdbbsg
  !
  lsbeamr = 0
  call btfusn_init   ! sets lsbeamr(:) elements
  !
  if(nlfusn) call fusion
  !
  !-------------------------------------------
  !
  !  zero all 1d output arrays
  call nbzero
  !
  !  copy nth generation beam-beam info to past timestep arrays
  call nbcx00(1)
  !
  !  clear 2d profile output arrays
  call nb2d00
  !
  !                              call expert(l1fimain,l2fimain,3)
  !                               calls to expert no longer needed
  !                               jed 7/13/01
  !
  !
  !  dmc jan 1992 -- evaluate cross sections and log(lambda) data for
  !  heavy neutral beams and associated ion species
  !
  call hvsigma(ierr)
  if(ierr.ne.0) return

#ifdef __MPI
  ztime5=mpi_wtime()
#endif
  !
  !  evaluate fokker-planck coeffs
  !
  call nbi_fokker(ierr)
  if(ierr.ne.0) return
  !                              call expert(l1fimain,l2fimain,4)
  !                               calls to expert no longer needed
  !                               jed 7/13/01
  !
  !  evaluate atomic physics for stopping 1st gen and cx fast neutrals
  !  (dep and nutrav)
  !
  !#ifdef __DEBUG
  !      if(nbsteps.eq.2) then
  !         cie(1)='full'
  !         cie(2)='half'
  !         cie(3)='third'
  !      do ib=1,nsbeam
  !            call r8_grafx1(xizc(lcentr:ledge),
  !     >           bnb_ienx(lcentr:ledge,1,ib)*1.e-6,
  !     >           nzones,' ','cm-3',' ',
  !     >           'co_inj','n_fast vs x')
  !            call r8_grafx1(xizc(lcentr:ledge),
  !     >           bnb_ienx(lcentr:ledge,2,ib)*1.e-6,
  !     >           nzones,' ','cm-3',' ',
  !     >           'cnt_inj','n_fast vs x')
  !      end do
  !
  !
  !      do ib=1,nbeam
  !            call infolog('nbstart: beam# = ',ib)
  !            cib=' '
  !            write(cib,'(i5)') ib
  !         do ien=1,3
  !            call r8_grafx2(xizc(lcentr:ledge),
  !     >           ebb_relx_co(lcentr:ledge,ien,ib,lsbeama(ib))*1.e-3,
  !     >           rms_ebbx_co(lcentr:ledge,ien,ib,lsbeama(ib))*1.e-3,
  !     >           nzones,' ','kev',' ',
  !     >           'co_inj beam #'//trim(adjustl(cib))//', energy ='
  !     >           //trim(adjustl(cie(ien))),'e_rel & t_fast vs x')
  !         end do
  !         do ien=1,3
  !            call r8_grafx2(xizc(lcentr:ledge),
  !     >           ebb_relx_cnt(lcentr:ledge,ien,ib,lsbeama(ib))*1.e-3,
  !     >           rms_ebbx_cnt(lcentr:ledge,ien,ib,lsbeama(ib))*1.e-3,
  !     >           nzones,' ','kev',' ',
  !     >           'cnt e_rel#'//trim(adjustl(cib))//', energy ='
  !     >           //trim(adjustl(cie(ien))),'e_rel & t_fast vs x')
  !         end do
  !         do ien=1,3
  !            call r8_grafx2(xizc(lcentr:ledge),
  !     >           ebb_relx_fn(lcentr:ledge,ien,ib,2),
  !     >           rms_ebbx_fn(lcentr:ledge,ien,ib,2),
  !     >           nzones,' ','ev',' ',
  !     >           'fn2 e_rel & rms_e_rel','nbstart')
  !         end do
  !         do ien=1,3
  !            call r8_grafx2(xizc(lcentr:ledge),
  !     >           ebb_relx_fn(lcentr:ledge,ien,ib,3),
  !     >           rms_ebbx_fn(lcentr:ledge,ien,ib,3),
  !     >           nzones,' ','ev',' ',
  !     >           'fn3 e_rel & rms_e_rel','nbstart')
  !         end do
  !      end do
  !      end if
  !#endif
#ifdef __MPI
  mpi_env_comm=nbi_comm !pass mpi_comm_world to mpi_env_mod for adas310_fortran_driver
  loc_myid=myid
  num_procs=numprocs
#endif
  call cpu_time(ctime_nbatom_in)
  call wall_seconds(itime_nbatom_in)

  call nbatom(ierr)

  call cpu_time(ctime_nbatom_out)
  call wall_seconds(itime_nbatom_out)
  dctime_nbatom=ctime_nbatom_out-ctime_nbatom_in
  idtime_nbatom=itime_nbatom_out-itime_nbatom_in

#ifdef __MPI
  !wall time in nbatom
  idtime_nbatom_all=0
  call  mpi_reduce(idtime_nbatom,idtime_nbatom_all,1,MPI_INTEGER, &
       MPI_SUM,0,nbi_comm,mpierr)
  if(myidd.eq.0) then
     idtime_nbatom=idtime_nbatom_all
  end if
  !cpu time in nbatom
  dctime_nbatom_all=0.0_fp
  call  mpi_reduce(dctime_nbatom,dctime_nbatom_all,1,MPI_REAL8, &
       MPI_SUM,0,nbi_comm,mpierr)
  if(myidd.eq.0) then
     dctime_nbatom=dctime_nbatom_all
  end if
#endif
  if(myidd.eq.0) then

     call infolog(' nbstart: cpu time elapsed in nbatom= ', &
          dctime_nbatom)
     call infolog(' nbstart: wall time elapsed in nbatom=', &
          idtime_nbatom)

  end if


  call nbzero_ebb ! zeroed relative energy average
  if(ierr.ne.0) return
  !
#ifdef __MPI
  ztime6=mpi_wtime()
#endif
  !
  !-------------------------------------------
  !
  !  fast-ion orbit time step     (jfs  6 dec 94)
  if ( ldtnfr ) then
     dtn = 1.72e-9_fp* raxis
     ldtnfr = .false.       ! dmc -- restartability bugfix, do only once!
  end if
  !
  call cpu_time(ztim2)
#ifdef __MPI
#else
  cpbstart = cpbstart + (ztim2-ztim1)
#endif

  !     *********************************************
  !     ********************************************
  !      all parallel processing starts here
  !     *********************************************
  !     *********************************************

  stime1=0.0_fp
  stime2=0.0_fp

  !debug      call wchk_debug

  !     itest=2

  !     -----------------------------------
  !      read local state files
  !      and perform calculations
  !     ------------------------------------
  !     if (itest.ge.4) then


#ifdef __MPI
  !      call mpi_barrier(nbi_comm,mpierr)
  stime1=mpi_wtime()
#endif

  ztime7=0.0_fp
  ztime8=0.0_fp
  ztime9=0.0_fp
  ztime10=0.0_fp
  ztime11=0.0_fp


#ifdef __MPI
  stime2=mpi_wtime()
  !      call mpi_barrier(nbi_comm,mpierr)
#endif

#ifdef __MPI
  ttime12=mpi_wtime()

  if (myidd.eq.0) then
     cpbstart = cpbstart + (ttime12-ttime11)
  end if
  !     call infolog('nbstart...')
  !     call infolog('total time in nbstart= ',ttime12-ttime11)
  !     call infolog(' ')
  !     call infolog ('time through nbstart_mcgrid= ',ztime1-ttime11)
  !     call infolog('time for bmzini= ',ztime2-ztime1)
  !     call infolog('time saving files= ',ztime3-ztime2)
  !     call infolog('time through pre_sawnbi_setup= ',ztime4-ztime3)
  !     call infolog('time through hvsigma= ',ztime5-ztime4)
  !     call infolog('time through rdbbsg= ',ztime6-ztime5)
  !     call infolog('time through stime1= ',stime1-ztime6)
  !     call infolog ('total time for file read= ',stime2-stime1)
  !     call infolog ( '(')
  !     call infolog ('total nbi_ptcl_read= ',ztime7)
  !     call infolog('total duplicate arrays= ',ztime8)
  !     call infolog ( 'total ptcl_local_copy_out= ',ztime9)
  !     call infolog(')')
  !     call infolog('time through end of nbstart= ',ttime12-stime2)
  !     call infolog(' ')

#endif

  return

  !-----------------------------
contains

  subroutine t1warn(zval,zstr)
    real(fp) :: zval
    character(len=*) zstr            ! warning string

    ! generate warning if item not close to one...

    ztest=zval-1
    if(abs(ztest).gt.1.0e-3_fp) call warn(zstr)

  end subroutine t1warn

  subroutine warn(zstr)
    character(len=*) zstr            ! warning string

    iwarn=iwarn+1
    if(iwarn.eq.1) then
       call warnlog('nbstart: ===============================')
    end if

    call warnlog(' => caution: '//trim(zstr))

  end subroutine warn

  !-----------------------------
end subroutine nbstart
!
!-----------------------------------------------------------------------
subroutine nbstart_nznbmri
  ! compute quantities related to 2d irregular spatial grid for mc outputs
  ! and distribution function

  use nbi_com
  implicit none
  integer ::itest,i_xbmbnd

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

  if(nlsym2b) then
     nthzsm0=nth0
  else
     nthzsm0=2*nth0
  end if

  if(nzone_fb.gt.(nzones/2)) then
     nzone_fb=nzones
  else
     nzone_fb=max(8,nzone_fb) - 1
     do
        nzone_fb = nzone_fb + 1
        itest = nzones/nzone_fb
        if((itest*nzone_fb).eq.nzones) exit
        if(nzone_fb.gt.(nzones/2)) then
           nzone_fb=nzones
           exit
        end if
     end do
  end if

  nznbmri = nzone_fb
  i_xbmbnd = (xbmbnd -1.0_fp)*nzone_fb+0.000001_fp
  nznbmr = nznbmri + max(3,i_xbmbnd)

end subroutine nbstart_nznbmri

!-----------------------------------------------------------------------
!  this routine can be used to re-establish an xplasma 2d "mcgrid".
!  it may be appropriate to have nbstart call this on a restart test--
!  the driver can force this by setting id_mcgrid=-1 before calling
!  nbstart, using nbi_get_grid and nbi_update_grid calls.
!
!  this should *not* be called from a normal time dependent driver,
!  because in general the quantities nlsym2b and nznbmri will not have
!  been initialized.
!
!  xplasma needs to have been initialized before this is called.
!
subroutine nbstart_mcgrid(id)

  use nbi_com
  use logmod

  implicit none
  integer, intent(out) :: id

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

  integer, parameter :: inphi = 1
  integer :: iudsym,inth0

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

  call nbstart_udsym(iudsym)

  inth0=nthzsm0/iudsym

  if(inth0.eq.0) then
     call errorlog( &
          ' nbstart_mcgrid:  inth0 = 0, cannot initialize 2d grid.')

     id=0
     id_mcgrid=0
  else
     call nubeam_mcgrid_define(inphi,iudsym,inth0,nznbmri,nznbmr,id)
     id_mcgrid=id
  end if

  return
end subroutine nbstart_mcgrid

!-----------------------------------------
subroutine nbstart_udsym(iudsym)
  !  return iudsym=1 if nubeam treats equilibrium as up-down symmetric;
  !  return iudsym=2 if it is treated up-down asymmetric

  use nbi_com

  implicit none

  integer, intent(out) :: iudsym

  if(nlsym2b) then
     iudsym=1
  else
     iudsym=2
  end if

  return
end subroutine nbstart_udsym

!-----------------------------------------
subroutine wchk_debug
  !
  use nbi_com
  use nbi_mpi_com
  use logmod
  implicit none
  !
  !  debug code
  !
  integer :: i,isb,ind
  real(fp) :: zwsum
  integer :: find_ptcl_ind
  !----------------------------------
  do isb=1,nsbeam
     zwsum=zero
     do i=1,nuse(isb)
        ind=nbndex(i,isb)
        ind=find_ptcl_ind(numprocd,myidd,ind)
        zwsum=zwsum + wghtay(ind,isb)
     end do
     call infolog('wchk_debug:')
     call infolog('bntots(isb) = ',bntots(isb))
     call infolog('isb = ',isb)
     call infolog('zwsum = ',zwsum)

  end do
  !
  return
end subroutine wchk_debug
!
!-------------------------------------------------------
! 2d grid initialization; bmzini call can result in additional
! nubeam array allocations
!
subroutine nbstart_init(ierr)

  use nbi_com
  use nbi_idhdr
  use logmod
  !-----------------------------
  implicit none
  !
  integer, intent(out) :: ierr
  !
  integer :: id
  !
  !-----------------------------
  !  make sure a 2d mc grid is defined...
  !
  if(id_mcgrid.le.0) then

     if(nznbmri.eq.0) call nbstart_nznbmri

     if(id_mcgrid.eq.-1) nbi_idhdr_skip=.true.  ! (nubeam_test call)

     call nbstart_mcgrid(id)
  end if
  !
  call bmzini(ierr)
  if(ierr.ne.0) then
     call errorlog('nbstart_init: error in bmzini')

  end if

  return
end subroutine nbstart_init
