subroutine nbdrive_fbm_write(filename,ierr)
  use iso_c_binding, only: fp => c_double
  ! write GET_FBM style NetCDF file containing fbm and geometry data
  use nbdrive_common
  use nbi_types
  implicit none
  character(len=*), intent(in) :: filename   ! name of file to write; if blank,
                                          ! generate a default filename...
  integer,  intent(out) :: ierr  ! status code returned (0=OK)

  !----------------------------------------------
  external nbdrive_getrz

  integer :: ii,jj,inum,inume,inumeb,inuma,idum,ist_th,ist_ph,iwarn,istat
  integer :: idtmp,ilun

  type(nbitype_flr) :: zflr
  logical :: nlfbmflr

  real(fp) :: zconv
  character(len=20), dimension(:), allocatable :: zflbls
  character(len=40) :: zname
  !----------------------------------------------
  !  first remove old file -- if there is no data, there will be no file...

  call find_io_unit(ilun)
  open(unit=ilun,file=filename,status='old',iostat=istat)
  if(istat.eq.0) then
     write(lun_msgs,*) ' %nbdrive_fbm_write: removed old copy: '// &
          trim(filename)
     close(unit=ilun,status='delete',iostat=istat)
  end if

  ! get the fast ion labels from NUBEAM -- needed to form names of
  ! xplasma "mcgrid" items --

  call nbfi_numfast(inum)
  allocate(zflbls(inum))

  call nbfi_lbl(zflbls,inum,idum)

  if(.not.allocated(fbm)) then
     call alloc_data(ierr)
     if(ierr.ne.0) return
  end if

  ! set angle range controls

  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.
  end if

  ! get volumes

  zconv = 1.0e6_fp  ! m**3 -> cm**3
  call mcgrid_getobj(id_mcgrid,'ZONE_VOLUME',ist_th,ist_ph, &
       zconv,bmvol,nzones_mc,iwarn)
  if(iwarn.ne.0) then
     write(lun_msgs,*) ' %nbdrive_fbm_write: ZONE_VOLUME: not found.'
     write(lun_msgs,*) '  => fast ion distribution data not available.'
     return
  end if

  ! get 2d density profiles
  !   indx_fi, mapping to NUBEAM fast species index, was set up in prior
  !   call to nbdrive_load(...)

  iwarn = 0
  zconv = 1.0_fp  ! #ptcls/zone returned...
  do ii=1,nsfast

     jj=indx_fi(ii)  ! NUBEAM index
     zname = 'BDENS2_'//trim(zflbls(jj))

     call mcgrid_getobj(id_mcgrid,zname,ist_th,ist_ph, &
          zconv,bdens2(:,ii),nzones_mc,ierr)
     if(ierr.ne.0) then
        write(lun_msgs,*) ' %nbdrive_fbm_write: failed to fetch: '//trim(zname)
        write(lun_msgs,*) '  => fast ion distribution data not available (OK).'
        iwarn = iwarn + 1
        ierr = 0
     end if

  end do
  if(iwarn.ne.0) return

  ! get <Eperp> and <Epll>
  zconv = 1000.0_fp  ! KeV -> eV
  do ii=1,nsfast

     jj=indx_fi(ii)  ! NUBEAM index
     zname = 'EBA2PP_'//trim(zflbls(jj))

     call mcgrid_getobj(id_mcgrid,zname,ist_th,ist_ph, &
          zconv,eba2pp(:,ii),nzones_mc,ierr)
     if(ierr.ne.0) then
        write(lun_msgs,*) ' ?nbdrive_fbm_write: failed to fetch: '//trim(zname)
        exit
     end if

     zname = 'EBA2PL_'//trim(zflbls(jj))

     call mcgrid_getobj(id_mcgrid,zname,ist_th,ist_ph, &
          zconv,eba2pl(:,ii),nzones_mc,ierr)
     if(ierr.ne.0) then
        write(lun_msgs,*) ' ?nbdrive_fbm_write: failed to fetch: '//trim(zname)
        exit
     end if

  end do
  if(ierr.ne.0) return

  zname = 'FBM_PTCL_'//trim(zflbls(1))
  call eq_gfitem(zname,idtmp)
  if(idtmp.gt.0) then
     nlfbmflr = .FALSE.  ! if FBM_PTCL exists, we read the GC dist. fcn
  else
     call nbi_get_flr(zflr)
     nlfbmflr = zflr%nlfbmflr
  end if

  ! get distribution function data
  !   indx_fi, mapping to NUBEAM fast species index, was set up in prior
  !   call to nbdrive_load(...)

  zconv = 1.0e-9_fp  ! #/m^3/KeV/dSA -> #/cm^3/eV/dSA
  do ii=1,nsfast

     jj=indx_fi(ii)  ! NUBEAM index
     zname = 'FBM_'//trim(zflbls(jj))

     inume = nznbmea(ii)  ! #energy pts expected
     inuma = nznbma       ! #vpll/v pts expected

     ! nznbmea(ii) < nznbme possible; both the array sizes and the expected
     ! data sizes need to be specified in the call:

     call mcgrid_getarr2(id_mcgrid,zname,ist_th,ist_ph, &
          zconv,fbm(:,:,:,ii),nznbme,nznbma,nzones_mc, &
          inume,inuma,ierr)

     if(ierr.ne.0) then
        write(lun_msgs,*) ' ?nbdrive_fbm_write: failed to fetch: '//trim(zname)
        exit
     end if

  end do
  if(ierr.ne.0) return

  !------------------------------------------------
  !  gfbm_set calls... transfer data, to prep for NetCDF file write...

  call gfbm_set_species(nsfast,xzbeams,abeams)

  call gfbm_set_dims(nsfast,nznbmea,nznbma,nzones_fb,nzones_mc)

  call gfbm_set_udsym(nlsym)

!MG temprorarily disabeled, April 15,2011
!MG changes in: nbdrive_fbm_write, get_fbm_fnc, fbm_wnetcdf_fp
!  call gfbm_set_fgc(.NOT.nlfbmflr)
  call gfbm_set_fgc(.true.) !MG only fbm at gc will be written to netcdf file

  do ii=1,nsfast
     inumeb = nznbmea(ii)+1
     call gfbm_set_egridb1(ii,inumeb,efbmb(1:inumeb,ii))
  end do

  call gfbm_set_types(nsfast,nsbfp,nlfprod)

  call gfbm_set_nthzns(nzones_fb,nthzns)

  call gfbm_set_bmvol(nzones_mc,bmvol)

  call gfbm_set_bdens2(nzones_mc,nsfast,bdens2)

  call gfbm_set_eba2(nzones_mc,nsfast,eba2pp,eba2pl)

  call gfbm_set_fbm(nznbme,nznbma,nzones_mc,nsfast,fbm)

  ! transfer deposition data, if available...

  call gfbm_set_track_data

  ! write NetCDF file:

  call fbm_wnetcdf_fp(nbdrive_getRZ,runid,filename,ierr)

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

end subroutine nbdrive_fbm_write

subroutine nbdrive_getrz(xi,theta,rr,zz)
  use iso_c_binding, only: fp => c_double

  ! (rho,theta) -> (R,Z) map -- real(fp), return values in cm
  ! will be declared external for fbm_write call

  use nbdrive_common
  implicit none

  real(fp), intent(in) :: xi,theta
  real(fp), intent(out) :: rr,zz

  !--------------------------------
  integer :: iregion,ierr
  real(fp) :: rho(1),th(1),phi(1),rrr(1),zzz(1)
  !--------------------------------

  rho(1) = xi
  th(1) = theta
  phi(1) = 0.0_fp

  call eq_rzget(1, rho,th,phi, rrr,zzz, iregion, ierr)
  if(ierr.ne.0) then
     write(lun_msgs,*) ' ?? nbdrive_getrz: eq_rzget returned ierr= ',ierr
     write(lun_msgs,*) '    (xi,theta) input: ',xi,theta
     call bad_exit
  end if

  rr = 1.0e2_fp*rrr(1)
  zz = 1.0e2_fp*zzz(1)

end subroutine nbdrive_getrz
