!-----------------------------------------------------------------------
!  depbox -- use dep to compute fast neutral density in a 3d box aligned
!     on the axis of a select neutral beam
!
!     new dmc june 1992
!     mod dmc july 1992:  support multiple box calculations in one run
!                         if memory permits.
!
subroutine depbox(ierr)
  !
  use nbi_com
  use nbi_mpi_com
  use logmod
  use nbi_halo_mod_def
  use tube_mod, only: deptube_from_box

  !
  !-----------------------------------------------------------------------
  !
  !  nbbox = number of boxes (one for each beam / energy fraction combo)
  !  nbsbox(...) = index to a neutral beams
  !  nbebox(...) = energy fraction indices
  !
  !  nadrbox(...) = output data address for this combo  (see addn0(...))
  !
  !    the neutral density in front of the indicated beam for the
  !    indicated energy fraction will be computed
  !
  !-----------------------------------------------------------------------
  !
  !  1.  zero out the output array
  !
  !============
  implicit none
  integer :: isizt,isiztt,idep0
  integer :: ierr,ii,ib,is,mpierr,ig
  !============
  real(fp) :: zfac
  integer :: nsplt_in,i_nbs_set,i_nbs_set_3d
  character(len=300) :: filename
  character(len=9) :: cmyid,cnumprocs

  if(myidd.eq.0)then
     write(nonlin,*)' in depbox'
  end if
  call enterlog("depbox") ! all returns through exitlog()
  !============
  ierr=0
  !
  isizt=nxbox*nybox*nlbox
  isiztt=isizt*nbbox

  !      resetr call commented out and replaced w/ f90 assign jed 8/01
  !      call resetr(boxn0,isiztt,0.0)


  idep0=mdepbox             !it supports mpi
  mdepbox=-mdepbox          ! this signals dep that calls are from depbox
  ! cleanup
  boxn0=0.0_fp
  xninjs_box=0.0_fp
  xninja_box=0.0_fp
  nsample_box=0
  ! recalculate xninj, xninjs and xninja for beam in box  subset
  do ii=1,nbbox
     ib=nbsbox(ii)
     if(nbsbox(ii).eq.0) cycle
     zfac=0 !get 0 in output if an error in input data
     if(nbebox(ii).eq.1) zfac=ffulla(ib)
     if(nbebox(ii).eq.2) zfac=fhalfa(ib)
     if(nbebox(ii).eq.3) zfac=1.0_fp-ffulla(ib)-fhalfa(ib)
     xninja_box(ii)=xninja(ib)*zfac !beam neutral current in set boxes
     do is=1,nsbeam
        if(lsbeama(ib).ne.is) cycle
        xninjs_box(is)=xninjs_box(is)+xninja_box(ii) !beam neutral current by species for box subset
     end do
  end do
  xninj_box=sum(xninja_box) !#/sec

  if(xninj_box.eq.0.0_fp) then
     mdepbox=idep0          ! restore ndepbox
     if(nlfbmflr) then
        fbm_ptcl=0.0_fp     ! clean it up before the orbiting
     else
        fbm=0.0_fp
     end if
     call track_halo_dalloc ! deallocate itrackh%
     return                 ! there is no particle current, nothing need to be saved
  end if
  !
  i_trk0=1
  nbi_den0mn=0.0_fp
  do while(i_trk0.le.idep0) ! idep0 is how many 0 generation neutral tracks need to be deposited
     ! to a box/s
     !
     call nbidep(0,ierr)    !0 generation only
     i_trk0=i_trk0+1        !0 generation neutral trk deposited to a box not only in plasma region
     if(ierr.ne.0) then
        boxn0=0.0_fp
        mdepbox=idep0       ! restore ndepbox
        if(nlfbmflr) then
           fbm_ptcl=0.0_fp  ! clean it up before the orbiting
        else
           fbm=0.0_fp
        end if
        call track_halo_dalloc ! deallocate itrackh%
        return
     end if
     !
  end do
  !
  boxn0=boxn0*xninj_box/ndepbox_use !weight is taken into acount for a 0 generation neutrals
  !
  if((mdepbox.lt.0).and.(levmod_halo.eq.2))then
     !     weight normalization for a 0 generation neutrals
     do ii=1,idep0
        itrackh(ii)%weight=xninj_box/ndepbox_use
     end do

     if(nsplt.gt.0) then ! goto 3d halo, else skip and count only 0 generation of neutrals
        !mg 10/16/2015           il_3dhalo_den=.true.
        !mg this work has been finished yet, need to compute outputs that frantic
        !mg does nubeam_frantic,
        !mg so i will set  il_3dhalo_den=.false.
        !mg to use frantic
        !mg some day i'll do it.
        il_3dhalo_den=.false.
        !
        ! check if it is eligible  to use halo density, vol av halo e and
        ! vol av halo omega from 3d halo model (not frantic one, set in parms.for)
        ! all beams have to be set in 3d-boxes
        !

        if(nbeam.ne.maxval(nbsbox)) then
           il_3dhalo_den=.false.
        else
           ib=1
           do while (ib.le.nbeam)
              i_nbs_set=0
              if(ffulla(ib).gt.0.0_fp)then
                 i_nbs_set=i_nbs_set+1 !full energy current fraction
              end if
              if(fhalfa(ib).gt.0.0_fp)then
                 i_nbs_set=i_nbs_set+1 !half energy current fraction
              end if
              if(ffulla(ib)+fhalfa(ib).lt.1.0_fp)then
                 i_nbs_set=i_nbs_set+1 !third energy current fraction
              end if
              i_nbs_set_3d=count(nbsbox.eq.ib)
              if(i_nbs_set_3d.ne.i_nbs_set) then
                 il_3dhalo_den=.false.
                 exit
              end if
              ib=ib+1
           end do

        end if
        wmin_halo=xninj_box/(nsplt*ndepbox_use)
        if(il_3dhalo_den) then
           voln0=0.0_fp
           evol0=0.0_fp
           ovol0=0.0_fp
           voln02=0.0_fp
           evol02=0.0_fp
           ovol02=0.0_fp
        end if
        !     now generate a next generation neutral tracks for 3d halo model
        call depbox_all(ierr) !call for 3d halo model for all generation > 0
        if(ierr.ne.0) then
           call errorlog(' emergency exit depbox: depbox_all error')
           if(nlfbmflr) then
              fbm_ptcl=0.0_fp ! clean it up before the orbiting
           else
              fbm=0.0_fp ! clean it up before the orbiting
           end if
           call track_halo_dalloc ! deallocate itrackh%
           return
        end if
     end if
     ! netcdf file name for 3d neutrals data
     write(cmyid,'(i5)') myidd
     cnumprocs=' '
     write(cnumprocs,'(i5)') numprocd
     filename=trim(workpath)//trim(runid)//'_3dhalo_'// &
          trim(adjustl(cmyid))//'_'// &
          trim(adjustl(cnumprocs))//'.cdf'
     call nbi_3dhalo_write(filename)
     if(myidd.eq.0)then
        write(nonlin,*) &
             '   depbox: data for 3d halo are in file '// &
             trim(filename)
     end if
     nsplt_in=abs(nsplt) !mg no needs to fill tubes in
#ifdef __DEBUG
     nsplt_in=-1
#endif


     if(nbtube.gt.0.and.nsplt_in.le.0) then !eflux diagnostic on
        call deptube_from_box(ierr)
        if(ierr.ne.0)then
           call errorlog &
                (' deptube_from_box: error')
           return
        end if
     end if

  end if

  !
  !     reduce data for boxn0
#ifdef __MPI
  call beam_in_box_reduce(il_3dhalo_den,ierr)
  if(myidd.eq.0)then
     write(nonlin,*) &
          '   depbox: data for beam in box have been reduced'
  end if
#endif
  boxn0=boxn0/dvbox
  if(il_3dhalo_den) then
     do ig=1,mig
        do ii=1,mzones
           if(voln0(ig,ii).gt.0.0_fp) then
              evol0(ig,ii)=evol0(ig,ii)/voln0(ig,ii)
              ovol0(ig,ii)=ovol0(ig,ii)/voln0(ig,ii)
           else
              evol0(ig,ii)=0.0_fp
              ovol0(ig,ii)=0.0_fp
           end if
           if(dvolzc(ii).gt.0.0_fp.and.voln0(ig,ii).gt.0.0_fp)then
              nbi_den0mn(ig,ii,:)=nbi_den0mn(ig,ii,:) &
                   /dvolzc(ii)*1.e+06!#/m-3
              voln0(ig,ii)=voln0(ig,ii)/dvolzc(ii) !#/cm-3
           end if
        end do
        do ii=1,mimxbzf
           if(voln02(ig,ii).gt.0.0_fp)then
              evol02(ig,ii)=evol02(ig,ii)/voln02(ig,ii)
              ovol02(ig,ii)=ovol02(ig,ii)/voln02(ig,ii)
           else
              evol02(ig,ii)=0.0_fp
              ovol02(ig,ii)=0.0_fp
           end if
           if(bmvol(ii).gt.0.0_fp.and.voln02(ig,ii).gt.0.0_fp)then
              voln02(ig,ii)=voln02(ig,ii)/bmvol(ii) !#/cm-3
           end if
        end do
     end do
     call fast_halo_sigvn !update atomic physics as new voln02 computed
  end if

  mdepbox=idep0             ! restore ndepbox
  call track_halo_dalloc    ! deallocate itrackh%
  call exitlog("depbox") ! all returns through exitlog()
  if(myidd.eq.0)then
     write(nonlin,*)' exit depbox'
  end if
  return
end subroutine depbox
! 11sep2004 fgtok
