subroutine nubeam_comp_ps_ini(ss, ilun_msgs, &
     nbi_init, gas_init, inzon_nbi, inzon_gas)

  use plasma_state_mod
  use nbi_types
  implicit NONE

  ! check/initialize nubeam_comp_exec output state

  ! IF the state is already initialized and the grid sizes it contains
  ! do not match, write a warning message.

  type (plasma_state) :: ss

  integer, intent(in) :: ilun_msgs   ! I/O unit for messages

  logical, intent(in) :: nbi_init    ! .TRUE. to check/initialize NBI & FUS
  logical, intent(in) :: gas_init    ! .TRUE. to check/initialize GAS

  integer, intent(in) :: inzon_nbi   ! #zones for NBI & FUS, if known
  integer, intent(in) :: inzon_gas   ! #zones for GAS, if known

  !  (nbi_init).AND.(inzon_nbi==0) means: query NUBEAM for grid size
  !  (gas_init).AND.(inzon_gas==0) means: match NUBEAM grid size

  !-------------------------------
  type(nbitype_grid) :: zgrid
  integer :: jnzon_nbi,jnzon_gas
  integer :: ierr
  !-------------------------------

  if(nbi_init) then
     if(inzon_nbi.gt.0) then
        jnzon_nbi=inzon_nbi
     else
        call nbi_get_grid(zgrid)
        jnzon_nbi = zgrid%nzones
     endif

     if(jnzon_nbi.eq.0) then
        write(ilun_msgs,*) ' ?nubeam_comp_ps_ini: found #zones = 0 for NBI, FUS'
     else

        if(ss%nrho_nbi.eq.0) then
           ss%nrho_nbi = jnzon_nbi + 1
           call ps_alloc_plasma_state(ierr, state=ss)
           if(ierr.ne.0) then
              write(ilun_msgs,*) ' ?nubeam_comp_ps_ini: ps_alloc_plasma_state failure on: nrho_nbi.'
           else
              call genrho(ss%nrho_nbi,ss%rho_nbi)
           endif
        else
           if(ss%nrho_nbi.ne.(jnzon_nbi+1)) then
              write(ilun_msgs,*) ' ?nubeam_comp_ps_ini: #zones in PS(NBI) does not match input: ',ss%nrho_nbi-1,jnzon_nbi
           endif
        endif

        ! note: FUS and NBI grids expected to match...

        if(ss%nrho_fus.eq.0) then
           ss%nrho_fus = jnzon_nbi + 1
           call ps_alloc_plasma_state(ierr, state=ss)
           if(ierr.ne.0) then
              write(ilun_msgs,*) ' ?nubeam_comp_ps_ini: ps_alloc_plasma_state failure on: nrho_fus.'
           else
              call genrho(ss%nrho_fus,ss%rho_fus)
           endif
        else
           if(ss%nrho_fus.ne.(jnzon_nbi+1)) then
              write(ilun_msgs,*) ' ?nubeam_comp_ps_ini: #zones in PS(FUS) does not match input: ',ss%nrho_fus-1,jnzon_nbi
           endif
        endif

     endif
  endif

  if(gas_init) then
     if(inzon_gas.gt.0) then
        jnzon_gas = inzon_gas
     else
        jnzon_gas = max(ss%nrho_nbi,ss%nrho_fus) - 1
     endif

     if(jnzon_gas.le.0) then
        write(ilun_msgs,*) ' ?nubeam_comp_ps_ini: found #zones = 0 for GAS'

     else
        if(ss%nrho_gas.eq.0) then
           ss%nrho_gas = jnzon_gas + 1
           call ps_alloc_plasma_state(ierr, state=ss)
           if(ierr.ne.0) then
              write(ilun_msgs,*) ' ?nubeam_comp_ps_ini: ps_alloc_plasma_state failure on: nrho_gas.'
           else
              call genrho(ss%nrho_gas,ss%rho_gas)
           endif
        else
           if(ss%nrho_gas.ne.(jnzon_gas+1)) then
              write(ilun_msgs,*) ' ?nubeam_comp_ps_ini: #zones in PS(GAS) does not match input: ',ss%nrho_gas-1,jnzon_gas
           endif
        endif
     endif
  endif

CONTAINS
  subroutine genrho(inum,zrho)

    ! make evenly spaced rho grid

    integer, intent(in) :: inum  ! #bdys
    real*8, dimension(:) :: zrho ! rho grid out

    !---------------------------
    real*8, parameter :: ZERO = 0.0d0
    real*8, parameter :: ONE = 1.0d0
    integer :: irho
    !---------------------------

    zrho(1)=ZERO
    zrho(inum)=ONE
    do irho=2,inum-1
       zrho(irho) = (irho-1)*ONE/(inum-1)
    enddo

  end subroutine genrho

end subroutine nubeam_comp_ps_ini
