subroutine fbm_wnetcdf_fp(rymap,tr_runid,zfile,ierr)
  use iso_c_binding, only: fp => c_double
  use get_fbm_module                ! get_fbm info here
  use ezcdf
  
  implicit none
  
  
  external :: rymap  ! passed subroutine, see below...
  
  character(len=*), intent(in) :: tr_runid ! TRANSP runid
  character(len=*), intent(in) :: zfile ! name of file to write
  integer, intent(out) :: ierr      ! completion code (0=OK)
  !
  !  call rymap(xi,theta,R,Y): flux coordinate -> (R,Y) map
  !       (R,Y) must be returned in cm!!  All rymap arguments: real(fp) ::
  !
  !-----------------------------------------
  !  local: 2d grid info
  integer :: icdf,ins_save,iz,ii,ix,ict,ictx,itot,inum
  integer :: ith1,ith2,nth,isym,igc,ir,is,idepb,iptcl
  integer, dimension(:), allocatable :: index2d
  real(fp), dimension(:), allocatable :: x2d,th2d,r2d,z2d
  real(fp), dimension(:), allocatable :: wx2d,wth2d,wr2d,wz2d
  !  local: flux surface info
  integer :: nxsurf,nthsurf
  real(fp), dimension(:), allocatable :: xsurf,thsurf
  real(fp), dimension(:,:), allocatable :: rsurf,zsurf
  real(fp) :: zx,zth,zth1,zth2
  !  local: species info
  logical :: present(nsfast)
  integer :: imap(nsfast)
  integer :: ilfprod(nsfast),ilbeam(nsfast)
  character(len=10) :: slabel(nsfast),skeep(nsfast),slname
  real(fp) :: znum(nsfast),znumck,ztest
  character(len=3) :: isotop(4,2)           ! (A,Z)
  !
  integer :: imaxe,imaxeb,imaxk,imaxkb,ine,ina,idm_fbm(3)
  real(fp), dimension(:,:), allocatable :: zksid,zksidb, &
       zefbm,zefbmb
  integer, dimension(:,:), allocatable :: ilfusion_sce
  !
  real(fp) :: zta,zdta
  real(fp), dimension(:), allocatable :: zmvol
  real(fp), dimension(:,:,:), allocatable :: zfbm
  !
  real(fp) :: zarg1,zarg2
  !
  character(len=20) :: hint = "dN = f*bmvol*de*da/2"
  !
  character(len=10) :: runid
  integer :: ie,ia,ith
  integer :: ns_use !# of species used for beam deposition
  integer :: inz_track !# of zones for a deposition of a given species
  integer :: inz_p, inz_c
  real(fp) :: fbm_fnc
  real(fp), dimension(:,:,:), allocatable :: wrk_bd2
  real(fp), dimension(:), allocatable :: zindr,zindz,wk
  integer :: i_mult
  
  !
  data isotop/'H','D','T','?', '?','?','He3','He4'/
  
  !
  !-----------------------------------------
  !
  !  write NetCDF file containing TRANSP NBI information
  !
  !  scan species list; check that dist fcn is non-zero and density
  !  check is OK; set up short species label
  !
  runid=tr_runid
  !
  zta=ta
  zdta=dta
  !
  imaxk=nznbma
  imaxkb=imaxk+1
  imaxe=maxval(nznbmea(1:nsfast))
  imaxeb=imaxe+1
  allocate(zksid(imaxk,nsfast),zksidb(imaxkb,nsfast))
  allocate(zefbm(imaxe,nsfast),zefbmb(imaxeb,nsfast))
  !
  ins_save = ns
  itot = 0
  imap = 0
  idepb = 0
  do ns = 1,nsfast
     iz=xzbeams(ns) + 0.5_fp
     ia=abeams(ns) + 0.5_fp
     slabel(ns)='_'//trim(isotop(ia,iz))//'_'//fitype(ktype(ns))
     write(6,*) ' specie: ',slabel(ns)(2:)
     
     call get_fbm_ncheck0(zarg1,zarg2)
     write(6,'(1x,a,1pe11.4)') 'Total #of fast ions (fbm):  ', &
          zarg1
     write(6,'(1x,a,1pe11.4)') 'Total #of fast ions (bdens2):  ', &
          zarg2
     write(6,*) ' (these numbers should match roughly).'
     
     znum(ns)=zarg1
     znumck=zarg2
     if(nldep_beam(ns)) idepb = idepb + 1
     if((.not.nldep_beam(ns)).AND.(.not.nldep_fusion(ns)).AND. &
          (max(znum(ns),znumck).le.0.0_fp)) then
        present(ns)=.FALSE.
     else
        ztest=abs(znum(ns)-znumck)/max(0.1_fp,max(znum(ns),znumck))
        if(ztest.gt.0.1_fp) then
           write(6,*) ' ?density normalization check failed.'
        end if
        
        present(ns)=.TRUE.
        itot=itot+1
        skeep(itot)=slabel(ns)
        imap(itot)=ns
        
     end if
     !
     zksid(1:nznbma,ns)=xxksid(1:nznbma)
     zksidb(1,ns)=-1; zksidb(nznbma+1,ns)=1
     do ii=2,nznbma
        zksidb(ii,ns)=(zksid(ii-1,ns)+zksid(ii,ns))/2.0_fp
     end do
     !
     zefbmb(1:nznbmea(ns)+1,ns)=efbmb(1:nznbmea(ns)+1,ns)
     do ii=1,nznbmea(ns)
        zefbm(ii,ns)=(zefbmb(ii,ns)+zefbmb(ii+1,ns))/2.0_fp
     end do
     !
     ilbeam(ns)=0
     ilfprod(ns)=0
     if(ns.le.nsbeam) then
        if(nlfprod(ns)) then
           ilfprod(ns)=1
        else
           ilbeam(ns)=1
        end if
     else
        continue ! RF minority? -- not beam nor fusion product
     end if
     
  end do
  
  if(itot.eq.0) then
     write(6,*) ' ?fbm_wnetcdf: no valid distribution fcn data.'
     ierr=99
     return
  end if
  
  call cdf_open(icdf,trim(zfile),'w',ierr)
  if(ierr.ne.0) then
     write(6,*) ' ?fbm_wnetcdf: open failure: ',trim(zfile)
     return
  end if
  
  nxsurf=2*nznbmrf+1
  if(nlsym) then
     nthsurf=101
     isym=1
  else
     nthsurf=201
     isym=0
  end if
  
  if(nlfbm_gc) then
     igc=1
  else
     igc=0
  end if
  if(nlfbm_ptcl) then
     !MG temprorarily disabeled, April 15,2011
     !MG changes in: nbdrive_fbm_write, get_fbm_fnc, fbm_wnetcdf_fp
     iptcl=0
     !         iptcl=1
  else
     iptcl=0
  end if
  
  allocate(xsurf(nxsurf),thsurf(nthsurf))
  allocate(rsurf(nthsurf,nxsurf),zsurf(nthsurf,nxsurf))
  
  inum=nthzsm(lcentr+nznbmrf-1)
  idm_fbm(3)=inum                   ! no. of spatial zones in fbm
  allocate(index2d(inum))
  allocate(x2d(inum),th2d(inum),r2d(inum),z2d(inum),zmvol(inum))
  allocate(wx2d(inum*4),wth2d(inum*4),wr2d(inum*4),wz2d(inum*4))
  allocate(wrk_bd2(inum*4,inum*4,nfast))
  allocate(zindr(inum*4),zindz(inum*4),wk(inum*4))
  zindr=0.0_fp
  zindz=0.0_fp
  wrk_bd2=10.0_fp
  wk=0.0_fp
  wx2d=0.0_fp
  wth2d=0.0_fp
  wr2d=0.0_fp
  wz2d=0.0_fp
  
  do ii=1,inum
     index2d(ii)=ii
     zmvol(ii)=bmvol(ii)
  end do
  
  ict=0
  xsurf(1)=0.0_fp
  xsurf(nxsurf)=1.0_fp
  ictx=1
  do ix=1,nznbmrf
     zx=(ix-0.5_fp)/nznbmrf
     if(ix.gt.1) then
        ictx=ictx+1
        xsurf(ictx)=(zx+xsurf(ictx-1))/2.0_fp
     end if
     ictx=ictx+1
     xsurf(ictx)=zx
     ith1=nthzsm(lcentr+ix-2)+1
     ith2=nthzsm(lcentr+ix-1)
     nth=ith2-ith1+1
     do ith=1,nth
        zth=thbdy0+(ith-0.5_fp)*(thbdy1-thbdy0)/nth
        zth1=thbdy0+(ith-1)*(thbdy1-thbdy0)/nth
        zth2=thbdy0+ith*(thbdy1-thbdy0)/nth
        ict=ict+1
        zindr(ict)=dble(ict)
        zindz(ict)=dble(ict)
        zindr(ict+inum)=dble(ict+inum)
        zindz(ict+inum)=dble(ict+inum)
        zindr(ict+inum*2)=dble(ict+inum*2)
        zindz(ict+inum*2)=dble(ict+inum*2)
        zindr(ict+inum*3)=dble(ict+inum*3)
        zindz(ict+inum*3)=dble(ict+inum*3)
        x2d(ict)=zx
        th2d(ict)=zth
        wx2d(ict)=zx
        wth2d(ict)=zth
        wx2d(ict+inum)=xsurf(ictx-1)
        wth2d(ict+inum)=zth
        wx2d(ict+2*inum)=xsurf(ictx-1)
        wth2d(ict+2*inum)=zth1
        wx2d(ict+3*inum)=xsurf(ictx-1)
        wth2d(ict+3*inum)=zth2
        call rymap(x2d(ict),th2d(ict),r2d(ict),z2d(ict))
        wr2d(ict)=r2d(ict)
        wz2d(ict)=z2d(ict)
        call rymap(wx2d(ict+inum),wth2d(ict+inum),wr2d(ict+inum), &
             wz2d(ict+inum))
        call rymap(wx2d(ict+2*inum),wth2d(ict+2*inum), &
             wr2d(ict+2*inum),wz2d(ict+2*inum))
        call rymap(wx2d(ict+3*inum),wth2d(ict+3*inum), &
             wr2d(ict+3*inum),wz2d(ict+3*inum))
     end do
  end do
  if(allocated(bdens2)) then
     call tr_dsort(wr2d,zindr,inum*4,2)
     call tr_dsort(wz2d,zindz,inum*4,2)
     do ix=1,inum*4
        do iz=1,inum*4
           if(zindr(ix).eq.zindz(iz)) then
              ict=int(zindz(iz))
              i_mult=(ict-1)/inum
              ict=ict-i_mult*inum
              do is=1,nfast
                 wrk_bd2(ix,iz,is)=bdens2(ict,is)
              end do
           end if
        end do
     end do
     !
     !  plot the data
     !
     !         do  is=1,nfast
     ! 	   call r8_grf3f1(wr2d,wrk_bd2(1:inum*4,1:inum*4,is),wz2d,wk,
     !     >           inum*4,inum*4,inum*4,
     !     >     'cm','#/cm**3','cm', 'cm','#/cm**3','cm',
     !     >     'bdens2(R,Z)','#/cm**3',1)
     !
     !         end do
     
  end if
  ictx=ictx+1
  xsurf(ictx)=1
  
  ! define polar angle grid -- top half plasma only for up-down symmetric
  ! cases...
  
  do ith=1,nthsurf
     thsurf(ith)=(thbdy0*(nthsurf-ith)+thbdy1*(ith-1))/(nthsurf-1)
  end do
  
  ! fill in the equilibrium geometry
  
  do ix=1,nxsurf
     if(ix.eq.1) then
        call rymap(xsurf(1),thsurf(1), &
             rsurf(1,1),zsurf(1,1))
        rsurf(2:nthsurf,1)=rsurf(1,1)
        zsurf(2:nthsurf,1)=zsurf(1,1)
     else
        do ith=1,nthsurf
           call rymap(xsurf(ix),thsurf(ith), &
                rsurf(ith,ix),zsurf(ith,ix))
        end do
     end if
  end do
  
  !  define stuff not dependent on specie index
  call cdf_define(icdf,'TRANSP_RUNID',trim(runid))
  call cdf_define(icdf,'TOKAMAK',tdev)
  call cdf_define(icdf,'TIME',zta)
  call cdf_define(icdf,'DT_AVG',zdta)
  call cdf_define(icdf,'HINT',hint)
  call cdf_define(icdf,'N_SPECIES',itot)
  do ii=1,itot
     slname=' '
     write(slname,'("SPECIES_",i1)') ii
     call cdf_define(icdf,slname,trim(skeep(ii)(2:)))
  end do
  call cdf_define(icdf,'SYMFLAG',isym)
  call cdf_define(icdf,'GCFLAG',igc)
  call cdf_define(icdf,'PTCLFLAG',iptcl)
  call cdf_define(icdf,'N_ZONEROWS',nznbmrf)
  call cdf_define(icdf,'N_2DZONES',inum)
  call cdf_define(icdf,'NTHZSM',nthzsm(1:lcentr+nznbmrf-1))
  call cdf_define(icdf,'INDEX_2D',index2d)
  call cdf_define(icdf,'X2D',x2d)
  call cdf_define(icdf,'TH2D',th2d)
  call cdf_define(icdf,'R2D',r2d)
  call cdf_define(icdf,'Z2D',z2d)
  call cdf_define(icdf,'BMVOL',zmvol(1:inum))
  call cdf_define(icdf,'NXSURF',nxsurf)
  call cdf_define(icdf,'NTHSURF',nthsurf)
  call cdf_define(icdf,'XSURF',xsurf)
  call cdf_define(icdf,'THSURF',thsurf)
  call cdf_define(icdf,'RSURF',rsurf)
  call cdf_define(icdf,'ZSURF',zsurf)
  
  !  define stuff dependent on species
  do ns=1,nsfast
     if(.not.present(ns)) cycle
     call cdf_define(icdf,'NTOT'//trim(slabel(ns)),znum(ns))
     ine=nznbmea(ns)
     ina=nznbma
     call cdf_define(icdf,'NE'//trim(slabel(ns)),ine)
     call cdf_define(icdf,'NA'//trim(slabel(ns)),ina)
     call cdf_define(icdf,'E'//trim(slabel(ns)),zefbm(1:ine,ns))
     call cdf_define(icdf,'EB'//trim(slabel(ns)),zefbmb(1:ine+1,ns))
     call cdf_define(icdf,'A'//trim(slabel(ns)),zksid(1:ina,ns))
     call cdf_define(icdf,'AB'//trim(slabel(ns)),zksidb(1:ina+1,ns))
     !  define fbm this way, to prevent redundant copy of large array...
     idm_fbm(1)=ine
     idm_fbm(2)=ina
     call cdf_define(icdf,'F'//trim(slabel(ns)),idm_fbm,'R8')
  end do
  
  !  additional get_fbm_kernel module data & mapping to orig. indexing
  call cdf_define(icdf,'imap',imap)
  call cdf_define(icdf,'ndep0_use',ndep0_use)
  call cdf_define(icdf,'nmom',nmom)
  call cdf_define(icdf,'nshot',nshot)
  call cdf_define(icdf,'nzones',nzones)
  call cdf_define(icdf,'nznbma',nznbma)
  call cdf_define(icdf,'nznbmea',nznbmea(1:nsfast))
  
  if(allocated(abeams)) then
     call cdf_define(icdf,'abeams',abeams)
     call cdf_define(icdf,'xzbeams',xzbeams)
     call cdf_define(icdf,'ilfprod',ilfprod)
     call cdf_define(icdf,'ilbeam',ilbeam)
  end if
  if(allocated(bdens2)) then
     call cdf_define(icdf,'bdens2',bdens2(1:inum,1:nfast))
     call cdf_define(icdf,'eba2pl',eba2pl(1:inum,1:nfast))
     call cdf_define(icdf,'eba2pp',eba2pp(1:inum,1:nfast))
  end if
  if(nldeptrk) then
     call cdf_define(icdf,'iversion_deptrk',iversion_deptrk)
     call cdf_define(icdf,'efusion_sce',efusion_sce)
     call cdf_define(icdf,'fusion_sce',fusion_sce(1:inum,1:mireact))
     if(allocated(ilfusion_sce)) deallocate(ilfusion_sce)
     allocate(ilfusion_sce(mireact,nfast))
     do is=1,nfast
        do ir=1,mireact
           if(nlfusion_sce(ir,is)) then
              ilfusion_sce(ir,is)=1
           else
              ilfusion_sce(ir,is)=0
           end if
        end do
     end do
     call cdf_define(icdf,'ilfusion_sce',ilfusion_sce)
     
     if(idepb.gt.0) then
        call cdf_define(icdf,'nstep_gtime',nstep_gtime)
        call cdf_define(icdf,'ntracks',ntracks)
        call cdf_define(icdf,'nzmax_track',nzmax_track)
        call cdf_define(icdf,'track_gtime1',track_gtime1)
        call cdf_define(icdf,'track_gtime2',track_gtime2)
     end if
     ns_use=0
     inz_p= 1
     do ns=1,nsfast
        if(nldep_beam(ns)) then
           ns_use=ns_use+1
           inz_track=sum(nztrack(:,ns_use),dim=1)
           if(inz_track.eq.0) cycle
           inz_c = inz_p + inz_track -1
           call cdf_define(icdf,'nztrack'//trim(slabel(ns)), &
                nztrack(:,ns_use))
           call cdf_define(icdf,'nstat_track_xja'//trim(slabel(ns)), &
                nstat_track_xja(inz_p:inz_c))
           call cdf_define(icdf,'ntrack_nbsce'//trim(slabel(ns)), &
                ntrack_nbsce(:,ns_use))
           call cdf_define(icdf,'ntrack_nbien'//trim(slabel(ns)), &
                ntrack_nbien(:,ns_use))
           call cdf_define(icdf,'track_phiin'//trim(slabel(ns)), &
                track_phiin(:,ns_use))
           call cdf_define(icdf,'track_phiout'//trim(slabel(ns)), &
                track_phiout(:,ns_use))
           call cdf_define(icdf,'track_rin'//trim(slabel(ns)), &
                track_rin(:,ns_use))
           call cdf_define(icdf,'track_rout'//trim(slabel(ns)), &
                track_rout(:,ns_use))
           call cdf_define(icdf,'track_zin'//trim(slabel(ns)), &
                track_zin(:,ns_use))
           call cdf_define(icdf,'track_zout'//trim(slabel(ns)), &
                track_zout(:,ns_use))
           call cdf_define(icdf,'track_einj'//trim(slabel(ns)), &
                track_einj(:,ns_use))
           call cdf_define(icdf,'track_sinj'//trim(slabel(ns)), &
                track_sinj(:,ns_use))
           call cdf_define(icdf,'track_time'//trim(slabel(ns)), &
                track_time(:,ns_use))
           call cdf_define(icdf,'track_xl'//trim(slabel(ns)), &
                track_xl(inz_p:inz_c))
           call cdf_define(icdf,'track_pdep'//trim(slabel(ns)), &
                track_pdep(inz_p:inz_c))
           call cdf_define(icdf,'track_vpv_dep'//trim(slabel(ns)), &
                track_vpv_dep(inz_p:inz_c))
           call cdf_define(icdf,'track_dr_flr'//trim(slabel(ns)), &
                track_dr_flr(inz_p:inz_c))
           call cdf_define(icdf,'track_dz_flr'//trim(slabel(ns)), &
                track_dz_flr(inz_p:inz_c))
           call cdf_define(icdf,'track_de_flr'//trim(slabel(ns)), &
                track_de_flr(inz_p:inz_c))
           call cdf_define(icdf,'track_dvpv_flr'//trim(slabel(ns)), &
                track_dvpv_flr(inz_p:inz_c))
           inz_p = inz_c +1
        end if
     end do
  end if
  if(nllost) then
     call cdf_define(icdf,'gtime1_lost',gtime1_lost)
     call cdf_define(icdf,'gtime2_lost',gtime2_lost)
     call cdf_define(icdf,'TIMELST',gtime_lost(1:mlost_data))
     call cdf_define(icdf,'NBSCELST',ynbsce_lost(1:mlost_data))
     call cdf_define(icdf,'NBIENLST',ynbien_lost(1:mlost_data))
     call cdf_define(icdf,'WGHTLST',wght_lost(1:mlost_data))
     call cdf_define(icdf,'RMJION',rmjion_lost(1:mlost_data))
     call cdf_define(icdf,'ZYION',yyion_lost(1:mlost_data))
     call cdf_define(icdf,'ZELST',yze_lost(1:mlost_data))
     call cdf_define(icdf,'XKSIDLST',xksid_lost(1:mlost_data))
     call cdf_define(icdf,'GOOSELST',goose_lost(1:mlost_data))
     call cdf_define(icdf,'SPECIESLST',species_lost(1:mlost_data))
  end if
  if((nllost.or.nldeptrk).and. nbeam.gt.0) then
     call cdf_define(icdf,'NLCO',znlco(1:nbeam))
     call cdf_define(icdf,'PINJA',pinja(1:nbeam))
     call cdf_define(icdf,'XBZETA',xbzeta(1:nbeam))
     call cdf_define(icdf,'RTCENA',rtcena(1:nbeam))
     call cdf_define(icdf,'XLBAPA',xlbapa(1:nbeam))
     call cdf_define(icdf,'XLBTNA',xlbtna(1:nbeam))
     call cdf_define(icdf,'nsjdotb',nsjdotb)
     call cdf_define(icdf,'nsnccw',nsnccw)
  end if
  
  if(allocated(omegag)) then
     call cdf_define(icdf,'omegag',omegag)
     if(nlsym) then
        call cdf_define(icdf,'rmom0',rmom0)
        call cdf_define(icdf,'rmom',rmom)
        call cdf_define(icdf,'ymom',ymom)
     else
        call cdf_define(icdf,'rmc',rmc)
        call cdf_define(icdf,'ymc',ymc)
     end if
     call cdf_define(icdf,'xitr',xitr)
  end if
  
  !  addl labeling/attributes
  call cdf_setatt(icdf,'SYMFLAG','=1 means up-down symmetry')
  call cdf_setatt(icdf,'GCFLAG', &
       '=1 means guiding center F; 0 means ptcl F')
  call cdf_setatt(icdf,'N_ZONEROWS','no. of f_i zone rows')
  call cdf_setatt(icdf,'N_2DZONES','no. of f_i zones')
  call cdf_setatt(icdf,'INDEX_2D','f_i zone index')
  call cdf_setatt(icdf,'X2D','x of f_i zones')
  call cdf_setatt(icdf,'TH2D','theta of f_i zones')
  call cdf_setatt(icdf,'R2D','R of f_i zones','cm')
  call cdf_setatt(icdf,'Z2D','Z of f_i zones','cm')
  call cdf_setatt(icdf,'BMVOL','volume of f_i zones','cm^3')
  call cdf_setatt(icdf,'NXSURF','size of xsurf grid')
  call cdf_setatt(icdf,'NTHSURF','size of thsurf grid')
  call cdf_setatt(icdf,'XSURF','xsurf grid (for flux surfaces)')
  call cdf_setatt(icdf,'THSURF','thsurf grid (for flux surfaces)')
  call cdf_setatt(icdf,'RSURF','flux surface R locations','cm')
  call cdf_setatt(icdf,'ZSURF','flux surface Z locations','cm')

  !  attributes for stuff dependent on species
  do ns=1,nsfast
     if(.not.present(ns)) cycle
     call cdf_setatt(icdf,'NTOT'//trim(slabel(ns)), &
          'total number of ions')
     call cdf_setatt(icdf,'NE'//trim(slabel(ns)), &
          'number of energy bins')
     call cdf_setatt(icdf,'NA'//trim(slabel(ns)), &
          'number of pitch bins')
     call cdf_setatt(icdf,'E'//trim(slabel(ns)), &
          'energy bin ctrs','eV')
     call cdf_setatt(icdf,'EB'//trim(slabel(ns)), &
          'energy bin bdys','eV')
     call cdf_setatt(icdf,'A'//trim(slabel(ns)), &
          'pitch bin ctrs (v||/v)')
     call cdf_setatt(icdf,'AB'//trim(slabel(ns)), &
          'pitch bin bdys (v||/v)')
     call cdf_setatt(icdf,'F'//trim(slabel(ns)), &
          'distribution function','#/cm^3/eV/d(omega/4pi)')
  end do
  if(nllost) then
     
     call cdf_setatt(icdf,'TIMELST','time when ptcl lost')
     call cdf_setatt(icdf,'NBSCELST','beam# from which lost ptcl, 0 &
          means fusion product')
     call cdf_setatt(icdf,'NBIENLST','energy fraction of lost ptcl')
     call cdf_setatt(icdf,'WGHTLST','weight of lost ptcl')
     call cdf_setatt(icdf,'RMJION','major radius at what ptcl lost')
     call cdf_setatt(icdf,'ZYION','Z at what ptcl lost')
     call cdf_setatt(icdf,'ZELST','energy lost ptcl')
     call cdf_setatt(icdf,'XKSIDLST','Vll/V of lost ptcl')
     call cdf_setatt(icdf,'GOOSELST','1.00000 means prompt loss')
     call cdf_setatt(icdf,'SPECIESLST','species')
  end if
  if((nllost.or.nldeptrk).and. nbeam.gt.0) then
     
     call cdf_setatt(icdf,'NLCO',' beam co-, or counter-injector')
     call cdf_setatt(icdf,'PINJA','beams injected power')
     call cdf_setatt(icdf,'XBZETA','beams toroidal angle, deg')
     call cdf_setatt(icdf,'RTCENA', &
          'tangency major rad., cm, of centerline of beam')
     call cdf_setatt(icdf,'XLBAPA', &
          'distance beam source to aperture, cm')
     call cdf_setatt(icdf,'XLBTNA', &
          'distance beam source tot angency point, cm')
     call cdf_setatt(icdf,'nsjdotb', &
          '+1 if Ip and Bt are parallel, -1 if anti-pll')
     call cdf_setatt(icdf,'nsnccw', &
          '+1 if Bphi is CCW looking down on tokamak from above')
     
  end if
  
  !  write stuff not dependent on specie index
  call cdf_write(icdf,'TRANSP_RUNID',trim(runid))
  call cdf_write(icdf,'TOKAMAK',tdev)
  call cdf_write(icdf,'TIME',zta)
  call cdf_write(icdf,'DT_AVG',zdta)
  call cdf_write(icdf,'HINT',hint)
  call cdf_write(icdf,'N_SPECIES',itot)
  do ii=1,itot
     slname=' '
     write(slname,'("SPECIES_",i1)') ii
     call cdf_write(icdf,slname,trim(skeep(ii)(2:)))
  end do
  call cdf_write(icdf,'SYMFLAG',isym)
  call cdf_write(icdf,'GCFLAG',igc)
  call cdf_write(icdf,'PTCLFLAG',iptcl)
  call cdf_write(icdf,'N_ZONEROWS',nznbmrf)
  call cdf_write(icdf,'N_2DZONES',inum)
  call cdf_write(icdf,'NTHZSM',nthzsm(1:lcentr+nznbmrf-1))
  call cdf_write(icdf,'INDEX_2D',index2d)
  call cdf_write(icdf,'X2D',x2d)
  call cdf_write(icdf,'TH2D',th2d)
  call cdf_write(icdf,'R2D',r2d)
  call cdf_write(icdf,'Z2D',z2d)
  call cdf_write(icdf,'BMVOL',zmvol(1:inum))
  call cdf_write(icdf,'NXSURF',nxsurf)
  call cdf_write(icdf,'NTHSURF',nthsurf)
  call cdf_write(icdf,'XSURF',xsurf)
  call cdf_write(icdf,'THSURF',thsurf)
  call cdf_write(icdf,'RSURF',rsurf)
  call cdf_write(icdf,'ZSURF',zsurf)
  
  !  write stuff dependent on species
  do ns=1,nsfast
     if(.not.present(ns)) cycle
     call cdf_write(icdf,'NTOT'//trim(slabel(ns)),znum(ns))
     ine=nznbmea(ns)
     ina=nznbma
     call cdf_write(icdf,'NE'//trim(slabel(ns)),ine)
     call cdf_write(icdf,'NA'//trim(slabel(ns)),ina)
     call cdf_write(icdf,'E'//trim(slabel(ns)),zefbm(1:ine,ns))
     call cdf_write(icdf,'EB'//trim(slabel(ns)),zefbmb(1:ine+1,ns))
     call cdf_write(icdf,'A'//trim(slabel(ns)),zksid(1:ina,ns))
     call cdf_write(icdf,'AB'//trim(slabel(ns)),zksidb(1:ina+1,ns))
     
     allocate(zfbm(ine,ina,inum))
     !MG
     !MG writing fbm_ptcl has been temprorarily disabeled, April 15,2011
     do ie=1,ine
        do ia=1,ina
           do ith=1,inum
              zfbm(ie,ia,ith) = fbm_fnc(ie,ia,ith,ns)
           end do
        end do
     end do
     !     MG
     call cdf_write(icdf,'F'//trim(slabel(ns)),zfbm)
     deallocate(zfbm)
     
  end do
  
  !  additional get_fbm_kernel module data & mapping to orig. indexing
  call cdf_write(icdf,'imap',imap)
  call cdf_write(icdf,'ndep0_use',ndep0_use)
  call cdf_write(icdf,'nmom',nmom)
  call cdf_write(icdf,'nshot',nshot)
  call cdf_write(icdf,'nzones',nzones)
  call cdf_write(icdf,'nznbma',nznbma)
  call cdf_write(icdf,'nznbmea',nznbmea(1:nsfast))
  
  if(allocated(abeams)) then
     call cdf_write(icdf,'abeams',abeams)
     call cdf_write(icdf,'xzbeams',xzbeams)
     call cdf_write(icdf,'ilfprod',ilfprod)
     call cdf_write(icdf,'ilbeam',ilbeam)
  end if
  if(allocated(bdens2)) then
     call cdf_write(icdf,'bdens2',bdens2(1:inum,1:nfast))
     call cdf_write(icdf,'eba2pl',eba2pl(1:inum,1:nfast))
     call cdf_write(icdf,'eba2pp',eba2pp(1:inum,1:nfast))
  end if
  if(nldeptrk) then
     call cdf_write(icdf,'iversion_deptrk',iversion_deptrk)
     call cdf_write(icdf,'efusion_sce',efusion_sce)
     call cdf_write(icdf,'fusion_sce',fusion_sce(1:inum,1:mireact))
     call cdf_write(icdf,'ilfusion_sce',ilfusion_sce)
     
     if(idepb.gt.0) then
        call cdf_write(icdf,'nstep_gtime',nstep_gtime)
        call cdf_write(icdf,'ntracks',ntracks)
        call cdf_write(icdf,'nzmax_track',nzmax_track)
        call cdf_write(icdf,'track_gtime1',track_gtime1)
        call cdf_write(icdf,'track_gtime2',track_gtime2)
     end if
     
     ns_use=0
     inz_p= 1
     
     do ns=1,nsfast
        inz_track = 0
        if(nldep_beam(ns)) then
           ns_use=ns_use+1
           inz_track=sum(nztrack(:,ns_use),dim=1)
           if(inz_track.eq.0) cycle
           call cdf_write(icdf,'nztrack'//trim(slabel(ns)), &
                nztrack(:,ns_use))
           inz_c = inz_p + inz_track -1
           call cdf_write(icdf,'nstat_track_xja'//trim(slabel(ns)), &
                nstat_track_xja(inz_p:inz_c))
           if(iversion_deptrk.eq.3.) then
              call cdf_write(icdf,'ntrack_nbsce'//trim(slabel(ns)), &
                   ntrack_nbsce(:,ns_use))
              call cdf_write(icdf,'ntrack_nbien'//trim(slabel(ns)), &
                   ntrack_nbien(:,ns_use))
           end if
           call cdf_write(icdf,'track_phiin'//trim(slabel(ns)), &
                track_phiin(:,ns_use))
           call cdf_write(icdf,'track_phiout'//trim(slabel(ns)), &
                track_phiout(:,ns_use))
           call cdf_write(icdf,'track_rin'//trim(slabel(ns)), &
                track_rin(:,ns_use))
           call cdf_write(icdf,'track_rout'//trim(slabel(ns)), &
                track_rout(:,ns_use))
           call cdf_write(icdf,'track_zin'//trim(slabel(ns)), &
                track_zin(:,ns_use))
           call cdf_write(icdf,'track_zout'//trim(slabel(ns)), &
                track_zout(:,ns_use))
           call cdf_write(icdf,'track_einj'//trim(slabel(ns)), &
                track_einj(:,ns_use))
           call cdf_write(icdf,'track_sinj'//trim(slabel(ns)), &
                track_sinj(:,ns_use))
           call cdf_write(icdf,'track_time'//trim(slabel(ns)), &
                track_time(:,ns_use))
           call cdf_write(icdf,'track_xl'//trim(slabel(ns)), &
                track_xl(inz_p:inz_c))
           call cdf_write(icdf,'track_pdep'//trim(slabel(ns)), &
                track_pdep(inz_p:inz_c))
           call cdf_write(icdf,'track_vpv_dep'//trim(slabel(ns)), &
                track_vpv_dep(inz_p:inz_c))
           call cdf_write(icdf,'track_dr_flr'//trim(slabel(ns)), &
                track_dr_flr(inz_p:inz_c))
           call cdf_write(icdf,'track_dz_flr'//trim(slabel(ns)), &
                track_dz_flr(inz_p:inz_c))
           call cdf_write(icdf,'track_de_flr'//trim(slabel(ns)), &
                track_de_flr(inz_p:inz_c))
           call cdf_write(icdf,'track_dvpv_flr'//trim(slabel(ns)), &
                track_dvpv_flr(inz_p:inz_c))
           inz_p = inz_c +1
        end if
     end do
  end if
  if(nllost) then
     ! if there are lost orbits write these data
     call cdf_write(icdf,'gtime1_lost',gtime1_lost)
     call cdf_write(icdf,'gtime2_lost',gtime2_lost)
     call cdf_write(icdf,'TIMELST',gtime_lost(1:mlost_data))
     call cdf_write(icdf,'NBSCELST',ynbsce_lost(1:mlost_data))
     call cdf_write(icdf,'NBIENLST',ynbien_lost(1:mlost_data))
     call cdf_write(icdf,'WGHTLST',wght_lost(1:mlost_data))
     call cdf_write(icdf,'RMJION',rmjion_lost(1:mlost_data))
     call cdf_write(icdf,'ZYION',yyion_lost(1:mlost_data))
     call cdf_write(icdf,'ZELST',yze_lost(1:mlost_data))
     call cdf_write(icdf,'XKSIDLST',xksid_lost(1:mlost_data))
     call cdf_write(icdf,'GOOSELST',goose_lost(1:mlost_data))
     call cdf_write(icdf,'SPECIESLST',species_lost(1:mlost_data))
  end if
  
  if((nllost.or.nldeptrk).and. nbeam.gt.0) then
     ! write  data for the direction of beam injection, and injected power
     
     call cdf_write(icdf,'NLCO',znlco(1:nbeam))
     call cdf_write(icdf,'PINJA',pinja(1:nbeam))
     call cdf_write(icdf,'XBZETA',xbzeta(1:nbeam))
     call cdf_write(icdf,'RTCENA',rtcena(1:nbeam))
     call cdf_write(icdf,'XLBAPA',xlbapa(1:nbeam))
     call cdf_write(icdf,'XLBTNA',xlbtna(1:nbeam))
     call cdf_write(icdf,'nsjdotb',nsjdotb)
     call cdf_write(icdf,'nsnccw',nsnccw)
  end if
  
  if(allocated(omegag)) then
     call cdf_write(icdf,'omegag',omegag)
     if(nlsym) then
        call cdf_write(icdf,'rmom0',rmom0)
        call cdf_write(icdf,'rmom',rmom)
        call cdf_write(icdf,'ymom',ymom)
     else
        call cdf_write(icdf,'rmc',rmc)
        call cdf_write(icdf,'ymc',ymc)
     end if
     call cdf_write(icdf,'xitr',xitr)
  end if
  
  call cdf_close(icdf)
  
  deallocate(index2d,x2d,th2d,r2d,z2d,zmvol)
  deallocate(xsurf,thsurf,rsurf,zsurf)
  deallocate(zksid,zksidb,zefbm,zefbmb)
  
  ns=ins_save
  return
end subroutine fbm_wnetcdf_fp
