subroutine nbdep_debug_plot
  !
  ! plot debug data after dep call with extra debug output captured
  !
 
  use nbi_com
  use nbdep_debug
 
  !--------------------------------------------------
 
      implicit none
      integer, parameter :: r8=selected_real_kind(12,100)
  real*8, dimension(:,:), allocatable :: xdatzb,ydatzb,xdatzc,ydatzc
  integer ibrys,i,inum,ibtyp,isel,idcod,iert
  integer, parameter :: m9=9
  integer iself(m9)
  real*8 zself(m9)
  character*10 xuns,yuns,lbsel(m9)
  character*17 kylab(5)
  character*60 mnlab(3)
 
  logical lchar,lxout
 
  !--------------------------------------------------
 
  mnlab(1)=runid
  write(mnlab(2),1001) 'nubeam (nbltes) centerline test ib=', &
       nbdep_set_beam,' ien=',nbdep_set_ien
1001 format(a,i2,a,i1)
  kylab=' '
  xuns='cm'
  kylab(1)='   data items   '
 
  if(.not.nbdep_debug_active) then
     write(nonlin,*) ' %nbdep_debug_plot: nbdep_debug module inactive.'
     return
  endif
 
  if(nzcrtk.eq.0) then
     write(nonlin,*) ' %nbdep_debug_plot: track missed plasma.'
     return
  endif
 
  ibrys=nbdep_nzbdys
  allocate(xdatzb(ibrys,m9),ydatzb(ibrys,m9)); xdatzb=0; ydatzb=0
  allocate(xdatzc(nzcrtk,m9),ydatzc(nzcrtk,m9)); xdatzc=0; ydatzc=0
 
  do i=1,m9
     xdatzc(1:nzcrtk,i)=lz_nbdep(1:nzcrtk)
     xdatzb(1:ibrys,i)=lb_nbdep(1:ibrys)
     iself(i)=i
     zself(i)=i
  enddo
  iself(1)=-1  ! no square scale
 
  do
 
     if(lxout(0)) write(nonlin,1005) ndepmod,xdepmod
 1005 format(/' ndepmod=',i1,'  xdepmod=',1pe11.4/ &
           ' enter plot option #, or "q" to quit:'/ &
           ' 1 -- {r(len),z(len)}'/ &
           ' 2 -- {xi(len),th(len),phi(len)}'/ &
           ' 3 -- {fbeam(len) -- survival fraction'/ &
           ' 4 -- eb_plasma, ti, te vs. len'/ &
           ' 5 -- ne, ion densities, vs. len'/ &
           ' 6 -- mean free path, vs. len'/ &
           ' 7 -- delta(fbeam) total & by atomic processes, vs. len'/ &
           ' 8 -- <sigma*v>{preact,inferred}:electron impact ionization'/ &
           ' 9 -- <sigma*v>{preact,inferred}:ion impact ionization'/ &
           ' 10-- <sigma*v>{preact,inferred}:impurity impact ionization'/ &
           ' 11-- <sigma*v>{preact,inferred}:h-charge-exchange'/ &
           ' 12-- <sigma*v>{preact,inferred}:he-charge-exchange'/ &
           ' 13-- excited states correction factor, vs. len'/ &
           )
     call uread('nbdep_debug_plot option:$')
     if(lchar('q')) exit
 
     isel=idcod(iert)
     if(iert.ne.0) cycle
 
     if(isel.eq.1) then
        ibtyp=1
        inum=2
        lbsel(1)='r'
        lbsel(2)='z'
        yuns='cm'
        ydatzb(1:ibrys,1)=r_nbdep(1:ibrys)
        ydatzb(1:ibrys,2)=y_nbdep(1:ibrys)
        mnlab(3)='(r,z) vs. beam path length'
 
     else if(isel.eq.2) then
        ibtyp=1
        inum=3
        lbsel(1)='xi'
        lbsel(2)='th'
        lbsel(3)='phi'
        yuns=' '
        ydatzb(1:ibrys,1)=xi_nbdep(1:ibrys)
        ydatzb(1:ibrys,2)=th_nbdep(1:ibrys)
        ydatzb(1:ibrys,3)=phi_nbdep(1:ibrys)
        mnlab(3)='(xi=sqrt(norm.tor.flux),theta,phi) vs. beam path length'
 
     else if(isel.eq.3) then
        ibtyp=1
        inum=1
        lbsel(1)='fbeam'
        yuns=' '
        ydatzb(1:ibrys,1)=fbeam_nbdep(1:ibrys)
        mnlab(3)='surviving beam fraction, vs. beam path length'
 
     else if(isel.eq.4) then
        ibtyp=2
        inum=3
        lbsel(1)='eb'
        lbsel(2)='ti'
        lbsel(3)='te'
        ydatzc(1:nzcrtk,1)=ebplas_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,2)=ti_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,3)=te_nbdep(1:nzcrtk)
        yuns='ev'
        mnlab(3)='eb(plasma frame), ti, te vs. beam path length'
 
     else if(isel.eq.5) then
        ibtyp=2
        inum=6
        lbsel(1)='ne'
        lbsel(2)='nh+nd+nt'
        lbsel(3)='nhe+nhe3'
        lbsel(4)='nli'
        lbsel(5)='nimpurity'
        lbsel(6)='nfast'
        ydatzc(1:nzcrtk,1)=ne_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,2)=nhs_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,3)=nhes_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,4)=nlis_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,5)=nimp_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,6)=nfas_nbdep(1:nzcrtk)
        yuns='#/cm**3'
        mnlab(3)='electron and ion densities vs. beam path length'
 
     else if(isel.eq.6) then
        ibtyp=2
        inum=1
        lbsel(1)='lamda(mfp)'
        ydatzc(1:nzcrtk,1)=lmfp_nbdep(1:nzcrtk)
        yuns='cm'
        mnlab(3)='mean free path vs. beam path length'
 
     else if(isel.eq.7) then
        ibtyp=2
        inum=9
        lbsel(1)='dfb(tot)'
        lbsel(2)='dfb(ie)'
        lbsel(3)='dfb(ii)'
        lbsel(4)='dfb(ix)'
        lbsel(5)='dfb(lor.)'
        lbsel(6)='dfb(cxh)'
        lbsel(7)='dfb(cxhe)'
        lbsel(8)='dfb(bbii)'
        lbsel(9)='dfb(bbcx)'
        ydatzc(1:nzcrtk,1)=delf_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,2)=delfie_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,3)=delfii_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,4)=delfix_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,5)=delflo_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,6)=delfcxh_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,7)=delfcxhe_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,8)=delfbbii_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,9)=delfbbcx_nbdep(1:nzcrtk)
        yuns=' '
        mnlab(3)='delta(fb)/dl vs. beam path length'
 
     else if(isel.eq.8) then
        ibtyp=2
        inum=2
        lbsel(1)='<sig*v>ie'
        lbsel(2)='preact ck'
        ydatzc(1:nzcrtk,1)=svie_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,2)=svie2_nbdep(1:nzcrtk)
        yuns='cm**3/sec'
        mnlab(3)='<sigma*v>ie vs. beam path length'
 
     else if(isel.eq.9) then
        ibtyp=2
        inum=2
        lbsel(1)='<sig*v>ii'
        lbsel(2)='preact ck'
        ydatzc(1:nzcrtk,1)=svii_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,2)=svii2_nbdep(1:nzcrtk)
        yuns='cm**3/sec'
        mnlab(3)='<sigma*v>ii vs. beam path length'
 
     else if(isel.eq.10) then
        ibtyp=2
        inum=2
        lbsel(1)='<sig*v>ix'
        lbsel(2)='preact ck'
        ydatzc(1:nzcrtk,1)=svix_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,2)=svix2_nbdep(1:nzcrtk)
        yuns='cm**3/sec'
        mnlab(3)='<sigma*v>impurity ioniz. vs. beam path length'
 
     else if(isel.eq.11) then
        ibtyp=2
        inum=2
        lbsel(1)='<sig*v>cxh'
        lbsel(2)='preact ck'
        ydatzc(1:nzcrtk,1)=svcxh_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,2)=svcxh2_nbdep(1:nzcrtk)
        yuns='cm**3/sec'
        mnlab(3)='<sigma*v>cx(h) vs. beam path length'
 
     else if(isel.eq.12) then
        ibtyp=2
        inum=2
        lbsel(1)='<sg*v>cxhe'
        lbsel(2)='preact ck'
        ydatzc(1:nzcrtk,1)=svcxhe_nbdep(1:nzcrtk)
        ydatzc(1:nzcrtk,2)=svcxhe2_nbdep(1:nzcrtk)
        yuns='cm**3/sec'
        mnlab(3)='<sigma*v>cx(he) vs. beam path length'
 
     else if(isel.eq.13) then
        ibtyp=2
        inum=1
        lbsel(1)='exc.factor'
        ydatzc(1:nzcrtk,1)=excit_nbdep(1:nzcrtk)
        yuns=' '
        mnlab(3)='excitation correction factor vs. beam path length'
 
     else
        call uermsi('invalid plot option',isel,0)
        cycle
     endif
 
     if(ibtyp.eq.1) then
        call r8_grfmm1(xdatzb,ydatzb,ibrys,ibrys,1, &
             zself,inum,iself,inum, &
             xuns,yuns,mnlab,kylab,lbsel)
     else
        call r8_grfmm1(xdatzc,ydatzc,nzcrtk,nzcrtk,1, &
             zself,inum,iself,inum, &
             xuns,yuns,mnlab,kylab,lbsel)
     endif
 
  enddo
 
  deallocate(xdatzb,xdatzc,ydatzb,ydatzc)
 
end subroutine nbdep_debug_plot
! 11sep2004 fgtok -s r8_precision.sub nubeam_test_r8.sub "r8con.csh conversion"
