subroutine nbdrive_summary_2d(ierr)

  use nbdrive_common
  implicit NONE

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

  !------------------------------
  !  this routine reads in 2d outputs from NUBEAM
  !  and then generates some summary output for a 
  !  somewhat arbitrarily chosen zone on the 2d grid.
  !
  !  NUBEAM uses the irregular 2d MCgrid, recently documented
  !  in XPLASMA.
  !
  !  this routine shows a (not necessarily very realistic) example
  !  of using 2d outputs from NUBEAM.
  !------------------------------

  real*8 rmin,rmax,zmin,zmax,zrho_bdy,ztol,zconv,zone_volume
  real*8 zrho(1),zth(1),zphi(1),raxis(1),zaxis(1),zR(1),zZ(1)
  integer izon,idum,inzons,ist_th,ist_ph,inume,inuma
  integer i,ii,j,k,iz,ia,nreact,ityp1,ityp2
  real*8 zvpov,ze,zde,zepll,zeprp,zv,zvpll,zdn

  !  array for scalar function on 2d grid: f(rho,theta)
  real*8, dimension(:), allocatable :: zf2d

  !  array for fast ion distribution function f(E,vpll/v,rho,theta) on 2d grid:
  real*8, dimension(:,:,:), allocatable :: zfbm

  real*8 :: zdens_fbm(nsfast),zeperp_fbm(nsfast),zepll_fbm(nsfast), &
       zvpll_fbm(nsfast)

  integer izmap  ! MCgrid zone index mapping function

  character*20 zflbl  ! fast specie label
  character*20 zthlbl ! thermal specie label
  character*32 zlbl   ! general purpose item label

  !------------------------------
  ! preliminaries:  find rho on axis (rho on bdy -> zdum)
  !    find (R,Z) of mag. axis: (Raxis,Zaxis)
  !    find (Rmin,Rmax) and (Zmin,Zmax) of plasma bdy
  !    find (rho,theta) of
  !       (Raxis,Zaxis) + [(Rmax,Zmax)-(Rmin,Zmin)]/10
  !       (use XPLASMA calls)
  !
  !    Use this (rho,theta) to select a single zone out of the 
  !    2d grid.  All test output will use data from this spatial 
  !    zone only.
  !    

  call eq_rholim(zrho(1),zrho_bdy)  ! rho on axis, rho at bdy

  ! should always be 0 & 1; code may exists that assumes these
  ! limits for the rho coordinate, without making the eq_rholim call.

  !  vectors of length 1 ...
  zth=0; zphi=0

  !  (R,Z) in meters here ... XPLASMA units always MKS

  call eq_rzget(1,zrho,zth,zphi,Raxis,Zaxis,idum,ierr)
  if(ierr.ne.0) go to 9999

  call eq_glimrz(zrho_bdy,Rmin,Rmax,Zmin,Zmax,ierr)
  if(ierr.ne.0) go to 9999
  
  zR(1) = Raxis(1) + (Rmax-Rmin)/10
  zZ(1) = Zaxis(1) + (Zmax-Zmin)/10

  ztol = 1.0d-5  ! inverse map accuracy tolerance (meters)
  call eq_inv(1,zR,zZ,zphi,ztol,zrho,zth,idum,ierr)
  if(ierr.ne.0) go to 9999

  ! zrho(1),zth(1) now give the flux coordinates corresponding to (zR(1),zZ(1))

  write(lun_summary,1001) zR(1)*100.0d0,zZ(1)*100.0d0,zrho(1),zth(1)
1001 format(/ &
          ' NUBEAM 2d output sampled at [R,Z] = [',1pd10.3,',',1pd10.3,'] cm'/&
          '        flux coords: (rho,theta)   = (',1pd10.3,',',1pd10.3,')')

  !  this is the spatial zone to be used for test output:
  !    (integer function izmap shows an example of code for mapping from
  !    (rho,theta) to zone index in a MCgrid.
  !    NOTE: izmap makes an assumption about the theta range coverage,
  !          consistent with the ist_th setting below.

  izon = izmap(id_mcgrid,zrho(1),zth(1))

  write(lun_summary,1002) izon
1002 format(' MCgrid zone index: ',i4/ &
          ' ---> outputs refer to this zone, until noted.')

  !------------------------------
  !  get size of irregular 2d MCgrid
  !  just the size of the part inside the plasma (inzons)
  !  the size including an extrapolated region is returned & ignored (idum)

  call mcgrid_getnumz(id_mcgrid,idum,inzons)

  !  allocate temporary arrays to receive NUBEAM 2d output data

  allocate(zf2d(inzons),zfbm(nznbme,nznbma,inzons))

  !  also set the angle coordinate range controls
  !  this user option applies when data over a MCgrid is fetched.

  ist_ph=0           ! no phi variation (axisymmetry)

  if(nlsym) then
     !  updown symmetric equilibrium
     ist_th=0        ! data will cover poloidal angle range [0,pi]
  else
     ist_th=-1       ! data will cover poloidal angle range [-pi,pi]
     !           (ist_th=0) here would lead to range [0,2pi] being covered.
  endif

  !-------------------------------------------------------------------
  !  in this code, information is read into temporary arrays; only
  !  the information for zone <izon> is saved across calls.  In a real
  !  driver, separate 2d arrays might be allocated for each object.

  !------------------------------
  !  get zone volumes  -- convert to cm**3

  zconv = 1.0d6  ! cm**3 per m**3
  call mcgrid_getobj(id_mcgrid,'ZONE_VOLUME',ist_th,ist_ph, &
       zconv,zf2d,inzons,ierr)
  if(ierr.ne.0) go to 9999

  zone_volume = zf2d(izon)

  write(lun_summary,1005) 'MCgrid zone volume',zone_volume,'cm**3'
1005 format(3x,a,': ',1pd10.3,1x,a)

  !------------------------------
  !  read distribution function
  !  stored in XPLASMA as ptcls/(dsa)/KeV/m**3
  !  convert units to ptcls/(dsa)/ev/cm**3  => 10**-9 conversion factor
  !    (dsa) means "delta-solid-angle" (in velocity space) which works
  !    out to (1/nznbma), for the vpll/v = cos(pitch-angle) grid [-1,1]
  !    evenly divided into NZNBMA pitch zones.
  !
  !    #ptcls in energy/pitch/space zone (i,j,k) is
  !     zfbm(i,j,k)*dE(i)*dVol(k)/NZNBMA
  !
  !  velocity space grids (E,vpll/v) are in the lab-stationary reference frame.

  zdens_fbm=0
  zeperp_fbm=0
  zepll_fbm=0
  zvpll_fbm=0

  do i=1,nsfast
     inume=nznbmea(i)  ! no of energy grid pts expected
     inuma=nznbma      ! no of pitch (vpll/v) grid pts expected

     call nbdrive_mapfast(i,ii,zflbl)  ! note use of label in object name 

     zconv = 1.0d-9    ! fbm units conversion
     call mcgrid_getarr2(id_mcgrid,'FBM_'//trim(zflbl),ist_th,ist_ph, &
          zconv,ZFBM,nznbme,nznbma,inzons, &  ! conversion & array dimensions
          inume,inuma,ierr)                   ! no. of v-space pts expected

     ! for the chosen zone, add up the #ptcls, pll energy, perp energy,
     ! and toroidal momentum; then back out density, <Epll>, <Eperp> and
     ! <vpll>.  These will be compared with outputs directly available
     ! from NUBEAM.

     do j=1,nznbma
        zvpov= -1.0d0 + 2.0d0*(j-0.5d0)/nznbma
        do k=1,inume
           ze=(efbmb(k+1,i)+efbmb(k,i))/2  ! eV (grid set up earlier)
           zde=(efbmb(k+1,i)-efbmb(k,i))
           zepll=zvpov*zvpov*ze
           zeprp=(1.0d0-zvpov*zvpov)*ze
           zv=1.384d6*sqrt(ze/abeams(i))   ! cm/sec
           zvpll=zvpov*zv

           zdn = zfbm(k,j,izon)*zde*zone_volume/nznbma

           zdens_fbm(i)=zdens_fbm(i)+zdn
           zeperp_fbm(i)=zeperp_fbm(i)+zdn*zeprp
           zepll_fbm(i)=zepll_fbm(i)+zdn*zepll
           zvpll_fbm(i)=zvpll_fbm(i)+zdn*zvpll

        enddo
     enddo

     if(zdens_fbm(i).gt.0.0d0) then
        zeperp_fbm(i)=zeperp_fbm(i)/zdens_fbm(i)  ! <Eperp> from fbm
        zepll_fbm(i)=zepll_fbm(i)/zdens_fbm(i)    ! <Epll>  from fbm
        zvpll_fbm(i)=zvpll_fbm(i)/zdens_fbm(i)    ! <vpll>  from fbm
        zdens_fbm(i)=zdens_fbm(i)/zone_volume     ! local density from fbm
     endif

  enddo

  !  the above sums will be checked against data provided by NUBEAM
  !  as 2d scalar functions, below; match will confirm consistency
  !  of fbm use with its normalization inside NUBEAM...

  write(lun_summary,'(1x)')
  !------------------------------
  !  get bdens2 data

  do i=1,nsfast

     call nbdrive_mapfast(i,ii,zflbl)  ! note use of label in object name 

     zconv = 1.0d0                     ! #ptcls/zone returned
     call mcgrid_getobj(id_mcgrid,'BDENS2_'//trim(zflbl), &
          ist_th,ist_ph,zconv,zf2d,inzons,ierr)
     if(ierr.ne.0) go to 9999

     write(lun_summary,1035) zflbl,'n(BDENS2)=',zf2d(izon)/zone_volume, &
          'n(FBM)=',zdens_fbm(i),'ptcls/cm3'

  enddo

  !------------------------------
  !  get (lab frame) <Eperp> data

  do i=1,nsfast

     call nbdrive_mapfast(i,ii,zflbl)  ! note use of label in object name 

     zconv = 1000.0d0                  ! KeV converted to eV
     call mcgrid_getobj(id_mcgrid,'EBA2PP_'//trim(zflbl), &
          ist_th,ist_ph,zconv,zf2d,inzons,ierr)
     if(ierr.ne.0) go to 9999

     write(lun_summary,1035) zflbl,'Eperp(EBA2PP)=',zf2d(izon), &
          'Eperp(FBM)=',zeperp_fbm(i),'eV'

  enddo

  !------------------------------
  !  get (lab frame) <Epll> data

  do i=1,nsfast

     call nbdrive_mapfast(i,ii,zflbl)  ! note use of label in object name 

     zconv = 1000.0d0                  ! KeV converted to eV
     call mcgrid_getobj(id_mcgrid,'EBA2PL_'//trim(zflbl), &
          ist_th,ist_ph,zconv,zf2d,inzons,ierr)
     if(ierr.ne.0) go to 9999

     write(lun_summary,1035) zflbl,'Epll(EBA2PL)=',zf2d(izon), &
          'Epll(FBM)=',zepll_fbm(i),'eV'

  enddo

  !------------------------------
  !  get (lab frame) <v_toroidal> data; compare to vpll from FBM summation...

  do i=1,nsfast

     call nbdrive_mapfast(i,ii,zflbl)  ! note use of label in object name 

     zconv = 100.0d0                   ! m/sec converted to cm/sec
     call mcgrid_getobj(id_mcgrid,'VBTOR2_'//trim(zflbl), &
          ist_th,ist_ph,zconv,zf2d,inzons,ierr)
     if(ierr.ne.0) go to 9999

     !  note that an exacte match is not expected here; v_toroidal
     !  and v_pll are not exactly the same...

     write(lun_summary,1035) zflbl,'v_toroidal(VBTOR2)=',zf2d(izon), &
          'vpll(FBM)=',zvpll_fbm(i),'cm/sec'

  enddo

1035 format(1x,a,1x,a,1pd10.3,1x,a,1pd10.3,1x,a)

  !-----------------------------------------------------------------------
  !  CX fast neutral densities, average energy & toroidal velocity
  !  these are statistics on neutral densities due to CX neutralization 
  !  of partially-slowed-down fast ion species
  !     (neutral density due to initial beam injection is not included).

  write(lun_summary,'(1x)')
  !------------------------------
  do i=1,nsfast

     call nbdrive_mapfast(i,ii,zflbl)  ! note use of label in object name 

     zconv = 1.0d0                     ! #/zone fetched
     call mcgrid_getobj(id_mcgrid,'BN0CX2_'//trim(zflbl), &
          ist_th,ist_ph,zconv,zf2d,inzons,ierr)
     if(ierr.ne.0) go to 9999

     write(lun_summary,1038) zflbl,'cx fast neutral density =', &
          zf2d(izon)/zone_volume,'[neutral atoms]/cm3'

  enddo

  !------------------------------
  do i=1,nsfast

     call nbdrive_mapfast(i,ii,zflbl)  ! note use of label in object name 

     zconv = 1000.0d0                  ! KeV -> eV
     call mcgrid_getobj(id_mcgrid,'BE0CX2_'//trim(zflbl), &
          ist_th,ist_ph,zconv,zf2d,inzons,ierr)
     if(ierr.ne.0) go to 9999

     write(lun_summary,1038) zflbl,'cx fast neutral <energy> =', &
          zf2d(izon),'eV'

  enddo

  !------------------------------
  do i=1,nsfast

     call nbdrive_mapfast(i,ii,zflbl)  ! note use of label in object name 

     zconv = 100.0d0                   ! m/sec -> cm/sec
     call mcgrid_getobj(id_mcgrid,'BV0CX2_'//trim(zflbl), &
          ist_th,ist_ph,zconv,zf2d,inzons,ierr)
     if(ierr.ne.0) go to 9999

     write(lun_summary,1038) zflbl,'cx fast neutral <v_toroidal> =', &
          zf2d(izon),'cm/sec'

  enddo

1038 format(1x,a,1x,a,1pd10.3,1x,a)
  !------------------------------
  ! thermal neutral sources and sinks
  
  do i=1,nsfast

     call nbdrive_mapfast(i,ii,zflbl)  ! note use of label in object name

     do j=1,ngmax

        iz=backz(j) + 0.1d0    ! Z
        ia=aplasm(j) + 0.1d0   ! A
        call namels(iz,ia,zthlbl)  ! returns "H","D","T","He3","He4"
        !  as appropriate for the Z & A values

        write(lun_summary,'(/3x,a)') &
             'interaction of '//trim(zflbl)//' and '// &
             trim(zthlbl)//' neutral atoms:'

        zlbl = 'SRBN02_'//trim(zthlbl)//'0_'//trim(zflbl)

        !  example: SRBN02_H0_D_MCBEAM refers to the hydrogen atom
        !    source due to action of NUBEAM simulation of the D neutral beam

        !  data retrieved is #atoms/sec in each zone!

        !  passing in -id_mcgrid means: set output data to zero without
        !  complaint, if the object does not exist.  This can happen e.g.
        !  with D beams and a He4 target species, because the D neutrals
        !  do not have enough electrons to create He4 neutrals via charge
        !  exchange

        zconv=1.0d0
        call mcgrid_getobj(-id_mcgrid,zlbl, &
             ist_th,ist_ph,zconv,zf2d,inzons,ierr)

        write(lun_summary,1038) zlbl,'neutral source',zf2d(izon)/zone_volume, &
             'atoms/cm3/sec'

        !  thermal neutral sink rate functions (1/sec) are also provided on the
        !  2d grid  (conversion factor = 1.0 retained)

        write(lun_summary,'("   stopping rates for thermal neutrals:")')

        zlbl(1:7)='SNBN02_'   ! sink due to charge exchange
        call mcgrid_getobj(-id_mcgrid,zlbl, &
             ist_th,ist_ph,zconv,zf2d,inzons,ierr)

        write(lun_summary,1038) zlbl,'<nfast*sigma*v>(cx)',zf2d(izon),'1/sec'

        zlbl(1:7)='SNB0I2_'   ! sink due to impact ionization
        call mcgrid_getobj(-id_mcgrid,zlbl, &
             ist_th,ist_ph,zconv,zf2d,inzons,ierr)

        write(lun_summary,1038) zlbl,'<nfast*sigma*v>(ii)',zf2d(izon),'1/sec'
     enddo
  enddo

  !------------------------------
  !  2d fusion rate data -- loop over reactions
  !
  !  library routines (associated with NUBEAM) contain a catalog of 
  !  possible reactions and associated labels used for data identification
  !  purposes.  This catalog spans not just the different types of fusion
  !  reactions (e.g. DD, DT) but also the possible fast ion - thermal target
  !  and fast ion - fast ion reagent combinations.

  write(lun_summary,'(/" fusion rate data on 2d grid:"/)')

  call catfus_num(nreact)  ! no. of possible reactions

  do i=1,nreact

     call catfus_mtype(i,ityp1)  ! available fast ion type, 1st reagent
     ityp2=0                     ! 2nd reagent: thermal target

     call catfus_name(i,ityp1,ityp2,zlbl)
     if(zlbl.eq.'error') then
        !  thermal target should always be allowed
        write(lun_msgs,*) ' ?nbdrive_summary_2d: 2d fusion rate catalog error.'
        write(lun_msgs,*) '  i, ityp1, ityp2 = ',i,ityp1,ityp2
     else

        zlbl = trim(zlbl)//'_2d'

        zconv=1.0d0
        call mcgrid_getobj(-id_mcgrid,zlbl, &
             ist_th,ist_ph,zconv,zf2d,inzons,ierr)

        if(sum(zf2d).gt.0.0d0) then
           write(lun_summary,1038) zlbl,'fusion',zf2d(izon)/zone_volume, &
                'reactions/cm3/sec'
        endif
     endif

     ityp2=1                     ! 2nd reagent: non-thermal beam ion target
     call catfus_name(i,ityp1,ityp2,zlbl)
     if(zlbl.eq.'error') then
        continue  ! not an allowed combination
                  ! (for beam-beam, D->He3 and He3->D are the same)
     else

        zlbl = trim(zlbl)//'_2d'

        zconv=1.0d0
        call mcgrid_getobj(-id_mcgrid,zlbl, &
             ist_th,ist_ph,zconv,zf2d,inzons,ierr)

        if(sum(zf2d).gt.0.0d0) then
           write(lun_summary,1038) zlbl,'fusion',zf2d(izon)/zone_volume, &
                'reactions/cm3/sec'
        endif
     endif
  enddo

  !------------------------------
  
  write(lun_summary,1099) izon
 1099 format(/' ---> end of output for MCgrid zone: ',i4)

9999 continue
  if(allocated(zf2d)) deallocate(zf2d,zfbm)

end subroutine nbdrive_summary_2d
