subroutine nubeam_mcgrid_define(inphi,iudsym,inth0,inznbmri,inznbmr,id)
 
  !  define a 2d irregular "mc" polar grid -- attribute ownership "nubeam"

  !  fewer zones near the axis, more zones towards the edge, each
  !  zone has roughly the same volume... see xplasma mcgrid documentation

  use xplasma_obj_instance
 
  implicit none
 
  integer, intent(in) :: inphi      ! no. of phi zones, 1 for axisymmetry
  integer, intent(in) :: iudsym     ! =1: updown symmetry; =2: no symmetry
  integer, intent(in) :: inth0      ! no. of zones spanning theta=[0,pi] @axis
  integer, intent(in) :: inznbmri   ! no. of radial zones spanning [axis,bdy]
  integer, intent(out) :: inznbmr   ! total no. of radial zones
  !  inznbmr = inznbmri + inznbmri/3 -- extension of grid beyond plasma
  !    boundary; usage is optional.
 
  integer, intent(out) :: id        ! grid id, positive integer; 0 means error
 
  !-----------------------------------------
  integer :: ierr
  integer :: loc_nphi,loc_udsym,loc_nth0,loc_nznbmri,loc_nznbmr

  logical :: idiff
  !-----------------------------------------

  idiff = .true.
  call mcgrid_find(id)
  if(id.gt.0) then
     do
        call mcgrid_getsym(id,loc_udsym)
        if(loc_udsym.ne.iudsym) exit

        call mcgrid_getnuma(id,loc_nth0,loc_nphi)
        if(loc_nphi.ne.inphi) exit
        if(loc_nth0.ne.inth0) exit

        call mcgrid_getnumr(id,loc_nznbmr,loc_nznbmri)
        if(loc_nznbmri.ne.inznbmri) exit

        inznbmr=loc_nznbmr
        idiff = .false.
        exit
     enddo
  endif

  if(idiff) then
     call mcgrid_named_define('nubeam',inphi,iudsym,inth0,inznbmri,inznbmr,id)
  endif
  return

end subroutine nubeam_mcgrid_define
