subroutine r8_boxn0_setnum(miboxm,nbsbox,nbebox, &
       xboxhw,yboxhw,xlbox1,xlbox2, &
       nxbox,nybox,nlbox, &
       nonlin, nbeam, xzbeama, &
       dvbox,nbbox,nadrbox,ierr)  ! output args on this line
 
  use iso_c_binding, only: fp => c_double
  implicit none
 
  !  *** note no modules *** all local variables
  ! because called from outside and inside nubeam

  integer, intent(in) :: miboxm    ! nbsbox nbebox nadrbox array dim.
 
  integer, intent(in) :: nbsbox(miboxm)    ! beam sce #s
  integer, intent(in) :: nbebox(miboxm)    ! energy fraction #s
 
  real(fp), intent(in) :: xboxhw,yboxhw,xlbox1,xlbox2    ! box dimensions
 
  integer, intent(in) :: nxbox,nybox,nlbox ! box grid sizes
  integer, intent(in) :: nonlin     ! file for error messages
  integer, intent(in) :: nbeam      ! number of beams
  real(fp),    intent(in) :: xzbeama(*) ! atomic number of beams
 
  !  output...
 
  real(fp), intent(out) :: dvbox               ! size of one cell in box
  integer, intent(out) :: nbbox            ! no. of boxes
  integer, intent(out) :: nadrbox(miboxm)  ! addresses of boxes
  integer, intent(out) :: ierr             ! status code; 0=ok
 
  !$r8real_input: xboxhw,yboxhw,xlbox1,xlbox2,xzbeama
  !$r8real_output: dvbox
  !----------------------------------
  !  local...
 
  integer ibx,iadr0,isizt,iminbx,iermax
  real(fp) zmindl
 
  !----------------------------------
 
  ierr=0
 
  dvbox=0.0
  nadrbox=0
  nbbox=0
  do ibx=1,miboxm
     if(nbsbox(ibx).eq.0) exit
     nbbox=ibx
  enddo
  if(nbbox.eq.0) return

  iadr0=0
  iermax=0
  do ibx=1,nbbox
     call r8_boxn0_errck(ibx,nbsbox(ibx),nbebox(ibx),nonlin,nbeam,xzbeama,iermax)
  enddo
  if(iermax.gt.0) then
     nbbox=0
  endif
 
  iminbx=min(nxbox,nybox,nlbox)
  if(iminbx.le.0) then
     call boxn0_errmsg('?r8_boxn0_setnum: min(nxbox,nybox,nlbox).le.0',nonlin)
     nbbox=0
  else
     isizt=nxbox*nybox*nlbox
  endif
 
  zmindl=min(xboxhw,yboxhw,(xlbox2-xlbox1))
  if(zmindl.le.0.0) then
     call boxn0_errmsg( &
          '?r8_boxn0_setnum: min(xboxhw,yboxhw,(xlbox2-xlbox1)).le.0',nonlin)
     nbbox=0
  endif
 
  if(nbbox.eq.0) then
     ierr=1
     call boxn0_errmsg( &
          '?r8_boxn0_setnum: beam in box neutral density calculation disabled.',nonlin)
     return
  endif
 
  dvbox=(2.0_fp*xboxhw)*(2.0_fp*yboxhw)*(xlbox2-xlbox1)/(isizt)
 
  iadr0=0
  isizt=nxbox*nybox*nlbox
  do ibx=1,nbbox
     nadrbox(ibx)=iadr0+1
     iadr0=iadr0+isizt
  enddo
 
end subroutine r8_boxn0_setnum
 
subroutine boxn0_errmsg(zstr, nonlin)
 
  implicit none
 
  character*(*) zstr
  integer       nonlin
 
  write(nonlin,*) trim(zstr)
 
  return
 
end subroutine boxn0_errmsg
 
subroutine r8_boxn0_errck(ibx,ibs,ibe,nonlin,nbeam,xzbeama,iermax)
 
  use iso_c_binding, only: fp => c_double
  implicit none
 
  integer, intent(in) :: ibx        ! box no.
  integer, intent(in) :: ibs        ! beam sce no.
  integer, intent(in) :: ibe        ! energy fraction no.
 
  integer, intent(inout) :: iermax  ! error count (may be incremented)
  integer, intent(in) :: nonlin     ! file for error messages
  integer, intent(in) :: nbeam      ! number of beams
  real(fp),    intent(in) :: xzbeama(*) ! atomic number of beams
  ! -----------------------
  integer imaxe
  ! -----------------------
 
  if((ibs.le.0).or.(ibs.gt.nbeam)) then
     write(nonlin,*) ' ?r8_boxn0_errck:  error in 3d beam in box calculation:'
     write(nonlin,*) '  box no. ',ibx,' refers to non-existant source no. ',ibs
     iermax=iermax+1
  else
     if(abs(xzbeama(ibs)-1.0_fp).le.0.05_fp) then
        imaxe=3     ! 3 energy fractions for h isotopes
     else
        imaxe=1     ! just 1 for all others
     endif
     if((ibe.le.0).or.(ibe.gt.imaxe)) then
        write(nonlin,*) ' ?r8_boxn0_errck:  error in 3d beam in box calculation:'
        write(nonlin,*) '  box no. ',ibx, &
             'refers to non-existant energy fraction: ',ibe
        write(nonlin,*) '  sce no. = ',ibs,' zbeam(sce)=',xzbeama(ibs)
        iermax=iermax+1
     endif
  endif
 
  return
end subroutine r8_boxn0_errck
