subroutine get_fbm_ncheck0(znum,znum2)
  use iso_c_binding, only: fp => c_double
  use get_fbm_module
  implicit none

  !  compute total #ptcls for species (ns), two ways

  real(fp), intent(out) :: znum,znum2  ! results from fbm; results from bdens2

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

  real(fp) :: densf(nznbmrf),densb(nznbmrf)

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

  call get_fbm_ncheck1(densf,densb,znum,znum2)

end subroutine get_fbm_ncheck0

subroutine get_fbm_ncheck1(densf,densb,znum,znum2)
  use iso_c_binding, only: fp => c_double

  use get_fbm_module
  implicit none

  !  compute 1d density profiles total #ptcls for species (ns), two ways

  real(fp) :: densb(nznbmrf),densf(nznbmrf)  ! from fbm; from bdens2
  real(fp), intent(out) :: znum,znum2  ! results from fbm; results from bdens2

  !--------------------
  real(fp) :: vol(nznbmrf)
  integer :: ix,ith,ie,ia
  real(fp) :: da,de
  real(fp) :: zsum,zdnw
  real(fp) :: fbm_fnc
  integer :: ith1,ith2,nth
  !--------------------

  densb=0.0_fp
  densf=0.0_fp
  vol=0.0_fp
  znum=0.0_fp
  znum2=0.0_fp
! 
  da=1.0_fp/nznbma
!
  do ix=1,nznbmrf
     ith1=nthzsm(lcentr+ix-2)+1
     ith2=nthzsm(lcentr+ix-1)
     nth=ith2-ith1+1
     !
     !  sum over theta.  Note that bdens2 COMMON array is in #, not #/cm3
     !
     zsum=0.0_fp
     do ith=ith1,ith2
        vol(ix)=vol(ix)+bmvol(ith)
! bdens2
        densb(ix)=densb(ix)+bdens2(ith,ns)
        !
        !  sum over energy & pitch
        !
        do ia=1,nznbma
           do ie=1,nznbmea(ns)
              zdnw=fbm_fnc(ie,ia,ith,ns)
              de=efbmb(ie+1,ns)-efbmb(ie,ns)
              zsum=zsum+ &
                   zdnw*dE*dA*bmvol(ith)
           end do
        end do
     end do
     densf(ix)=zsum
!
     znum=znum+densf(ix)
     znum2=znum2+densb(ix)
     densf(ix)=densf(ix)/vol(ix)
     densb(ix)=densb(ix)/vol(ix)
  end do

end subroutine get_fbm_ncheck1
