subroutine nbdrive_summary(ierr)

  use nbi_types
  use nbdrive_common
  implicit NONE

  !------------------------------
  !  generate summary output from nubeam_driver (nbdrive)...
  !  this is hand coded to demonstrate select NUBEAM output capability.
  !
  !  the subroutine nbdrive_details is written (mostly) by a Python
  !  code generator and produces a comprehensive survey of NUBEAM output
  !------------------------------

  integer,intent(out) :: ierr    ! exit status code (0=OK)

  !------------------------------
  character*20 zflbl  ! fast specie label

  !  to output rho = [0,0.25] and = [0,1] volume and area averages:
  real*8 rho_int(2),vol_int(2),ar_int(2),f_int(2)

  !  for explicit averaging of parameter profiles
  integer, parameter :: inz=20, inzb=21   ! inz should be divisible by 4
  real*8 zxizc(inz),zxizb(inzb),zdvol(inz),zvol(inzb),zf(inz)

  !  miscellaneous local variables...
  integer i,ii,ij,idata
  real*8 zsum

  !------------------------------
  !  output data types returnable from NUBEAM...

  type (nbotype_fusion) :: zfusion
  type (nbotype_density) :: zdensity
  type (nbotype_deposition) :: zdeposition
  type (nbotype_n0_fast) :: zn0_fast
  type (nbotype_trap_fraction) :: ztrap_fraction
  type (nbotype_erngfi_output) :: zerngfi_output
  type (nbotype_average_energies) :: zaverage_energies
  type (nbotype_excited_states) :: zexcited_states
  type (nbotype_mc_statistics) :: zmc_statistics
  type (nbotype_rotation) :: zrotation
  type (nbotype_compression) :: zcompression
  type (nbotype_powers) :: zpowers
  type (nbotype_currents) :: zcurrents
  type (nbotype_neutral_sources) :: zneutral_sources
  type (nbotype_recapture) :: zrecapture
  type (nbotype_neutral_sinks) :: zneutral_sinks
  type (nbotype_sources) :: zsources
  type (nbotype_fokker_planck) :: zfokker_planck
  type (nbotype_radial_current) :: zradial_current
  type (nbotype_torques) :: ztorques
  type (nbotype_hi_z_beams) :: zhi_z_beams
  type (nbotype_power_balance) :: zpower_balance
  type (nbotype_momentum_balance) :: zmomentum_balance
  type (nbotype_torque_work) :: ztorque_work
  type (nbotype_particle_balance) :: zparticle_balance
  type (nbotype_nclass) :: znclass
  !------------------------------

  ierr = 0

  !------------------------------
  !  set up some local quantities to aid in profile summary output
 
  rho_int(1)=0.25d0
  rho_int(2)=1.00d0

  call eq_volume(2,rho_int,0,vol_int,ierr)
  if(ierr.ne.0) return
  ! vol_int(1) contains volume of plasma in region (rho) [0,0.25d0]
  ! vol_int(2) contains volume of plasma in region (rho) [0,1]
  vol_int = vol_int*1.0d6  ! **** in cm**3

  call eq_area(2,rho_int,0,ar_int,ierr)
  if(ierr.ne.0) return
  ! ar_int(1) contains x-sect. area of plasma in region (rho) [0,0.25d0]
  ! ar_int(2) contains x-sect. area of plasma in region (rho) [0,1]
  ar_int = ar_int *1.0d4   ! in cm**2

  do i=1,inzb
     zxizb(i) = (i-1)*1.0d0/(inz)
  enddo
  call eq_volume(inzb,zxizb,0,zvol,ierr)
  if(ierr.ne.0) return

  zvol=zvol*1.0d6  ! **** volumes now in cm**3

  do i=1,inz
     zxizc(i)=(zxizb(i)+zxizb(i+1))/2
     zdvol(i)=(zvol(i+1)-zvol(i))
  enddo

  !------------------------------
  !  some general points: 
  !    (a) nbdrive species labels were generated in nbdrive_nubeam_checknaml
  !
  !    (b) fast species indexing is derived from the neutral beam species
  !        plus (namelist) selected fusion product species.  In general,
  !        the driver code's indexing might not be the same as NUBEAM's
  !        indexing and so: a mapping routine should be used.
  !
  !    (c) all structures for NUBEAM scalar and 1d profile (flux surface
  !        averages) are defined by the contents of nubeam/nbspec.dat ... 
  !        1d profiles are distinguished from scalars by the presence of
  !        the array index "mj".  XPLASMA ids are stored as integer scalar
  !        or array members of a compound structure; the indexing is as 
  !        defined by nbspec.dat with "mj" removed.  For example:
  !
  !          "DA  bdeps(mj,mibs)  ! n/sec (integrated profile) -- deposition"
  !
  !        causes an array of XPLASMA ids io_bdeps(mibs) to be defined as
  !        a member of the output data structure corresponding to the block
  !        in which this record is found ("deposition").
  !
  !        Array indices are defined near the top of nbspec.dat; "mibs"
  !        refers to fast ion species; i.e. one profile (bdeps) is defined
  !        for each fast ion specie in the simulation.
  !
  !        There are a limited number of spatially 2d outputs available--
  !        see subroutine nbdrive_summary_2d.f90 ...
  !
  !    (d) for some control settings it may be normal for certain profiles
  !        to remain undefined.  In this case, the XPLASMA id is zero.  User
  !        code will need to check for this to avoid producing warning
  !        messages from XPLASMA interpolation routines receiving 0 as 
  !        function id.
  !
  !    (e) many profiles are in integrated form-- i.e. total number of
  !        ptcls (energy, power, torque, whatever) vs. rho.  Conversion
  !        to a density requires differencing between two rho values and
  !        dividing out the corresponding volume-- or cross-sectional
  !        area, in the case of toroidal driven current profiles.  The
  !        form and units of output profiles are described in 
  !        nubeam/nbspec.dat.
  !
  !    (f) some profiles (slowing down times, average energies) are 
  !        parameters, not in integrated form.  See nubeam/nbspec.dat.
  !
  !    (g) the set of profiles provided basically constitute the union
  !        of physicists' requests over a 25 year period.  By no means
  !        would a transport code need all of them, but it will require
  !        at least the following:
  !           -- fast ion densities (affects Zeff/quasi-neutrality)
  !           -- fast ion energy content (affects MHD equilibrium)
  !           -- electron and thermal ion source/sinks
  !           -- low energy "thermal" neutral atom sources and sink rates
  !              (a thermal neutral gas transport model is needed to
  !              complete the computation of ptcl/momentum/power transport 
  !              effects associated with the low energy neutral atoms 
  !              created from neutral beams due to charge exchange).
  !           -- momentum sources (torques)
  !           -- power sources.
  !
  !------------------------------
  !  fast ion ptcl and energy density
  !  2d profile information also exists.

  call nbo_get_density(zdensity)

  idata = 0
  do i=1,nsfast
     write(lun_summary,1001) trim(label_fastions(i))
1001 format(/1x,a,':')
     call nbdrive_mapfast(i,ii,zflbl)
     call denso('density','n/cm3',zdensity%io_bdenss(ii))
     call denso('parallel energy density','Jles/cm3',zdensity%io_udenspl(ii))
     call denso('perp energy density','Jles/cm3',zdensity%io_udenspl(ii))
     if( zdensity%io_bdenss(ii) .gt. 0 ) then
        idata = idata + 1
     endif
  enddo

  if(idata.eq.0) then
     write(lun_summary,*) ' '
     write(lun_summary,*) &
          ' --> CAUTION: no fast ions found in nbdrive_summary.f90... '
     write(lun_summary,*) ' '
  endif

  !------------------------------
  !  fusion rates
  !  these are totals; 1d flux surface averages and 2d emission profile 
  !  information exists as well.

  call nbo_get_fusion(zfusion)

  !  scalar info:
  write(lun_summary,1005) zfusion%btneut,zfusion%bbntot

1005 format(/1x,'beam-target neutrons (total): ',1pd10.3/ &
             1x,'  beam-beam neutrons (total): ',1pd10.3)

  !------------------------------
  !  deposition

  call nbo_get_deposition(zdeposition)

  do i=1,nsfast
     write(lun_summary,1001) trim(label_fastions(i))
     call nbdrive_mapfast(i,ii,zflbl)
     call denso('deposition source','n/cm3/sec',zdeposition%io_bdeps(ii))
     call denso('deposition source','n/cm3/sec',zdeposition%io_sdep_oa(ii))
  enddo

  !------------------------------
  !  trapping fraction info

  call nbo_get_trap_fraction(ztrap_fraction)
  
  !  write scalar data -- based on whole fast ion population...

  do i=1,nsfast
     call nbdrive_mapfast(i,ii,zflbl)
     write(lun_summary,1009) trim(label_fastions(i)),ztrap_fraction%tbtrap(ii)
  enddo
1009 format(/1x,a,' trapping fraction: ',1pd10.3)

  !------------------------------
  !  average energies

  call nbo_get_average_energies(zaverage_energies)

  do i=1,nsfast
     write(lun_summary,1001) trim(label_fastions(i))
     call nbdrive_mapfast(i,ii,zflbl)
     call paramo('ion parallel energy','KeV', &
          zaverage_energies%io_ebapl2s(ii))
     call paramo('ion perpendicular energy','KeV', &
          zaverage_energies%io_ebapp2s(ii))
  enddo

  !------------------------------
  !  rotation

  call nbo_get_rotation(zrotation)

  do i=1,nsfast
     write(lun_summary,1001) trim(label_fastions(i))
     call nbdrive_mapfast(i,ii,zflbl)
     call paramo('lab frame angular velocity','1/sec',zrotation%io_omegabs(ii))
  enddo

  !------------------------------
  !  powers -- profiles are available summed over species and separate
  !            for each species; sums are output here.

  call nbo_get_powers(zpowers)

  write(lun_summary,1015)
1015 format(/' from all fast ions:')

  call denso('electron heating (beams)','W/cm3',zpowers%io_pbe)
  call denso('ion heating (beams)','W/cm3',zpowers%io_pbi)

  call denso('beam thermalization power','W/cm3',zpowers%io_pbth)

  call denso('electron heating (fusion products)','W/cm3',zpowers%io_pfe)
  call denso('ion heating (fusion products)','W/cm3',zpowers%io_pfi)

  !------------------------------
  !  driven currents

  call nbo_get_currents(zcurrents)

  write(lun_summary,1015)
  call curdenso('beam driven current','Amps/cm2',zcurrents%io_curjbeam)
  call curdenso('fusion product driven current','Amps/cm2', &
       zcurrents%io_curjfusi)

  !------------------------------
  !  thermal neutral sources (charge exchange)
  !  these are flux surface averages; 2d profile information exists as well.
  !     ** profiles are separated by target plasma species **

  call nbo_get_neutral_sources(zneutral_sources)

  write(lun_summary,1015)
  do i=1,ngmax
     ij=index(label_thions(i),' ')-1
     call denso( &
          label_thions(i)(1:ij)//' halo neutral source','n/cm3/sec', &
          zneutral_sources%io_sorbn0(i))
  enddo

  !------------------------------
  !  neutral sinks (charge exchange and impact ionization)
  !  these are flux surface averages; 2d profile information exists as well.
  !     ** profiles are separated by target plasma species **

  !  note these are rate profiles -- 1/sec -- ready for use by
  !  neutral gas models.

  call nbo_get_neutral_sinks(zneutral_sinks)

  write(lun_summary,1015)
  do i=1,ngmax
     ij=index(label_thions(i),' ')-1
     call paramo( &
          ' beam charge exch. sink rate for '//label_thions(i)(1:ij), &
          '1/sec', zneutral_sinks%io_sinbn0(i))
     call paramo( &
          ' beam impact ioniz sink rate for '//label_thions(i)(1:ij), &
          '1/sec', zneutral_sinks%io_sinb0i(i))
  enddo

  !------------------------------
  !  electron and ion sources

  call nbo_get_sources(zsources)

  write(lun_summary,1015)
  call denso('electron source','n/cm3/sec',zsources%io_sorbe)

  do i=1,ngmax
     ij=index(label_thions(i),' ')-1
     call denso(label_thions(i)(1:ij)//' thermal ion source','n/cm3/sec', &
          zsources%io_sorbh(i))
  enddo

  !------------------------------
  !  collisional timescales

  call nbo_get_fokker_planck(zfokker_planck)

  do i=1,nsfast
     write(lun_summary,1001) trim(label_fastions(i))
     call nbdrive_mapfast(i,ii,zflbl)
     call paramo('90 degree p.a. scatter time','sec', &
          zfokker_planck%io_taubpa0(ii))
     call paramo('slowing down time','sec', &
          zfokker_planck%io_taubsl0(ii))
  enddo
  !------------------------------
  !  torques; only the sum is fetched; separate information for
  !  different momentum transfer processes are also available.

  call nbo_get_torques(ztorques)

  write(lun_summary,1015)
  call denso('net fast ion torque to plasma','Nt-m/cm3',ztorques%io_tqbsum)

  !------------------------------
  !  power balance information

  write(lun_summary,1019)
1019 format(' ... rough power balance ...')

  call nbo_get_power_balance(zpower_balance)
  call nbo_get_compression(zcompression)
  call nbo_get_torque_work(ztorque_work)

  do i=1,nsfast
     write(lun_summary,1001) trim(label_fastions(i))
     call nbdrive_mapfast(i,ii,zflbl)
     write(lun_summary,1088) '+injected power (W): ',zpower_balance%pinjs(ii)
     zsum = zpower_balance%pinjs(ii)
     write(lun_summary,1088) '+OH -> beam ions:    ',zpower_balance%bpohs(ii)
     zsum = zsum + zpower_balance%bpohs(ii)

     zf(1)=0
     if(zcompression%io_pbcprss(ii).gt.0) then
        call eq_rgetf(1,rho_int(2:2),zcompression%io_pbcprss(ii),0,zf(1:1), &
             ierr)
     endif
     zsum = zsum + zf(1)
     write(lun_summary,1088) '+compression input:  ',zf(1)

     !  still have zpowers -- get total electron and ion heating from
     !  this fast specie.  The power profiles in XPLASMA are already 
     !  volume integrals of the power densities, as mentioned before.

     zf(1:2) = 0.0d0

     if(zpowers%io_pbes(ii).gt.0) then
        call eq_rgetf(1,rho_int(2:2),zpowers%io_pbes(ii),0,zf(1:1),ierr)
     endif

     if(zpowers%io_pbis(ii).gt.0) then
        call eq_rgetf(1,rho_int(2:2),zpowers%io_pbis(ii),0,zf(2:2),ierr)
     endif

     zsum = zsum - zf(1) - zf(2)

     write(lun_summary,1088) '-electron heating:   ',zf(1)
     write(lun_summary,1088) '-ion heating:        ',zf(2)

     write(lun_summary,1088) '-"internal" cx loss: ',zpower_balance%bpcxins(ii)
     zsum = zsum - zpower_balance%bpcxins(ii)
     write(lun_summary,1088) '-"external" cx loss: ',zpower_balance%bpcxots(ii)
     zsum = zsum - zpower_balance%bpcxots(ii)
     write(lun_summary,1088) '-shine-through:      ',zpower_balance%bpshins(ii)
     zsum = zsum - zpower_balance%bpshins(ii)
     write(lun_summary,1088) '-bad orbit loss:     ',zpower_balance%bpls(ii)
     zsum = zsum - zpower_balance%bpls(ii)
     write(lun_summary,1088) '-thermalization:     ',zpower_balance%bpths(ii)
     zsum = zsum - zpower_balance%bpths(ii)
     write(lun_summary,1088) '-ripple loss:        ',zpower_balance%bprpls(ii)
     zsum = zsum - zpower_balance%bprpls(ii)
     write(lun_summary,1088) '-fishbone loss:      ',zpower_balance%bpfshs(ii)
     zsum = zsum - zpower_balance%bpfshs(ii)
!MG should use io_pbaepms
!MG     write(lun_summary,1088) '-AEP(MHD)activity    ',sum(zpower_balance%pbaepms(:,:,ii))
!MG      zsum = zsum - sum(zpower_balance%pbaepms(:,:,ii))

     zf(1)=ztorque_work%bpjxbs(ii)
     zf(2)=ztorque_work%bptcols(ii) + ztorque_work%bptdfbs(ii) + &
          ztorque_work%bpthros(ii) + ztorque_work%bpthsfs(ii)

     write(lun_summary,1088) '-torque work (JxB):  ',zf(1)
     write(lun_summary,1088) '-torque work (other):',zf(2)

     zsum = zsum - zf(1) - zf(2)

     write(lun_summary,1088) '-d/dt(f.i. energy):  ',zpower_balance%bpsts(ii)
     zsum = zsum - zpower_balance%bpsts(ii)
     write(lun_summary,1088) '   ->residual:       ',zsum
  enddo

1088 format(4x,a,1x,1pd10.3)

  !---------------------------------------------------------------------
  !  a separate routine is provided to show example of access to 2d
  !  output

  if(idata.gt.0) then
     call nbdrive_summary_2d(ierr)
  else
     write(lun_summary,*) ' '
     write(lun_summary,*) ' %nbdrive_summary_2d call skipped: no fast ions...'
     write(lun_summary,*) ' '
  endif

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

  contains

    subroutine denso(zlbl,zuns,id)
      character*(*), intent(in) :: zlbl   ! item label
      character*(*), intent(in) :: zuns   ! physical units label
      integer, intent(in) :: id           ! XPLASMA id of profile

      if(id.eq.0) then
         write(lun_summary,2000) trim(zlbl),trim(zuns)
2000     format(3x,a,1x,a,': no NUBEAM profile available.')
         return
      endif

      call eq_rgetf(2,rho_int,id,0,f_int,ierr)
      if(ierr.ne.0) then
         write(lun_msgs,*) ' ?nbdrive_summary -- data fetch error: ', &
              trim(zlbl)
      else
         write(lun_summary,2001) &
              ' rho <= 0.25 average  ',trim(zlbl),f_int(1)/vol_int(1), &
              trim(zuns)
         write(lun_summary,2001) &
              ' whole plasma average ',trim(zlbl),f_int(2)/vol_int(2), &
              trim(zuns)
      endif

2001  format(2x,a,a,': ',1pd10.3,1x,a)

    end subroutine denso

    subroutine curdenso(zlbl,zuns,id)
      character*(*), intent(in) :: zlbl   ! item label
      character*(*), intent(in) :: zuns   ! physical units label
      integer, intent(in) :: id           ! XPLASMA id of profile

      if(id.eq.0) then
         write(lun_summary,2000) trim(zlbl),trim(zuns)
2000     format(3x,a,1x,a,': no NUBEAM profile available.')
         return
      endif

      call eq_rgetf(2,rho_int,id,0,f_int,ierr)
      if(ierr.ne.0) then
         write(lun_msgs,*) ' ?nbdrive_summary -- data fetch error: ', &
              trim(zlbl)
      else
         write(lun_summary,2001) &
              ' rho <= 0.25 average  ',trim(zlbl),f_int(1)/ar_int(1), &
              trim(zuns)
         write(lun_summary,2001) &
              ' whole plasma average ',trim(zlbl),f_int(2)/ar_int(2), &
              trim(zuns)
      endif

      write(lun_summary,2001) 'Total ',trim(zlbl),f_int(2),'Amps.'

2001  format(2x,a,a,': ',1pd10.3,1x,a)

    end subroutine curdenso

    subroutine paramo(zlbl,zuns,id)
      character*(*), intent(in) :: zlbl   ! item label
      character*(*), intent(in) :: zuns   ! physical units label
      integer, intent(in) :: id           ! XPLASMA id of profile

      ! output some averages from a fast ion parameter profile.

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

      integer inz25,inz25p1,i
      real*8 zsum,zans25,zans

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

      if(id.eq.0) then
         write(lun_summary,2000) trim(zlbl),trim(zuns)
2000     format(3x,a,1x,a,': no NUBEAM profile available.')
         return
      endif

      call eq_rgetf(inz,zxizc,id,0,zf,ierr)
      if(ierr.ne.0) then
         write(lun_msgs,*) ' ?nbdrive_summary -- data fetch error: ', &
              trim(zlbl)
      else
         inz25 = inz/4
         inz25p1 = inz25+1
         zsum=0
         do i=1,inz
            zsum=zsum+zdvol(i)*zf(i)
            if(i.eq.inz25) then
               zans25=zsum/zvol(inz25p1)  ! vol avg to rho=0.25
            endif
         enddo
         zans = zsum/zvol(inzb)           ! whole plasma vol avg
         write(lun_summary,2001) &
              ' rho <= 0.25 average  ',trim(zlbl),zans25, &
              trim(zuns)
         write(lun_summary,2001) &
              ' whole plasma average ',trim(zlbl),zans, &
              trim(zuns)
      endif
2001  format(2x,a,a,': ',1pd10.3,1x,a)

    end subroutine paramo

end subroutine nbdrive_summary
