!----------------------------------------------------------------------
!  outptb2  -- finishing touches on output of m.c. beam model
!
subroutine outptb2(ierr)
  !
  use xstraln_calls
  use nbi_com
  use physconst_mod, only : evptocm_sec
  use nbi_mpi_com
  use xplasma_obj_instance
  use netcdf
  !
  ! this is only set in case of a "severe" error:
  implicit none
  integer, intent(out) :: ierr   ! status code (0=ok)
  !-------------
  !
  !    dmc july 1998 -- to support higher resolution radial grid in rest
  !    of code, smoothing of radial profiles is now supported in the mc
  !    code.  ** profiles renormalized in outptb2 should also be smoothed
  !    in subroutine orbsmoo2 **
  !    see also:  outptb, orbsmoo.
  !    and see also:  depnrm, depsmoo.
  !
  real(fp) :: zjac(2,2)
  integer :: nk
  integer :: jl(nsbeam)
  character(len=150) :: fname       !netcdf file name for fbm_ptcl
  integer :: icdf
  integer :: i_sb
  integer :: id_mimxbe,id_mimxba,id_mimxbzf,id_mgyro
  integer, dimension(:), allocatable :: id_fbm_ptcl, id_e

  !

  !-------------
  !
  !  in the process of generalizing the montecarlo beam code for multiple
  !  beam ion species, this new routine was created.  much of the code
  !  comes from the old routine outptb
  !
  !  here are the renormalization subroutines in the m.c. beam code:
  !
  !   >depnrm renormalizes monte carlo sums computed while depositing
  !    new fast ions; depnrm is called from inside a fast ion species
  !    loop after the deposition loop for said species is complete.
  !
  !   >outptb renormalizes monte carlo sums computed while following
  !    the orbits of a single ion species (i.e. outptb is inside the
  !    fast ion species loop, called after the orbit loop for that
  !    species)
  !
  !   >outptb2, operating outside the fast ion species loop (cf subroutine
  !    orball), computes sums over species, and, renormalizes quantities
  !    that are summed inside the orbit loop with an explicit species
  !    index.  also, final auxilliary calculations, such as cross-species
  !    beam-beam fusion, are called out of outptb2.
  !
  !   outputs normalized by depnrm and outptb are summed or otherwise
  !   used in outptb2.
  !
  !--------------------------------------------
  !
  !  table of neutron producting fusion reactions
  !  see the comments in routine fusion.fnb for description of
  !  reaction indexing
  !
  integer,parameter :: itreact=16
  integer :: ineuts(itreact)
  !
  !  0's for computed fusion reactions that do not produce neutrons
  !  also a 0 for the fusion-product burnup neutrons category (#9)
  !  which is counted separately cf sumntx.for (dmc 17 oct 1991)
  !
  !--------------------------------------------
  character(len=32) :: zlbls(mibs),zlbl  ! fast ion labels
  character(len=32) :: zlbls_nb          ! beam line label
  character(len=80) :: label
  character(len=32) :: zlblobj           ! mc 2d gridded output label
  character(len=10) :: zlbth(mig)
  character(len=2) :: zrnum
  integer :: ilbls(mibs),ilbth(mig),ilbl
  !
  real(fp) :: bsumint,zwnorm
  real(fp) :: zdti,zbpf,ztqbi
  real(fp) :: zv2av,zvzav,zf,zn0
  real(fp) :: zeav,zeav2,zde2,znt
  integer :: isb,isb2,iz,j,jp1,i,inum,ia,jerr,idpro
  integer :: ig,ien,ib,ireact
  integer :: ie,iplot,iadr,iadw,iadrf,iadrw
  integer :: iinc,inx,ix,ityp
  real(fp) :: zconv,zconve,zwsum,zws1
  real(fp) :: zsumjb,zsumjbco
  !
  real(fp) :: zbzxr,zr2i,zri,zb2,get_j_from_jdotb,xjbfacnb
  !
  integer :: id_vpv
  real(fp) :: zxpa(nznbma+1)
  character(len=32) :: nbnc_lbl
  real(fp) :: evptocm_sec2
  !
  !----------------------------------------------------------------------
  !
  data ineuts/1,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0/
  !
  !----------------------------------------------------------------------
  !  fast ion name strings...
  !
  ierr=0
  !
  zlbls=' '
  call nbfi_lbl(zlbls,mibs,inum)
  if(inum.ne.nsfast) then
     write(nonlin,*) ' ?outptb2: expected: inum = nsfast!'
     write(nonlin,*) '  inum = ',inum,' nsfast = ',nsfast
     write(nonlin,*) '  inconsistency in number of fast species!'
     ierr=1
     return
  end if
  do i=1,inum
     ilbls(i)=len_trim(zlbls(i))
  end do
  !
  do ig=1,ng
     ia=aplasm(ig)+0.5_fp
     iz=backz(ig)+0.5_fp
     call namels(iz,ia,zlbth(ig))
     ilbth(ig)=len_trim(zlbth(ig))
  end do
  evptocm_sec2= evptocm_sec*evptocm_sec
  !
  !----------------------------------------------------------------------
  !
  !  timestep normalization
  !
  zdti=1.0_fp/delt
  zbpf=xe*zdti
  !
  !-----------------------------------------------------------------
  !  sum scalar outputs over all fast ion species
  !
  bperr=0.0_fp
  bpherr=0.0_fp
  !
  do isb=1,nsbeam
     !
     !  renormalize quantities summed in orball subroutine
     !  compute gain terms
     !
     !  energy
     enolds(isb)=enolds(isb)*xe
     bpsts(isb)=(enolds(isb)-ebcons(isb))*zdti
     ebcons(isb)=enolds(isb)
     !  momentum
     bpholds(isb)=bpholds(isb)*1.0e-7_fp
     bphsts(isb)=(bpholds(isb)-bphcons(isb))*zdti
     bphcons(isb)=bpholds(isb)
     !
     !  following quantities:  see also outptb
     !
     !  orbit code error checks
     bpherr=bpherr+bpherrs(isb)
     bperr=bperr+bperrs(isb)
     !  momentum balance check
     bphcks(isb)=bphdeps(isb)+bphohbs(isb) &
          -bph0was(isb)-bphjxbs(isb)-bphcols(isb)-bphdfbs(isb) &
          -bphths(isb)-bphorbs(isb)-bphsts(isb)-bphfshs(isb)
     !

  end do
  !
  !  zero neutral density sums...
  !
  bn0tot=0.0_fp; bn00t=0.0_fp; bn00ts=0.0_fp
  !
  !  zero out beam target fusion total
  btneut=0.0_fp
  rbneut=0.0_fp

  !
  !--------------------------------------------------------------------
  !  profile sums
  !
  !
  bdens_nb = 0.0_fp
  pbe_nb = 0.0_fp
  pbi_nb = 0.0_fp
  pbtot_nb = 0.0_fp
  bdepb_nb = 0.0_fp
  tqtot_nb = 0.0_fp
  !
  do j=lcentr,ledge
     jp1=j+1
     !
     ztqbi=0.0_fp
     pbi(j)=0.0_fp
     pbe(j)=0.0_fp
     pfi(j)=0.0_fp
     pfe(j)=0.0_fp
     pbcprs(j)=0.0_fp
     pbth(j)=0.0_fp
     sorbe(j)=0.0_fp
     sorbth(j)=0.0_fp
     !
     tqbcx(j)=0.0_fp
     tqjxb(j)=0.0_fp
     !
     tqbi(j)=0._fp
     tqbe(j)=0._fp
     tqbth(j)=0._fp
     tqbxf(j)=0._fp
     tqrpl(j)=0._fp
     tqohb(j)=0._fp
     tqdifb(j)=0.0_fp
     !
     tirorb(jp1)=0.0_fp
     tirabd(jp1)=0.0_fp
     tirfsh(jp1)=0.0_fp
     tirrip(jp1)=0.0_fp
     !
     do ig=1,ng
        sorbh(j,ig)=0.0_fp
        sorbn0(j,ig)=0.0_fp
        sinbn0(j,ig)=0.0_fp
        sinb0i(j,ig)=0.0_fp
     end do
     !  pbe,pbi zeroed in outptb
     do isb=1,nsbeam
        !
        if(nlfprod(isb)) then
           pfi(j)=pfi(j)+pbis(j,isb)+pbdfbis(j,isb)
           pfe(j)=pfe(j)+pbes(j,isb)+pbdfbes(j,isb)
        else
           pbi(j)=pbi(j)+pbis(j,isb)+pbdfbis(j,isb)
           pbe(j)=pbe(j)+pbes(j,isb)+pbdfbes(j,isb)
        end if
        pbcprs(j)=pbcprs(j)+pbcprss(j,isb)
        pbth(j)=pbth(j)+pbths(j,isb)+pbthr2s(j,isb)
        !
        sorbe(j)=sorbe(j)+sorbes(j,isb)
        sorbth(j)=sorbth(j)+sorbths(j,isb)
        !
        do ig=1,ng
           sorbh(j,ig)=sorbh(j,ig)+sorbhs(j,ig,isb)
           sorbn0(j,ig)=sorbn0(j,ig)+sorbn0s(j,ig,isb)
           ! dmc:  convert sinbn0 & sinb0i to <nb*sigma*v> (1/sec)
           !   by dividing out the neutral density (if there is no neutral
           !   density the item must be zero).  also divide out the volume!
           zn0=walln0x(ig,j)+voln0(ig,j)
           if(zn0.gt.0.0_fp) then
              sinbn0s(j,ig,isb)=sinbn0s(j,ig,isb)/(dvolzc(j)*zn0)
              sinb0is(j,ig,isb)=sinb0is(j,ig,isb)/(dvolzc(j)*zn0)
           else
              sinbn0s(j,ig,isb)=0.0_fp
              sinb0is(j,ig,isb)=0.0_fp
           end if
           sinbn0(j,ig)=sinbn0(j,ig)+sinbn0s(j,ig,isb)
           sinb0i(j,ig)=sinb0i(j,ig)+sinb0is(j,ig,isb)

        end do
        !
        tqohb(j)=tqohb(j)+tqohbs(j,isb)
        tqbxf(j)=tqbxf(j)+tqbxfs(j,isb)*abs(tqbis(j,isb))
        ztqbi=ztqbi+abs(tqbis(j,isb))
        tqbth(j)=tqbth(j)+tqbths(j,isb)
        tqbi(j)=tqbi(j)+tqbis(j,isb)
        tqbe(j)=tqbe(j)+tqbes(j,isb)

        tqjxb(j)=tqjxb(j)+tqjxbs(j,isb)
        tqrpl(j)=tqrpl(j)+tqrpls(j,isb)
        tqbcx(j)=tqbcx(j)+tqbcxs(j,isb)
        tqdifb(j)=tqdifb(j)+tqdifbs(j,isb)
        !
        tqbsum(j)=tqbe(j)+tqbi(j)+tqbth(j)+tqjxb(j)+tqbcx(j)+tqdifb(j)
        !
        !  radial current sums
        !
        tirorb(jp1)=tirorb(jp1)+tirorbs(jp1,isb)
        tirabd(jp1)=tirabd(jp1)+tirabds(jp1,isb)
        tirfsh(jp1)=tirfsh(jp1)+tirfshs(jp1,isb)
        tirrip(jp1)=tirrip(jp1)+tirrips(jp1,isb)
        !
        !  normalize nth generartion neutral density and
        !	avg energy into past timestep array
        !
        if(bn0cx(j,isb).gt.0._fp) then
           zv2av=be0cx(j,isb)/bn0cx(j,isb)
           zvzav=bv0cx(j,isb)/bn0cx(j,isb)
           be0cx(j,isb)=(zv2av-zvzav*zvzav)*abeams(isb)/evptocm_sec2
           bv0cx(j,isb)=zvzav
           !  excited states stuff
           if(nsigexc_use.gt.0) then
              excitr(j,isb)=excitr(j,isb)/bn0cx(j,isb)
           else
              excitr(j,isb)=1.0_fp
           end if
           !
           bn0cx(j,isb)=bn0cx(j,isb)*zdti
        else
           bv0cx(j,isb)=0.0_fp
           be0cx(j,isb)=0.0_fp
        end if
        !

     end do
     !

     if(ztqbi.gt.0.0_fp) then
        tqbxf(j)=tqbxf(j)/ztqbi
     else
        tqbxf(j)=0.0_fp
     end if
     !
     !	beam-beam reaction rate arrays
     !
     do ien=1,3
        do ib=1,nbeam
           bn00t(ien,j)=bn00t(ien,j)+bn00(ib,ien,j)
           ! tot 1st gen neutral density by beam species -- for output
           do isb=1,nsbeam
              if(lsbeama(ib).eq.isb) then

                 bn0tot(j,isb)=bn0tot(j,isb)+bn00(ib,ien,j)
                 bn00ts(ien,j,isb)=bn00ts(ien,j,isb)+bn00(ib,ien,j)

              end if

           end do

        end do

     end do
     !
     !  beam target fusion products output
     !    (see bbneut for beam beam fusion)
     !
     btntnt(j)=0.0_fp
     rbntnt(j)=0.0_fp
     do ireact=1,itreact
        if(btntws(j,ireact).gt.0.0_fp) then
           btntrs(j,ireact)=btntrs(j,ireact)/btntws(j,ireact)
        else
           btntrs(j,ireact)=0.0_fp
        end if
        btntns(j,ireact)=btntns(j,ireact)*zdti
        btnta(ireact)=btnta(ireact)+btntns(j,ireact)
        if(ireact.le.8) then
           btntnt(j)=btntnt(j)+ineuts(ireact)*btntns(j,ireact)
           btneut=btneut+ineuts(ireact)*btntns(j,ireact)
        end if
        if(ireact.ge.11) then
           rbntnt(j)=rbntnt(j)+ineuts(ireact)*btntns(j,ireact)
           rbneut=rbneut+ineuts(ireact)*btntns(j,ireact)
        end if


     end do
     !
     !----------------------------------------------------------------------
     !  beam driven current
     !  dmc generalized apr 1985 - now evaluating <j dot b>driven * dv
     !
     !  d. mc cune impliment neoclassical beam current option
     !    collisionality, trapping fraction formulae from "resis.for"
     !    a,b, k0 coefficients from hirschman, phys. fluids vol 21
     !    no. 8 aug 1978.
     !	put into fokker.for, for self-consistent e-field effect
     !	on fast ions    --  rjg dec. 27, 1982
     !	 --> nmcurb=0 ==> no beam current
     !
     !
     !	note that nmcurb=0 gives inconsistent oh power balance
     !	since some of i*v is going to beam ions, but magdif won't know
     !
     pljbbm(j)=0.0_fp
     pljbfi(j)=0.0_fp
     !
     if(nmcurb.gt.0) then
        !
        !  this output has volume normalization removed, in case
        !  plasma changes size  ** to get tesla*amps/cm2, divide by dvol **
        !
        !  **see subroutine mgpout, magnetics library**
        !  where the toroidal driven current density is computed
        !
        zbzxr = fbz(j,2)*bzxr
        zr2i = gr2izc(j)
        zri = grizc(j)
        zb2 = gb2zc(j)
        !
        xjbfac(j)=0
        zwsum=0
        zsumjbco=0
        zsumjb=0
        do isb=1,nsbeam
           ucurjs(j,isb)= &
                get_j_from_jdotb(zbzxr,zr2i,zb2,zri,pljbbms(j,isb))
           curjs(j,isb)=ucurjs(j,isb)*xjbfacs(j,isb) ! shielded

           zws1=max(1.0e-50_fp,abs(ucurjs(j,isb)))
           zwsum=zwsum+zws1
           xjbfac(j)=xjbfac(j)+zws1*xjbfacs(j,isb)

           if(nlfprod(isb)) then
              pljbfi(j)=pljbfi(j)+pljbbms(j,isb)*xjbfacs(j,isb) ! shielded
           else
              zsumjb=zsumjb+pljbbms(j,isb)
              zsumjbco=zsumjbco+pljbcct(j,isb)
              pljbbm(j)=pljbbm(j)+pljbbms(j,isb)*xjbfacs(j,isb) ! shielded
           end if

        end do
        xjbfac(j)=xjbfac(j)/max(1.0e-50_fp,zwsum)
        !
        curjbeam(j)=get_j_from_jdotb(zbzxr,zr2i,zb2,zri,pljbbm(j))
        curjfusi(j)=get_j_from_jdotb(zbzxr,zr2i,zb2,zri,pljbfi(j))
        !
        ucurjco(j)= &
             get_j_from_jdotb(zbzxr,zr2i,zb2,zri,zsumjbco)
        !! dv * j dot b (unshielded, driven by co sces) -> da*j
        ucurjctr(j)= &
             get_j_from_jdotb(zbzxr,zr2i,zb2,zri,zsumjb-zsumjbco)
        !! dv * j dot b (unshielded, driven by ctr sces) -> da*j
        !
     end if
     !
     !  loop over beams & energy fractions
     !
     do ib=1,nbeam
        do ien=1,3
           znt=bnb_ien(j,ien,ib)
           if(znt.gt.0.0_fp) then
              upb_ien(j,ien,ib)=upb_ien(j,ien,ib)/znt
           end if
           bnb_ien(j,ien,ib)=bnb_ien(j,ien,ib)*zdti
           pbe_ien(j,ien,ib)=pbe_ien(j,ien,ib)*zbpf
           pbi_ien(j,ien,ib)=pbi_ien(j,ien,ib)*zbpf
           !
           bdens_nb(j,ib)=bdens_nb(j,ib)+bnb_ien(j,ien,ib)
           bdepb_nb(j,ib)=bdepb_nb(j,ib)+bdepb(j,ien,ib)
           pbe_nb(j,ib)=pbe_nb(j,ib)+pbe_ien(j,ien,ib)
           pbi_nb(j,ib)=pbi_nb(j,ib)+pbi_ien(j,ien,ib)
           !
        end do
        !
        pbth_nb(j,ib) = pbth_nb(j,ib)*zbpf
        pbtot_nb(j,ib)=pbe_nb(j,ib)+pbi_nb(j,ib)+pbth_nb(j,ib)
        !  dv*<j dot b> -> da*<j dot b>
        pljbbm_nb(j,ib)=pljbbm_nb(j,ib)*zbpf* &
             dareazc_r8(j)/dvolzc_r8(j)

        if(nmcurb.gt.0) then
           ucurjnb(j,ib)= &   !unshielded 
                get_j_from_jdotb(zbzxr,zr2i,zb2,zri,pljbbm_nb(j,ib))
        end if
        xjbfacnb=0.0_fp
        do isb=1,nsbeam
           if(.not.nlfprod(isb).and.lsbeama(ib).eq.isb)then
              xjbfacnb=xjbfacs(j,isb)
           end if
        end do

        if(nmcurb.gt.0) then
           curjnb(j,ib)=ucurjnb(j,ib)*xjbfacnb !shielded
        end if


     end do

     !
     !  end of 1d radial profile renorm (j) loop
     !

  end do
  sdep_oauc(:) = sdep_oa(lep1,:) !orbit averaged unconfined"
  ! torque normalization
  tqjxb_nb=(tqjxb_nb+tqjxbnb_ntr)*zdti*1.0e-7_fp
  tqjxb_nb_dep=tqjxb_nb_dep*zdti*1.0e-7_fp
  tqbco_nb=tqbco_nb*zdti*1.0e-7_fp
  tqbth_nb=tqbth_nb*zdti*1.0e-7_fp
  tqjxb_nb=tqjxb_nb+tqjxb_nb_dep
  tqtot_nb=tqbco_nb+tqjxb_nb+tqbth_nb

  !
  !  pseudo-density and dn/dt outside the plasma.  reusing defunct common
  !  variables xysbcx/xysbcxs to hold this output calculation -- used in
  !  the global fast ion particle balance (dmc feb 1992)
  do isb=1,nsbeam
     xysbcxs(isb)=(bdenss(lep1,isb)-bdenssp(lep1,isb))*zdti
  end do
  !
  !----------------------------------------------------
  !  compute beam-beam fusion neutrons
  !  (beam dist fcns have been flipped and normalized in outptb + flipr)
  !

  if(nlfbm) then
     call bbneut(ierr)
     if(ierr.ne.0) return
  end if

  !
  !--------------------------------------
  !  2d profile data normalization:  dmc mar 1985



  ! calculate bdens2, eba2pl, eba3pp, vbtor2 at gc
  call param_fbm


  !
  jl(1:nsbeam)=0
  do iz=1,nfbzns
     !  average perp and pll energies
     do isb=1,nsbeam
        ! at guiding center
        if(bdens2(iz,isb).gt.0.0_fp) then
           vbtor2(iz,isb)=vbtor2(iz,isb)*btiltc(2,iz)/bdens2(iz,isb)
           eba2pp(iz,isb)=eba2pp(iz,isb)/bdens2(iz,isb)
           eba2pl(iz,isb)=eba2pl(iz,isb)/bdens2(iz,isb)
        else

           jl(isb)=jl(isb)+1

           vbtor2(iz,isb)=0.0_fp
           eba2pp(iz,isb)=0.0_fp
           eba2pl(iz,isb)=0.0_fp
        end if
        bdens2(iz,isb)=bdens2(iz,isb)*zdti

        !
        vbtr2p(iz,isb)=vbtor2(iz,isb)  ! for nbatom, next step
        eb2ppp(iz,isb)=eba2pp(iz,isb)  ! ditto
        eb2plp(iz,isb)=eba2pl(iz,isb)  ! ditto
        bdns2p(iz,isb)=bdens2(iz,isb)  ! ditto
        !
        !  unshielded current density *dv
        xcurb2(iz,isb)=xcurb2(iz,isb)*xe*zdti
        xcrb2c(iz,isb)=xcrb2c(iz,isb)*xe*zdti
        !
        !  sce fcns
        do ig=1,ng
           srbn02(iz,ig,isb)= &
                (srbn02(iz,ig,isb)+srbn02_dep(iz,ig,isb))*zdti ! ptcls/sec/zone
           srbn02_dep(iz,ig,isb)=srbn02_dep(iz,ig,isb)*zdti
           ! dmc:  convert sinbn0 & sinb0i to <nb*sigma*v> (1/sec)
           !   by dividing out the neutral density (if there is no neutral
           !   density the item must be zero).  also divide out the volume!
           zn0=walln02(ig,iz)+voln02(ig,iz)
           if(zn0.gt.0.0_fp) then
              snbn02(iz,ig,isb)=snbn02(iz,ig,isb)*zdti/(zn0*bmvol(iz))
              snb0i2(iz,ig,isb)=snb0i2(iz,ig,isb)*zdti/(zn0*bmvol(iz))
           else
              snbn02(iz,ig,isb)=0.0_fp
              snb0i2(iz,ig,isb)=0.0_fp
           end if

        end do

        !  nth generation cx neutral density, energy
        if(bn0cx2(iz,isb).gt.0.0_fp) then
           zv2av=be0cx2(iz,isb)/bn0cx2(iz,isb)
           zvzav=bv0cx2(iz,isb)/bn0cx2(iz,isb)
           bv0cx2(iz,isb)=zvzav
           bvrhocx2(iz,isb)=bvrhocx2(iz,isb)/bn0cx2(iz,isb)
           bvthcx2(iz,isb)=bvthcx2(iz,isb)/bn0cx2(iz,isb)
           bvpllcx2(iz,isb)=bvpllcx2(iz,isb)/bn0cx2(iz,isb)
           balpha_gyro(iz,isb)=balpha_gyro(iz,isb) &
                /bn0cx2(iz,isb)
           be0cx2(iz,isb)=(zv2av-zvzav*zvzav-bvrhocx2(iz,isb)* &
                bvrhocx2(iz,isb)-bvthcx2(iz,isb)*bvthcx2(iz,isb)) &
                *abeams(isb)/evptocm_sec2
           be0cx2(iz,isb)=max(0.0_fp,be0cx2(iz,isb))
           !  excited state stuff
           if(nsigexc_use.gt.0) then
              do isb2=1,nsbeam
                 hexbbxrc(iz,isb,isb2)= &
                      hexbbxrc(iz,isb,isb2)/bn0cx2(iz,isb)
                 hexbbxrp(iz,isb,isb2)=hexbbxrc(iz,isb,isb2)
              end do
           else
              do isb2=1,nsbeam
                 hexbbxrc(iz,isb,isb2)=1.0_fp
                 hexbbxrp(iz,isb,isb2)=hexbbxrc(iz,isb,isb2)
              end do
           end if
           bn0cx2(iz,isb)=bn0cx2(iz,isb)*zdti
        else
           !  excited state stuff
           do isb2=1,nsbeam
              hexbbxrc(iz,isb,isb2)=1.0_fp
              hexbbxrp(iz,isb,isb2)=1.0_fp
           end do
           bv0cx2(iz,isb)=0.0_fp
           be0cx2(iz,isb)=0.0_fp
        end if
        !

     end do

     !
     !  beam target neutrons -- 2d grid --
     !
     do ireact=1,nreact
        btntn2(iz,ireact)=btntn2(iz,ireact)*zdti
     end do
     !

  end do

  !
  !  1st gen neutrals (normalized in depnrm) and reaction rate coeffs
  ! sn00i2 and sn00x2 are  <nb*sigma*v> (1/sec)
  do isb=1,nsbeam
     do iz=1,nfbznsi
        do ien=1,3
           do ib=1,nbeam
              sn00i2(ib,ien,iz,isb)=sn00i2(ib,ien,iz,isb)*zdti
              sn00x2(ib,ien,iz,isb)=sn00x2(ib,ien,iz,isb)*zdti
           end do
        end do
     end do
  end do
  if(levmod_halo.eq.2) then
     svn3vcx_h=svn3vcx_h*zdti
     svn3vcx_he=svn3vcx_he*zdti
     svn3vii_h=svn3vii_h*zdti
     svn3vii_he=svn3vii_he*zdti
  end if
  !
  !  dmc 18 sep 1995:  trapping fraction calculations
  !
  do isb=1,nsbeam
     !  average over profiles, whole population
     zwnorm=bsumint(wtrap(1:mj,isb))
     if(zwnorm.gt.0.0_fp) then
        tbtrap(isb)=bsumint(btrap(1:mj,isb))/zwnorm
     end if
     do j=lcentr,lep1
        !  for whole population
        if(wtrap(j,isb).eq.0.0_fp) then
           btrap(j,isb)=0.0_fp
        else
           btrap(j,isb)=btrap(j,isb)/wtrap(j,isb)
        end if
        !  for energy bins
        do ie=1,nerngfi
           if(wtrape(j,ie,isb).eq.0.0_fp) then
              btrape(j,ie,isb)=0.0_fp
           else
              btrape(j,ie,isb)=btrape(j,ie,isb)/wtrape(j,ie,isb)
           end if
        end do               ! ie
     end do                          ! j
  end do                             ! isb
  !
  ! normalization of the average relative energy [ev] for beam-beam reaction
  !co-injected beam populations
  do isb=1,nsbeam
     do j=1,mimxbzf
        if(bnb_ien2(j,1,isb).gt.0.0_fp) then
           do ib=1,nbeam
              do ien=1,3
                 ebb_rel_co(j,ien,ib,isb)= &
                      ebb_rel_co(j,ien,ib,isb)/ &
                      bnb_ien2(j,1,isb)
                 ebb_rel2_co(j,ien,ib,isb)= &
                      ebb_rel2_co(j,ien,ib,isb)/ &
                      bnb_ien2(j,1,isb)
                 rms_ebb_co(j,ien,ib,isb)= sqrt(max(0._fp, &
                      ebb_rel2_co(j,ien,ib,isb)- &
                      ebb_rel_co(j,ien,ib,isb)* &
                      ebb_rel_co(j,ien,ib,isb)))
              end do
           end do
        else
           do ib=1,nbeam
              do ien=1,3
                 ebb_rel_co(j,ien,ib,isb)=0.0_fp
                 ebb_rel2_co(j,ien,ib,isb)=0.0_fp
                 rms_ebb_co(j,ien,ib,isb)=0.0_fp
              end do
           end do
        end if
     end do
  end do
  !     counter-injected beam populations
  do isb=1,nsbeam
     do j=1,mimxbzf
        if(bnb_ien2(j,2,isb).gt.0.0_fp) then
           do ib=1,nbeam
              do ien=1,3
                 ebb_rel_cnt(j,ien,ib,isb)= &
                      ebb_rel_cnt(j,ien,ib,isb)/ &
                      bnb_ien2(j,2,isb)
                 ebb_rel2_cnt(j,ien,ib,isb)= &
                      ebb_rel2_cnt(j,ien,ib,isb)/ &
                      bnb_ien2(j,2,isb)
                 rms_ebb_cnt(j,ien,ib,isb)= sqrt(max(0._fp, &
                      ebb_rel2_cnt(j,ien,ib,isb)- &
                      ebb_rel_cnt(j,ien,ib,isb)* &
                      ebb_rel_cnt(j,ien,ib,isb)))
              end do
           end do
        else
           do ib=1,nbeam
              do ien=1,3
                 ebb_rel_cnt(j,ien,ib,isb)=0.0_fp
                 ebb_rel2_cnt(j,ien,ib,isb)=0.0_fp
                 rms_ebb_cnt(j,ien,ib,isb)=0.0_fp
              end do
           end do
        end if
     end do
  end do
  !fusioun product populations
  do isb=1,nsbeam
     do j=1,mimxbzf
        if(bnb_fn(j,isb).gt.0.0_fp) then
           do ib=1,nbeam
              do ien=1,3
                 ebb_rel_fn(j,ien,ib,isb)= &
                      ebb_rel_fn(j,ien,ib,isb)/ &
                      bnb_fn(j,isb)
                 ebb_rel2_fn(j,ien,ib,isb)= &
                      ebb_rel2_fn(j,ien,ib,isb)/ &
                      bnb_fn(j,isb)
                 rms_ebb_fn(j,ien,ib,isb)= sqrt(max(0._fp, &
                      ebb_rel2_fn(j,ien,ib,isb)- &
                      ebb_rel_fn(j,ien,ib,isb)* &
                      ebb_rel_fn(j,ien,ib,isb)))
              end do
           end do
        else
           do ib=1,nbeam
              do ien=1,3
                 ebb_rel_fn(j,ien,ib,isb)=0.0_fp
                 ebb_rel2_fn(j,ien,ib,isb)=0.0_fp
                 rms_ebb_fn(j,ien,ib,isb)=0.0_fp
              end do
           end do
        end if
     end do
  end do
  ! fast ion density normalization
  bnb_ien2=bnb_ien2*zdti
  bnb_fn=bnb_fn*zdti
  !c
  !----------------------------------------------
  !  patch holes in some output
  !

  call orbsmoo2


  !
  !--------------------------------------------------------------------
  !  save 2d grid quantities via xplasma, for use by callers
  !
  !  in transp, these 2d grid profiles are received in
  !  source/trfpi/fi_finish.for -- changes should be synchronized
  !  between that routine and this section of outptb2 (dmc dec. 2001)
  !
  !--------------------------------------------------------------------
  if(norbtot.eq.0) goto 1000       ! skip if no orbits calculated.
  !
  !cc      call xplasma_author_set(s,'nubeam',ierr) -- handled in nbfinish now!
  !
  !  retrieve or save pa grid for fbm
  !
  call xplasma_gridid(s,'FBM_VPVGRID',id_vpv)
  if(id_vpv.eq.0) then
     zxpa(1)=-1.0_fp
     zxpa(nznbma+1)=1.0_fp
     do i=1,nznbma-1
        zxpa(i+1)=-1.0_fp + i*2.0_fp/nznbma
     end do
     call xplasma_create_grid(s,'FBM_VPVGRID',xplasma_vpv_coord, &
          zxpa,id_vpv,ierr)
     if(ierr.ne.0) then
        write(nonlin,*) ' ?unexpected xplasma error in outptb2:'
        write(nonlin,*) '  vpll/v grid creation attempt failed.'
        call xplasma_error(s,ierr,nonlin)
        goto 1000
     end if
  end if
  !
  do isb=1,inum
     zlbl=zlbls(isb)
     ilbl=ilbls(isb)
     if((index(zlbl,'_MCBEAM ').gt.0).or. &
          (index(zlbl,'_MCFUSN ').gt.0)) then
        !  save mc gridded output
        zconv=1.0_fp
        do ig=1,ng
           !
           !  neutral source (#/sec, each zone)
           !  example:  "SRBN02_D0_H_MCBEAM" -- d thermal neutral source due to
           !             h neutral beam
           !
           zlblobj='SRBN02_'//zlbth(ig)(1:ilbth(ig))//'0_'// &
                zlbl(1:ilbl)
           call mklabel
           call xplasma_mcgrid_putobj( &
                s,id_mcgrid,zlblobj,idpro,jerr, &
                data_1d=srbn02(1:mimxbzf,ig,isb), &
                label=label,units='#/sec',normcode=-1)
           call ckerr(jerr,zlblobj)
           !
           !  neutral sink coefficients:  cx then ii  <nb*sigma*v>  (1/sec)
           !
           zlblobj='SNBN02_'//zlbth(ig)(1:ilbth(ig))//'0_'// &
                zlbl(1:ilbl)
           call mklabel
           call xplasma_mcgrid_putobj( &
                s,id_mcgrid,zlblobj,idpro,jerr, &
                data_1d=snbn02(1:mimxbzf,ig,isb), &
                label=label,units='1/sec',normcode=0)
           call ckerr(jerr,zlblobj)
           !
           zlblobj='SNB0I2_'//zlbth(ig)(1:ilbth(ig))//'0_'// &
                zlbl(1:ilbl)
           call mklabel
           call xplasma_mcgrid_putobj( &
                s,id_mcgrid,zlblobj,idpro,jerr, &
                data_1d=snb0i2(1:mimxbzf,ig,isb), &
                label=label,units='#/sec',normcode=0)
           call ckerr(jerr,zlblobj)
        end do
        !
        !  no. of ptcls per 2d zone (#, each zone)
        !    (also put data for ptcls beyond the plasma; no volume normalization
        !    is available, however).
        !
        zlblobj='BDENS2_'//zlbl(1:ilbl)
        call mklabel
        call xplasma_mcgrid_putobj( &
             s,id_mcgrid,zlblobj,idpro,jerr, &
             data_1d=bdens2(1:mimxbzf,isb), &
             label=label,units='#',normcode=-1)
        call ckerr(jerr,zlblobj)
        !
        !  no. of neutral ptcls from cx per 2d zone (#, each zone)
        !    (also put data for ptcls beyond the plasma; no volume normalization
        !    is available, however).
        !
        zlblobj='BN0CX2_'//zlbl(1:ilbl)
        call mklabel
        call xplasma_mcgrid_putobj( &
             s,id_mcgrid,zlblobj,idpro,jerr, &
             data_1d=bn0cx2(1:mimxbzf,isb), &
             label=label,units='#',normcode=-1)
        call ckerr(jerr,zlblobj)
        !
        !  average energy perp, pll; convert to kev
        !
        zconv=0.001_fp
        !  <eperp>
        zlblobj='EBA2PP_'//zlbl(1:ilbl)
        call mklabel
        call xplasma_mcgrid_putobj( &
             s,id_mcgrid,zlblobj,idpro,jerr, &
             data_1d=zconv*eba2pp(1:mimxbzf,isb), &
             label=label,units='KeV',normcode=0)
        call ckerr(jerr,zlblobj)
        !  <epll>
        zlblobj='EBA2PL_'//zlbl(1:ilbl)
        call mklabel
        call xplasma_mcgrid_putobj( &
             s,id_mcgrid,zlblobj,idpro,jerr, &
             data_1d=zconv*eba2pl(1:mimxbzf,isb), &
             label=label,units='KeV',normcode=0)
        call ckerr(jerr,zlblobj)
        !  <eneutral>
        zlblobj='BE0CX2_'//zlbl(1:ilbl)
        call mklabel
        call xplasma_mcgrid_putobj( &
             s,id_mcgrid,zlblobj,idpro,jerr, &
             data_1d=zconv*be0cx2(1:mimxbzf,isb), &
             label=label,units='KeV',normcode=0)
        call ckerr(jerr,zlblobj)
        !
        !  average fast ion toroidal velocity; convert to m/s
        !
        zconv=0.01_fp
        !  <vtor>
        zlblobj='VBTOR2_'//zlbl(1:ilbl)
        call mklabel
        call xplasma_mcgrid_putobj( &
             s,id_mcgrid,zlblobj,idpro,jerr, &
             data_1d=zconv*vbtor2(1:mimxbzf,isb), &
             label=label,units='m/sec',normcode=0)
        call ckerr(jerr,zlblobj)
        !  <vneutral>
        zlblobj='BV0CX2_'//zlbl(1:ilbl)
        call mklabel
        call xplasma_mcgrid_putobj( &
             s,id_mcgrid,zlblobj,idpro,jerr, &
             data_1d=zconv*bv0cx2(1:mimxbzf,isb), &
             label=label,units='m/sec',normcode=0)
        call ckerr(jerr,zlblobj)
        !
        zconv=1.00_fp
        !
        !  unshielded current (amps, each zone)
        !
        zlblobj='XCURB2_'//zlbl(1:ilbl)
        call mklabel
        call xplasma_mcgrid_putobj( &
             s,id_mcgrid,zlblobj,idpro,jerr, &
             data_1d=xcurb2(1:mimxbzf,isb), &
             label=label,units='A',normcode=0)
        call ckerr(jerr,zlblobj)
        !
        !  unshielded current (amps, each zone) by co-injected ions
        !
        zlblobj='XCRB2C_'//zlbl(1:ilbl)
        call mklabel
        call xplasma_mcgrid_putobj( &
             s,id_mcgrid,zlblobj,idpro,jerr, &
             data_1d=xcrb2c(1:mimxbzf,isb), &
             label=label,units='A',normcode=0)
        call ckerr(jerr,zlblobj)
        !
        !  gc distribution function...
        !  distribution function...
        !    units: #/(dsa)/ev/cm^3  --> #/(dsa)/kev/m^3)  (mult. by 1.0e9 for mks)
        !    (dsa) means "delta-solid-angle" which works out to 1/nznbma
        !      for cos(pitch-angle) grid [vpll/v= -1 to +1] evenly divided into
        !      nznbma zones
        !
        zconv=1.0e9_fp
        zconve=1.0e-3_fp
        zlblobj='FBM_'//zlbl(1:ilbl)
        call mklabel

        call nbi_put_fbm(id_mcgrid,zlblobj,label,id_vpv, &
             zconve,efbmb(1,isb),zconv,fbm(1,1,1,isb), &
             nznbmea(isb),nznbma,mimxbe,mimxba,mimxbzf,jerr)

        call ckerr(jerr,zlblobj)
        !
     end if                          ! test if it is an mc ion specie
  end do                             ! loop over species
  !
  !    units: #/(dsa)/ev/cm^3  --> #/(dsa)/kev/m^3)  (mult. by 1.0e9 for mks)
  !    (dsa) means "delta-solid-angle" which works out to 1/nznbma
  !      for cos(pitch-angle) grid [vpll/v= -1 to +1] evenly divided into
  !      nznbma zones
  if(nlfbmflr.and.nltest_output)then
     fname=trim(workpath)//trim(runid)//'_fbm_ptcl.cdf'
     !overwrite netcdf  file for fbm_ptcl distribution function
     ierr = nf90_create(path=trim(fname),cmode= or(nf90_clobber, &
          nf90_64bit_offset),ncid=icdf)

     if(ierr.ne.0) then
        write(nonlin,*) &
             ' ?outptb2 error: cannot open file: '//trim(fname)
        return
     end if
     call check(nf90_def_dim(icdf,'mimxbe',mimxbe,id_mimxbe), &
          'mimxbe')
     call check(nf90_def_dim(icdf,'mimxba',mimxba,id_mimxba), &
          'mimxba')
     call check(nf90_def_dim(icdf,'mimxbzf',mimxbzf,id_mimxbzf), &
          'mimxbzf')
     call check(nf90_def_dim(icdf,'mgyro',mgyro,id_mgyro), &
          'mgyro')
     i_sb=0
     zconv=1.0e9_fp
     zconve=1.0e-3_fp
     if(allocated(id_fbm_ptcl))deallocate(id_fbm_ptcl)
     allocate(id_fbm_ptcl(inum))
     id_fbm_ptcl=0
     if(allocated(id_e))deallocate(id_e)
     allocate(id_e(inum))
     id_e=0
     do isb=1,inum
        zlbl=zlbls(isb)
        ilbl=ilbls(isb)
        if((index(zlbl,'_MCBEAM ').gt.0).or. &
             (index(zlbl,'_MCFUSN ').gt.0)) then
           i_sb=i_sb+1

           zlblobj='E_'//zlbl(1:ilbl)
           call check( nf90_def_var(icdf,zlblobj,nf90_double, &
                (/id_mimxbe/),id_e(i_sb)),zlblobj)
           call check( nf90_put_att(icdf,id_e(i_sb),zlblobj, &
                'fast ion distribution energy table for'// &
                trim(zlbl(1:ilbl))//',KeV'),zlblobj)


           zlblobj='FBM_PTCL_'//zlbl(1:ilbl)
           call check( nf90_def_var(icdf,zlblobj,nf90_double, &
                (/id_mimxbe,id_mimxba,id_mimxbzf,id_mgyro/), &
                id_fbm_ptcl(i_sb)),zlblobj)

           call check( nf90_put_att(icdf,id_fbm_ptcl(i_sb),zlblobj, &
                'flr fast ion distribution for '//trim(zlbl(1:ilbl)) &
                //',atom*m:-3*KeV:-1'),zlblobj)

           !
        end if               ! test if it is an mc ion specie
     end do                  ! loop over species
     ierr = nf90_enddef(icdf) ! out of def mode and goto data mode
     i_sb=0
     do isb=1,inum
        zlbl=zlbls(isb)
        ilbl=ilbls(isb)
        if((index(zlbl,'_MCBEAM ').gt.0).or. &
             (index(zlbl,'_MCFUSN ').gt.0)) then

           i_sb=i_sb+1
           zlblobj='E_'//zlbl(1:ilbl)
           call check( nf90_put_var(icdf,id_e(i_sb), &
                efbm(1:mimxbe,isb)*zconve),zlblobj)

           zlblobj='FBM_PTCL_'//zlbl(1:ilbl)
           call check( nf90_put_var(icdf,id_fbm_ptcl(i_sb), &
                fbm_ptcl(1:mimxbe,1:mimxba,1:mimxbzf,isb,1:mgyro)* &
                zconv),zlblobj)
        end if               ! test if it is an mc ion specie
     end do                  ! loop over species

     ierr = nf90_close( icdf)

     if(ierr.ne.0) then
        write(nonlin,*) ' ?outptb2: netcdf close failure: ', &
             trim(fname)
        return
     end if
     if(allocated(id_fbm_ptcl))deallocate(id_fbm_ptcl)
     if(allocated(id_e))deallocate(id_e)
  end if !nlfbmflr.and.nltest_output
  !
  !  fusion rate data on 2d grids  (#/sec)
  !

  zconv=1.0_fp
  do ireact=1,nreact

     call catfus_mtype(ireact,ityp)  ! 1:beam, 2:fusion-product, 3:rf tail

     if(sftot(ireact).gt.0.0_fp) then

        ! 2d fusion source profile
        ! for some reaction codes this includes a thermonuclear
        ! contribution; beam-beam and beam-target contributions are
        ! from preceding time step.  renormalization to data may have
        ! occurred.

        write(zrnum,'(i2.2)') ireact
        label = 'fusion source: catfus no. '//zrnum
        zlblobj = 'sfusion_'//zrnum
        call xplasma_mcgrid_putobj( &
             s,id_mcgrid,zlblobj,idpro,jerr, &
             data_1d=sfusion(1:mimxbzf,ireact), &
             label=label,units='#/sec',normcode=-1)
        call ckerr(jerr,zlblobj)
     end if

     if(sum(btntn2(1:nfbzns,ireact)).gt.0.0_fp) then

        ! thermal target

        call catfus_name(ireact,ityp,0,zlbl) ! beam-target
        if(zlbl.eq.'error') then
           write(nonlin,*) ' ?outptb2: fusion data handling error'
           ierr=1
        else
           label='fast ion -> target fusion: '//trim(zlbl)
           zlblobj=trim(zlbl)//'_2D'
           call xplasma_mcgrid_putobj( &
                s,id_mcgrid,zlblobj,idpro,jerr, &
                data_1d=btntn2(1:mimxbzf,ireact), &
                label=label,units='#/sec',normcode=-1)
           call ckerr(jerr,zlblobj)
        end if

     end if

     if(sum(bbnvol(1:nfbzns,ireact)).gt.0.0_fp) then

        ! fast-ion - fast-ion (beam) reaction

        call catfus_name(ireact,ityp,1,zlbl) ! beam-beam
        if(zlbl.eq.'error') then
           write(nonlin,*) ' ?outptb2: fusion data handling error'
           ierr=1
        else
           label='fast ion <-> fast ion fusion: '//trim(zlbl)
           zlblobj=trim(zlbl)//'_2D'
           call xplasma_mcgrid_putobj( &
                s,id_mcgrid,zlblobj,idpro,jerr, &
                data_1d=bbnvol(1:mimxbzf,ireact), &
                label=label,units='#/sec',normcode=-1)
           call ckerr(jerr,zlblobj)
        end if

     end if

  end do
  !
  !  finish up -- copy stuff for the nb state file
  !
  !  this code executed even if there were no orbits...
  !
1000 continue
  if(nlbout) then
     call lost_orbit_out(jerr)
     ierr=max(ierr,jerr)
  end if
  if((nlbout.or.nltrk_dep0).and.nbeam.gt.0) then
     call beam_info_out(jerr)
     ierr=max(ierr,jerr)
  end if

  !
  tfluxp = tflux
  tbm1 = tbm2                       ! close timestep
  !
  return
  !-----------------------------------------------------------------------
contains

  subroutine ckerr(jerr,zlbl)
    use iso_c_binding, only: fp => c_double

    integer, intent(in) :: jerr
    character(len=*), intent(in) :: zlbl

    if(jerr.ne.0) then
       write(nonlin,*) &
            ' ?outptb2: xplasma_mcgrid_putobj error on item: '// &
            trim(zlblobj)
       call xplasma_error(s,jerr,nonlin)
       ierr=1
    end if

  end subroutine ckerr

  subroutine mklabel
    character(len=12) :: zroot
    character(len=20) :: ztail

    integer :: indx

    indx = index(zlblobj,'_')
    zroot = zlblobj(1:indx-1)
    ztail = zlblobj(indx+1:)

    if(zroot.eq.'SRBN02') then
       indx=index(ztail,'_')
       label = 'neutral source '//ztail(1:indx-1)//' due to '// &
            trim(ztail(indx+1:))

    else if(zroot.eq.'SNBN02') then
       indx=index(ztail,'_')
       label = 'CX rate of thermal '//ztail(1:indx-1)// &
            ' due to '//trim(ztail(indx+1:))//' ions'

    else if(zroot.eq.'SNB0I2') then
       indx=index(ztail,'_')
       label = 'ionization rate of thermal '//ztail(1:indx-1)// &
            ' due to '//trim(ztail(indx+1:))//' ions'

    else if(zroot.eq.'BDENS2') then
       label = 'density of '//trim(ztail)//' ions'
    else if(zroot.eq.'BN0CX2') then
       label = 'density of '//trim(ztail)//' CX fast neutrals'

    else if(zroot.eq.'EBA2PP') then
       label = '<Eperp> of '//trim(ztail)//' ions'
    else if(zroot.eq.'EBA2PL') then
       label = '<Epll> of '//trim(ztail)//' ions'
    else if(zroot.eq.'BE0CX2') then
       label = '<E> of '//trim(ztail)//' CX fast neutrals'

    else if(zroot.eq.'VBTOR2') then
       label = 'toroidal velocity of '//trim(ztail)//' ions'
    else if(zroot.eq.'BV0CX2') then
       label = 'toroidal velocity of '//trim(ztail)// &
            ' CX fast neutrals'

    else if(zroot.eq.'XCURB2') then
       label = 'unshielded current of '//trim(ztail)//' ions'

    else if(zroot.eq.'XCRB2C') then
       label = 'unshielded current of co-injected '// &
            trim(ztail)//' ions'

    else if(zroot.eq.'FBM') then
       label = trim(ztail)//' GC distribution function'

    else if(zroot.eq.'FBM_PTCL') then
       label = trim(ztail)//' distribution function'

    else
       label = ' ??? '
       write(nonlin,*) &
            ' ?outptb2: root not recognized for labeling: ',zroot

    end if

  end subroutine mklabel

  subroutine check(status,varnam)
    integer, intent ( in) :: status
    character(len=*),intent(in) :: varnam ! dimension name input

    if(status /= nf90_noerr) then
       write(nonlin,*) ' ?nbi_wrdeptrk error:'// &
            trim(nf90_strerror(status))//' for var.= '//trim(varnam)
       call exit(status)
    end if
  end subroutine check



end subroutine outptb2

