subroutine rf_wavecode_output(pstate,ierr)

  !  fetch RF wavecode outputs looping over all frequencies and Nphi values

  use rfxqlo_interface
  use plasma_state_mod
  implicit NONE

  type (plasma_state) :: pstate   ! store outputs in this plasma state object
  integer, intent(out) :: ierr    ! completion code, 0=OK

  !--------------------------------------
  type (rfxqlo_dims) :: zdims
  type (rfxqlo_antenna) :: zant
  type (rfxqlo_rfgrid) :: zgrid
  type (rfxqlo_ion) :: zions    ! species data

  type(rfxqlo_outpro) :: zpro

  integer :: lunmsg = 6

  REAL*8, PARAMETER :: ZERO = 0.0d0

  integer i,j,id_rho,inrho,idum,iant,n_antenna,inphi,inions,itype
  integer :: ierloc,iersum,i1,i2,istate
  logical :: isym

  integer, dimension(:), allocatable :: jj
  real*8, dimension(:), allocatable :: zwkbuf,zpzon,zrho
  real*8 :: zcpl,zcpl_sum,zpext,zp1

  !--------------------------------------
  !  basic dimension info
  call rfxqlo_dim_get(zdims,ierr)
  if(ierr.ne.0) then
     write(lunmsg,*) ' ?rf_wavecode_output: rfxqlo_dim_get failure.'
     return
  endif

  !  # of antennas

  n_antenna = zdims%n_antenna
  if(n_antenna.eq.0) then
     ierr=2
     write(lunmsg,*) ' ?rf_wavecode_output: no ICRF antennas found.'
     return
  endif

  !--------------------------------------
  !  grid info

  call rfxqlo_rfgrid_get(zgrid,ierr)
  if(ierr.ne.0) then
     write(lunmsg,*) ' ?rf_wavecode_output: rfxqlo_rfgrid_get failure.'
     return
  endif

  !  "qlo" grid -- historically, grid used by FP code to model minority
  !  fast ion species; here, should match rho_icrf(...) in pstate

  id_rho = zgrid%id_xqlo
  call eq_ngrid(id_rho,inrho)

  if(pstate%nrho_icrf.eq.0) then
     pstate%nrho_icrf = inrho
     call ps_alloc_plasma_state(ierr, state=pstate)
     if(ierr.ne.0) then
        write(lunmsg,*) ' ?rf_wavecode_output: ps_alloc_plasma_state error #1.'
        return
     endif
     call eq_grid(id_rho,pstate%rho_icrf,pstate%nrho_icrf,idum,ierr)
     if(ierr.ne.0) then
        write(lunmsg,*) ' ?rf_wavecode_output: rho_icrf grid fetch error.'
        return
     endif
  else if(pstate%nrho_icrf.ne.inrho) then
     ierr=4
     write(lunmsg,*) ' ?rf_wavecode_output: QLO and state grid sizes mismatch:'
     write(lunmsg,*) '  inrho(QLO)=',inrho,' pstate%nrho_icrf=',pstate%nrho_icrf
     return
  else
     allocate(zrho(inrho))
     call eq_grid(id_rho,zrho,inrho,idum,ierr)
     if(ierr.ne.0) then
        write(lunmsg,*) ' ?rf_wavecode_output: rho_icrf grid fetch error.'
        deallocate(zrho)
        return
     endif
     if(maxval(abs(zrho-pstate%rho_icrf)).gt.1.0d-5) then
        write(lunmsg,*) ' ?rf_wavecode_output: pstate%rho_icrf and rfxqlo "xqlo" grids do not match.'
        ierr=5
        deallocate(zrho)
        return
     endif
  endif

  !--------------------------------------
  ! check number of antennas...

  if(pstate%nicrf_src.eq.0) then
     pstate%nicrf_src = n_antenna
     call ps_alloc_plasma_state(ierr, state=pstate)
     if(ierr.ne.0) then
        write(lunmsg,*) ' ?rf_wavecode_output: ps_alloc_plasma_state error #2.'
        return
     endif
  endif

  if(n_antenna.ne.pstate%nicrf_src) then
     write(lunmsg,*) &
          ' ?rf_wavecode_output: mismatch btw rfxqlo and pstate icrf info:'
     write(lunmsg,*) '  rfxqlo n_antenna=',n_antenna,'; pstate nicrf_src: ', &
          pstate%nicrf_src
     ierr = 99
     return
  endif

  !--------------------------------------
  ! check readiness of state...

  if(.not.allocated(pstate%vol)) then
     write(lunmsg,*) &
          ' ?rf_wavecode_output: state lacks equilibrium volumes.'
     ierr = 100
     return
  endif

  if(.not.allocated(pstate%picrf_nphi_srcs)) then
     write(lunmsg,*) &
          ' ?rf_wavecode_output: state lacks picrf_nphi_srcs array.'
     ierr = 101
     return
  endif

  !--------------------------------------
  ! proceed...

  allocate(zwkbuf(inrho),zpzon(inrho-1))

  !---------------
  !  Zero the state outputs

  !  Profiles:

  pstate%picrf_nphi_srcs = ZERO
  pstate%picrf_srcs = ZERO
  pstate%picrf_totals = ZERO
  pstate%picth = ZERO

  pstate%cdicrf_nphi = ZERO
  pstate%cdicrf = ZERO
  pstate%curich = ZERO

  !  Scalar sets:

  !  Powers lost outside

  pstate%picrf_ext = ZERO
  pstate%wt_nphi_ext = ZERO  ! Nphi spectrum

  !  Total power absorbed inside plasma

  pstate%picrf_abs = ZERO
  pstate%wt_nphi_abs = ZERO  ! Nphi spectrum

  !---------------
  !  species information

  call rfxqlo_ion_get(zions,ierr)
  if(ierr.ne.0) then
     write(lunmsg,*) ' ?rf_wavecode_output: rfxqlo_ion_get failure.'
     return
  endif

  !---------------
  !  OK loop over antennas

  inions=zions%n_ions
  allocate(jj(inions))
  do i=1,inions
     jj(i)=0
     do j=1,pstate%nspec_alla
        if(zions%name_ions(i).eq.pstate%alla_name(j)) then
           jj(i)=j
        endif
     enddo
     if(jj(i).eq.0) then
        write(lunmsg,*) ' ?rf_wavecode_output: not found in pstate: '// &
             zions%name_ions(i)
        ierr=1
     endif
  enddo
  if(ierr.ne.0) return

  do iant = 1,n_antenna

     call rfxqlo_select_antenna(iant,ierr)
     if(ierr.ne.0) write(lunmsg,*) ' ...error on antenna #',iant

     call rfxqlo_antenna_get(zant,ierr)
     if(zant%power .le. ZERO) cycle

     isym = zant%sym_nphi

     zpext = ZERO
     zcpl_sum = ZERO
     iersum = 0

     do inphi = 1,zant%num_nphi
        if(zant%wnphi_vac(inphi).le.ZERO) cycle
 
        call rfxqlo_select_nphi(inphi,ierr)
        if(ierr.ne.0) write(lunmsg,*) ' ...select error on Nphi index #', &
             inphi,' for antenna #',iant

        call rfxqlo_get_mode_power(zcpl,istate,ierr)
        if(ierr.ne.0) write(lunmsg,*) ' ...mode power error on Nphi index #', &
             inphi,' for antenna #',iant
        if(zcpl.le.ZERO) cycle

        zcpl_sum = zcpl_sum + zcpl

        call rfxqlo_outpro_fetch(zpro,ierloc)
        iersum = iersum + ierloc
        if(ierloc.ne.0) cycle

        zp1 = sum(zpro%ext_pwrs)
        zpext = zpext + zp1

        if(isym) then
           i1=zant%num_nphi+1-inphi  ! stored w/o symmetry assumption in state
           i2=zant%num_nphi+inphi
           pstate%wt_nphi_ext(i1,iant) = zp1/2
           pstate%wt_nphi_ext(i2,iant) = zp1/2
           pstate%wt_nphi_abs(i1,iant) = zcpl/2
           pstate%wt_nphi_abs(i2,iant) = zcpl/2
        else
           pstate%wt_nphi_ext(inphi,iant) = zp1
           pstate%wt_nphi_abs(inphi,iant) = zcpl
        endif


        do i=1,inions + 2

           call eq_rgetf(inrho,pstate%rho_icrf,zpro%id_pwrs(i),0,zwkbuf,ierloc)
           iersum = iersum + ierloc
           if(ierloc.ne.0) cycle

           if(i.gt.inions) then
              j=0
              itype=ps_electron
           else
              j=jj(i)
              itype=pstate%alla_type(j)
           endif

           zpzon = zwkbuf(2:inrho)-zwkbuf(1:inrho-1)  ! power, W/zone

           if(isym) then
              pstate%picrf_nphi_srcs(:,i1,iant,j)= &
                   pstate%picrf_nphi_srcs(:,i1,iant,j) + zpzon/2
              pstate%picrf_nphi_srcs(:,i2,iant,j)= &
                   pstate%picrf_nphi_srcs(:,i2,iant,j) + zpzon/2
           else
              pstate%picrf_nphi_srcs(:,inphi,iant,j)= &
                   pstate%picrf_nphi_srcs(:,inphi,iant,j) + zpzon
           endif

           pstate%picrf_srcs(:,iant,j) = pstate%picrf_srcs(:,iant,j) + zpzon
           pstate%picrf_totals(:,j) = pstate%picrf_totals(:,j) + zpzon
           if((itype.eq.ps_therm_ion).or.(itype.eq.ps_impurity)) then
              pstate%picth = pstate%picth + zpzon
           endif
        enddo

        if(zpro%id_jrfe.ne.0) then
           call eq_rgetf(inrho,pstate%rho_icrf,zpro%id_jrfe,0,zwkbuf,ierloc)
           iersum = iersum + ierloc
           if(ierloc.eq.0) then
              zpzon = zwkbuf(2:inrho)-zwkbuf(1:inrho-1)  ! current, A/zone
              if(isym) then
                 pstate%cdicrf_nphi(:,i1,iant)=zpzon/2
                 pstate%cdicrf_nphi(:,i2,iant)=zpzon/2
              else
                 pstate%cdicrf_nphi(:,inphi,iant)=zpzon
              endif
              pstate%cdicrf(:,iant)=pstate%cdicrf(:,iant) + zpzon
              pstate%curich=pstate%curich + zpzon
           endif
        endif

     enddo

     if(zcpl_sum.gt.ZERO) then
        pstate%picrf_abs(iant) = zcpl_sum
        pstate%wt_nphi_abs(:,iant) = pstate%wt_nphi_abs(:,iant)/zcpl_sum
     endif

     if(zpext.gt.ZERO) then
        pstate%picrf_ext(iant) = zpext
        pstate%wt_nphi_ext(:,iant) = pstate%wt_nphi_ext(:,iant)/zpext
     endif

  enddo

  if(iersum.gt.0) then
     write(lunmsg,*) ' ?rf_wavecode_output: errors occurred fetching'
     write(lunmsg,*) '  wavecode heating and current drive profiles.'
     ierr = 200
  else
     call ps_state_memory_update(ierr, state=pstate, eqcheck=.FALSE.)
  endif

  deallocate(zpzon,zwkbuf)

end subroutine rf_wavecode_output
