subroutine rfxqlo_state_input(pstate,imsym,ierr)

  ! initialize rfxqlo from Plasma State -- to the extent possible
  ! RFXQLO module is cleared and default values are set; then modified
  ! from state data.

  ! DMC May 2008 -- information settable from Plasma State is expected
  ! to evolve...

  ! DMC Jul 2011 -- add code to screen out species with Z & A values that
  ! conflict with RF minority species; at most <max_screen> may be so
  ! screened, and the most populous thermal specie is not eligible, nor 
  ! are impurities...

  use plasma_state_mod
  use rfxqlo_interface
  implicit NONE

  !-----------------------------
  !  passed:
  type (plasma_state) :: pstate

  logical, intent(in) :: imsym(*)  ! nphi symmetrization option
  !  (imsym: LOGICAL array of size #antennas, at least)
  !  imsym(j)=T means: treat antenna #j nphi spectrum as symmetric
  !                    pstate nphi and weighting must be consistent with this

  integer, intent(out) :: ierr

  integer :: ii,iirf,iisav,jj,inrho,inth,itype,iersum,ierloc,inum,inum2,inzon
  integer :: id_rho,id_fun
  logical :: iok
  real*8, dimension(:,:), allocatable :: zna,zta
  real*8, dimension(:), allocatable :: zna_bdy
  real*8, dimension(:),  allocatable :: zwk
  character*32 znam,znamu,znamu_save
  !-----------------------------
  !  local:

  integer :: lunmsg,iwarn

  type (rfxqlo_dims) :: zdims
  type (rfxqlo_vacuum) :: zvac
  type (rfxqlo_system) :: zsys
  type (rfxqlo_profile) :: zprof
  type (rfxqlo_antenna) :: zant
  type (rfxqlo_ion) :: zions
  type (rfxqlo_rfgrid) :: zgrid
  type (rfxqlo_options) :: zopt

  real*8 :: zg,zq,zraxis,zflux,zdvdrho,zmu0
  real*8, parameter :: cpi  = 3.1415926535897931D+00
  real*8, parameter :: c2pi = 6.2831853071795862D+00
  real*8 :: zedge 
  real*8 :: zedge_ne ! use this factor to prevent quasineutrality violation

  integer, parameter :: max_screen = 2
  integer :: num_screen,imajority,intot,indx,ispec

  logical, dimension(:), allocatable :: lscreen
  integer, dimension(:), allocatable :: iza,iaa
  real*8, dimension(:), allocatable :: zntot,zdvol,znwk,znesum
  real*8 :: znmax,znesum_bdy,znmini_bdy
  character*1 msgchar
  logical ::  nl_first

  !-----------------------------
  zedge=0.0_rspec
  iersum = 0
  call rfxqlo_get_lunmsg(lunmsg)

  if(pstate%nrho.eq.0) then
     iersum = iersum + 1
     write(lunmsg,*) ' ?rfxqlo_state_input: pstate%nrho = 0.'
  endif

  if(pstate%nrho_eq.eq.0) then
     iersum = iersum + 1
     write(lunmsg,*) ' ?rfxqlo_state_input: pstate%nrho_eq = 0.'
  endif

  if(iersum.gt.0) then
     ierr=1
     return
  endif

  !-----------------------------
  !  check / set up rfxqlo module dimensioning:
  !  dimensions: no. of antenna, max no. of Nphi/antenna, ...

  call rfxqlo_dim_get(zdims,iwarn)

  if(iwarn.ne.0) then
     ! set the dimensions now; they were not set previously...

     call rfxqlo_dim_dflts(zdims)

     zdims%n_antenna = pstate%nicrf_src
     zdims%max_nphi = pstate%max_num_nphi

     call rfxqlo_dim_set(zdims,ierloc)
     iersum = iersum + abs(ierloc)

  else
     ! check that dimensions match requirements of state data

     if(zdims%n_antenna.ne.pstate%nicrf_src) then
        write(lunmsg,*) ' ?rfxqlo_state_input: in rfxqlo, preset n_antenna =',&
             zdims%n_antenna
        write(lunmsg,*) '  does not match nicrf_src in plasma state data: ', &
             pstate%nicrf_src
        iersum = iersum + 1
     endif

     if(zdims%max_nphi.lt.pstate%max_num_nphi) then
        write(lunmsg,*) ' ?rfxqlo_state_input: in rfxqlo, preset max_nphi =',&
             zdims%max_nphi
        write(lunmsg,*) '  does not match maximum in plasma state data: ', &
             pstate%max_num_nphi
        iersum = iersum + 1
     endif
  endif

  if(iersum.ne.0) then
     ierr=1
     return
  endif

  !-----------------------------
  call rfxqlo_vacuum_dflts(zvac)

  zg = pstate%g_eq(pstate%nrho_eq)  ! vacuum R*B_phi
  zq = pstate%q_eq(pstate%nrho_eq)  ! q at edge
  zraxis = pstate%R_axis

  zvac%rmaj0 = zraxis
  zvac%btor0 = zg/zraxis
  zflux = pstate%phit(pstate%nrho_eq)
  zvac%tflux = zflux              ! enclosed toroidal flux

  call ps_intrp_1d_sc(1.0_rspec,pstate%id_vol,zdvdrho,ierloc, &
       ideriv=1, state=pstate)
  iersum = iersum + abs(ierloc)

  ! Ampere's law estimate of enclosed current

  zvac%placur = pstate%curt(pstate%nrho_eq)

  call rfxqlo_vacuum_set(zvac,ierloc)
  iersum = iersum + abs(ierloc)

  !-----------------------------
  call rfxqlo_system_dflts(zsys)

  zsys%time = pstate%t0

  zsys%run_id = pstate%runid
  zsys%run_label = pstate%runid//' loaded from state by rfxqlo_state_input.f90'

  call rfxqlo_system_set(zsys,ierloc)
  iersum = iersum + abs(ierloc)

  !-----------------------------
  call rfxqlo_profile_dflts(zprof)

  ! "zero" profile
  inzon = pstate%nrho
  allocate(zwk(inzon)); zwk = 0.0_rspec

  zedge=0.0_rspec
  call mkpro('rho',zwk,'_RFXQLO_ZERO',zedge)
  deallocate(zwk)

  zprof%xplasma_file = '(plasma state)'

  zprof%xnam_g = 'G_EQ'
  zprof%xnam_q = 'Q_EQ'

  zprof%xnam_te = 'RFX_TE'
  zedge=0.0_rspec
  call mkpro('rho',pstate%ts(:,0),zprof%xnam_te,zedge)

  zprof%xnam_ne = 'RFX_NE'
  zedge_ne=0.0_rspec
  call mkpro('rho',pstate%ns(:,0),zprof%xnam_ne,zedge_ne,pstate%ns_bdy(0))

  zprof%xnam_Tiav = 'RFX_TIAV'
  zedge=0.0_rspec
  call mkpro('rho',pstate%ti(:),zprof%xnam_tiav,zedge)

  zprof%xnam_vloop = 'V_loop'  ! use state profile "as is"

  zprof%xnam_Epot = 'Epot'     ! use state profile "as is"

  zprof%xnam_omega = 'RFX_OMEGA'
  zedge=0.0_rspec
  call mkpro('rho',pstate%omegat(:),zprof%xnam_omega,zedge)

  call rfxqlo_profile_set(zprof,ierloc)
  iersum = iersum + abs(ierloc)

  !-----------------------------
  do ii=1,pstate%nicrf_src
     call rfxqlo_select_antenna(ii,ierloc)
     if(ierloc.ne.0) exit

     call rfxqlo_antenna_dflts(zant)

     zant%freq = pstate%freq_ic(ii)  ! Hz
     zant%power = pstate%power_ic(ii)

     zant%dx_farsh = pstate%dx_fshield(ii)

     inum = pstate%nrz_antGeo(ii)
     zant%nantrz = inum
     if(inum.le.0) then
        iersum = iersum + 1
        write(lunmsg,*) ' ?rfxqlo_state_input: antenna# ',ii, &
             ' (R,Z) geometry has 0 points.'
        
     else
        zant%antRgeo(1:inum) = pstate%R_antgeo(1:inum,ii)
        zant%antZgeo(1:inum) = pstate%Z_antgeo(1:inum,ii)
     endif

     zant%sym_nphi = imsym(ii)
     inum = pstate%num_nphi(ii)
     if(inum.le.0) then
        iersum = iersum + 1
        write(lunmsg,*) ' ?rfxqlo_state_input: antenna# ',ii, &
             ' n_phi spectrum has 0 points.'

     else if(zant%sym_nphi) then
        call cksym(pstate%nphi(1:inum,ii),pstate%wt_nphi(1:inum,ii),iok)
        if(iok) then
           inum2=inum/2
           zant%num_nphi = inum2
           zant%nphi(1:inum2) = pstate%nphi(inum2+1:inum2+inum2,ii)
           zant%wnphi_vac(1:inum2) = 2*pstate%wt_nphi(inum2+1:inum2+inum2,ii)

        else
           iersum = iersum + 1
           write(lunmsg,*) ' ?rfxqlo_state_input: on antenna #',ii,':'
           write(lunmsg,*) '  symmetric n_phi spectrum expected but not found:'
           write(lunmsg,*) '    n_phi = ',pstate%nphi(1:inum,ii)
           write(lunmsg,*) '    wt_n_phi = ',pstate%wt_nphi(1:inum,ii)
        endif

     else
        zant%num_nphi = inum
        zant%nphi(1:inum) = pstate%nphi(1:inum,ii)
        zant%wnphi_vac(1:inum) = pstate%wt_nphi(1:inum,ii)
     endif

     call rfxqlo_antenna_set(zant,ierloc)
     if(ierloc.ne.0) exit
  enddo
  iersum = iersum + abs(ierloc)

  !-----------------------------
  ! deal with ions.  Allow "screening" of ions with (Z,A) values that
  ! conflict with RF minority ions -- per request of physics user groups...

  call rfxqlo_ion_dflts(zions)

  inrho = pstate%nrho - 1    ! #zones
  inth = pstate%nspec_tha    ! #thermal (nothing screened)
  intot = pstate%nspec_alla  ! #total (nothing screened)

  !  arrays...

  allocate(zna(inrho,inth),zta(inrho,inth),zdvol(inrho),zwk(inrho+1))
  allocate(zna_bdy(inth))
  allocate(lscreen(intot),iza(intot),iaa(intot),zntot(inth))
  allocate(znwk(inrho),znesum(inrho))

  !  zone volumes...
  call ps_intrp_1d(pstate%rho,pstate%id_vol,zwk,ierloc,state=pstate)
  if(ierloc.ne.0) then
     iersum=iersum+1
     zdvol=1.0_rspec
  else
     zdvol(1:inrho) = (zwk(2:inrho+1)-zwk(1:inrho))
  endif
  deallocate(zwk)

  !  get densities & temperatures...
  call ps_tha_fetch(ierloc,state=pstate,ns=zna,ts=zta)
  zna_bdy=0.0_rspec
  do ispec=1,pstate%nspec_th
     indx=pstate%s_to_sa(ispec)
     if(indx.gt.0) then
        zna_bdy(indx) = pstate%ns_bdy(ispec)
     endif
  enddo

  num_screen = 0

  znmax=0.0_rspec
  zntot=0.0_rspec
  if(ierloc.ne.0) then
     iersum=iersum+1
     zna=0.0_rspec
     zta=0.0_rspec
     imajority=1
  else
     ! set Z & A arrays;
     do ii=1,intot
        itype = pstate%alla_type(ii)

        lscreen(ii)=.FALSE.
        iza(ii) = pstate%q_alla(ii)/ps_xe + 0.1_rspec
        iaa(ii) = pstate%m_alla(ii)/ps_mp + 0.1_rspec

        if(.NOT.((itype.eq.ps_impurity).OR.(itype.eq.ps_therm_ion))) cycle

        ! compute Ntot of each thermal specie
        do jj=1,inrho
           zntot(ii)=zntot(ii)+zdvol(jj)*zna(jj,ii)
        enddo
        if(itype.eq.ps_impurity) cycle

        ! find largest Ntot of any non-impurity specie
        if(zntot(ii).gt.znmax) then
           znmax=zntot(ii)
           imajority=ii
        endif
     enddo
  endif

  ! find screened species, if any...
  znesum = 0.0_rspec
  znesum_bdy = 0.0_rspec
  goto 111 !MG 08/2019 ignore the following part

  nl_first=.TRUE.
  do iirf=1,intot
     if(pstate%alla_type(iirf).ne.ps_rf_minority) cycle

     do ii=1,intot
        if(ii.eq.iirf) cycle
        if(iza(ii).ne.iza(iirf)) cycle
        if(iaa(ii).ne.iaa(iirf)) cycle

        ! Z & A match

        itype = pstate%alla_type(ii)
        if(itype.eq.ps_impurity) then
           iersum = iersum + 1
           write(lunmsg,*) &
                ' ?rfxqlo_state_input: RF minority Z & A indistinguishable'
           write(lunmsg,*) &
                '  from impurity: '//trim(pstate%alla_name(iirf))//'; '// &
                trim(pstate%alla_name(ii))
           cycle
        endif
        if(itype.eq.ps_rf_minority) then
           iersum = iersum + 1
           write(lunmsg,*) &
                ' ?rfxqlo_state_input: two RF minorities Z & A values'
           write(lunmsg,*) &
                '  indistinguishable: '//trim(pstate%alla_name(iirf))//'; '// &
                trim(pstate%alla_name(ii))
           cycle
        endif
        if(ii.eq.imajority) then
           iersum = iersum + 1
           write(lunmsg,*) &
                ' ?rfxqlo_state_input: RF minority Z & A indistinguishable'
           write(lunmsg,*) &
                '  from thermal specie: '//trim(pstate%alla_name(iirf))// &
                '; '//trim(pstate%alla_name(ii))
           write(lunmsg,*) '  ...this is the majority thermal specie.'
           cycle
        endif

        if(num_screen.lt.max_screen) then
           msgchar = '%'
           num_screen = num_screen + 1
           lscreen(ii)=.TRUE.
           if(itype.eq.ps_beam_ion.and.nl_first) then
              num_screen = num_screen - 1
              lscreen(ii)=.FALSE.  !one beam ion sp is allowed
              nl_first=.FALSE.
           endif
        else
           msgchar = '?'
           iersum = iersum + 1
           !  lscreen(ii)=.FALSE.  (already set)
        endif

        write(lunmsg,*) ' '//msgchar// &
             'rfxqlo_state_input: RF minority Z & A indistinguishable'

        jj=pstate%alla_index(ii)
        if(itype.eq.ps_therm_ion) then
           write(lunmsg,*) &
                '  from thermal specie: '//trim(pstate%alla_name(iirf))// &
                '; '//trim(pstate%alla_name(ii))
           id_fun=0
        else if(itype.eq.ps_beam_ion) then
           write(lunmsg,*) &
                '  from beam specie: '//trim(pstate%alla_name(iirf))//'; '// &
                trim(pstate%alla_name(ii))
           if(allocated(pstate%id_nbeami)) then
              id_fun = pstate%id_nbeami(jj)
           else
              id_fun=-111 !there is no beam fast ion source yet
              !recover screening result
              if(.not.lscreen(ii)) then
                 num_screen = num_screen + 1
                 lscreen(ii)=.TRUE.
                 nl_first=.TRUE.
              endif
           endif
        else if(itype.eq.ps_fusion_ion) then
           write(lunmsg,*) &
                '  from fusion product: '//trim(pstate%alla_name(iirf))// &
                '; '//trim(pstate%alla_name(ii))
           if(allocated(pstate%id_nfusi)) then
              id_fun = pstate%id_nfusi(jj)
           else
              id_fun=-111 !there is no fusion source yet
           endif
        endif
        if(.NOT.lscreen(ii)) write(lunmsg,*) ' ...too many Z & A matches...'

        if(id_fun.eq.0) then
           znesum = znesum + iza(ii)*zna(:,ii)
           znesum_bdy = znesum_bdy + iza(ii)*zna_bdy(ii)
        elseif(id_fun.gt.0) then
           call ps_rho_rezone(pstate%rho,id_fun,znwk,ierloc, state=pstate)
           if(ierloc.ne.0) then
              iersum = iersum + 1
           else
              znesum = znesum + iza(ii)*znwk(:)
           endif
        endif
     enddo
  enddo
111 continue  !MG 08/2019
  imajority=0 !MG 08/2019
  zions%n_ions = intot - num_screen

  iisav = 0
  do ii=1,intot
     if(lscreen(ii)) cycle
     iisav = iisav + 1

     znam = pstate%alla_name(ii)
     zions%name_ions(iisav) = znam
     zions%a_ions(iisav) = pstate%m_alla(ii)/ps_mp
     zions%z_ions(iisav) = pstate%q_alla(ii)/ps_xe

     itype = pstate%alla_type(ii)
     if((itype.eq.ps_therm_ion).or.(itype.eq.ps_impurity)) then
        znamu = 'RFX_'//trim(znam)//'_N'
        zions%xnam_ni(iisav)=znamu
        if(ii.eq.imajority) then
           call mkpro('rho',zna(:,ii) + znesum(:)/iza(ii),znamu,zedge_ne,zna_bdy(ii) + znesum_bdy/iza(ii))
        else
           call mkpro('rho',zna(:,ii),znamu,zedge_ne,zna_bdy(ii))
        endif

        znamu = 'RFX_'//trim(znam)//'_T'
        zions%xnam_ti(iisav)=znamu
        zions%xnam_tpll(iisav)=znamu
        zedge=0.0_rspec
        call mkpro('rho',zta(:,ii),znamu,zedge)
        znamu_save = znamu
     else
        jj=pstate%alla_index(ii)
        if(itype.eq.ps_beam_ion) then
           if(pstate%nrho_nbi.eq.0) then
              zions%xnam_ni(iisav)='_RFXQLO_ZERO'
              zions%xnam_ti(iisav)=znamu_save
              zions%xnam_tpll(iisav)=znamu_save
           else
              allocate(zwk(pstate%nrho_nbi-1))

              znamu = 'RFX_'//trim(znam)//'_N'
              zions%xnam_ni(iisav)=znamu
              call mkpro('rho_nbi',pstate%nbeami(:,jj),znamu,zedge_ne,pstate%nbeami_bdy(jj))

              zwk = (2.0_rspec/3.0_rspec)* &
                   (pstate%eperp_beami(:,jj) + pstate%epll_beami(:,jj))
              znamu = 'RFX_'//trim(znam)//'_T'
              zions%xnam_ti(iisav)=znamu
              zedge=0.0_rspec
              call mkpro('rho_nbi',zwk,znamu,zedge)

              zwk = (2.0_rspec)*pstate%epll_beami(:,jj)
              znamu = 'RFX_'//trim(znam)//'_TPLL'
              zions%xnam_tpll(iisav)=znamu
              zedge=0.0_rspec
              call mkpro('rho_nbi',zwk,znamu,zedge)

              deallocate(zwk)
           endif

        else if(itype.eq.ps_fusion_ion) then
           if(pstate%nrho_fus.eq.0) then
              zions%xnam_ni(iisav)='_RFXQLO_ZERO'
              zions%xnam_ti(iisav)=znamu_save
              zions%xnam_tpll(iisav)=znamu_save
           else
              allocate(zwk(pstate%nrho_fus-1))

              znamu = 'RFX_'//trim(znam)//'_N'
              zions%xnam_ni(iisav)=znamu
              call mkpro('rho_fus',pstate%nfusi(:,jj),znamu,zedge_ne,pstate%nfusi_bdy(jj))

              zwk = (2.0_rspec/3.0_rspec)* &
                   (pstate%eperp_fusi(:,jj) + pstate%epll_fusi(:,jj))
              znamu = 'RFX_'//trim(znam)//'_T'
              zions%xnam_ti(iisav)=znamu
              zedge=0.0_rspec
              call mkpro('rho_fus',zwk,znamu,zedge)

              zwk = (2.0_rspec)*pstate%epll_fusi(:,jj)
              znamu = 'RFX_'//trim(znam)//'_TPLL'
              zions%xnam_tpll(iisav)=znamu
              zedge=0.0_rspec
              call mkpro('rho_fus',zwk,znamu,zedge)

              deallocate(zwk)
           endif

        else if(itype.eq.ps_rf_minority) then
           if(pstate%nrho_icrf.eq.0) then
              zions%xnam_ni(iisav)='_RFXQLO_ZERO'
              zions%xnam_ti(iisav)=znamu_save
              zions%xnam_tpll(iisav)=znamu_save
           else
              allocate(zwk(pstate%nrho_icrf-1))

              znamu = 'RFX_'//trim(znam)//'_N'
              zions%xnam_ni(iisav)=znamu
              if(allocated(pstate%nmini_bdy)) then
                znmini_bdy=pstate%nmini_bdy(jj)
              else
                znmini_bdy=0.0d0
              endif
              call mkpro('rho_icrf',pstate%nmini(:,jj),znamu,zedge_ne,znmini_bdy)

              zwk = (2.0_rspec/3.0_rspec)* &
                   (pstate%eperp_mini(:,jj) + pstate%epll_mini(:,jj))
              znamu = 'RFX_'//trim(znam)//'_T'
              zions%xnam_ti(iisav)=znamu
              zedge=0.0_rspec
              call mkpro('rho_icrf',zwk,znamu,zedge)

              zwk = (2.0_rspec)*pstate%epll_mini(:,jj)
              znamu = 'RFX_'//trim(znam)//'_TPLL'
              zions%xnam_tpll(iisav)=znamu
              zedge=0.0_rspec
              call mkpro('rho_icrf',zwk,znamu,zedge)

              deallocate(zwk)
           endif

        endif
     endif
  enddo

  call rfxqlo_ion_set(zions,ierloc)
  iersum = iersum + abs(ierloc)

  deallocate(zna,zta,zdvol,zntot,znwk,znesum)
  deallocate(lscreen,iza,iaa)
  !-----------------------------
  call rfxqlo_rfgrid_dflts(zgrid)

  call eq_ganum('RHO_ICRF',id_rho)
  zgrid%xnam_xqlo="RHO_ICRF"
  zgrid%id_xqlo=id_rho

  call rfxqlo_rfgrid_set(zgrid,ierloc)
  iersum = iersum + abs(ierloc)

  !-----------------------------
  call rfxqlo_options_dflts(zopt)
  call rfxqlo_options_set(zopt,ierloc)
  iersum = iersum + abs(ierloc)

  !-----------------------------
  if(iersum.gt.0) then
     ierr=1
  else
     ierr=0
  endif

CONTAINS

  subroutine mkpro(rhonam,frhozc,znam,edge_ratio,frhozb)

    !  make a Hermite profile -- 0 slope on axis
    !     linear extrap to edge but avoid sign change
    !
    ! factor 'edge_ratio' can be applied at the boundary condition
    ! frhob(isize+1) to save quasineutrality
    !
    ! for example: aplly 'edge_ratio' computed for electron density
    !              to all ion densities  

    character*(*), intent(in) :: rhonam
    real*8, dimension(:), intent(in) :: frhozc
    character*(*), intent(in) :: znam
    real*8, intent(inout) :: edge_ratio ! to k
     real*8, intent(in), optional :: frhozb !value at zone boundary
   !---------------------
    integer :: isize,id_rho,id,iii
    real*8, dimension(:), allocatable :: frhob
    real*8, parameter :: zbc=0.0_rspec
    !---------------------

    isize=size(frhozc)
    allocate(frhob(isize+1))

    frhob(1)=(9.0_rspec*frhozc(1) - frhozc(2))/8.0_rspec
    if(present(frhozb)) then
       frhob(isize+1)=frhozb
       edge_ratio=0.0_rspec
    else
       frhob(isize+1)=(3.0_rspec*frhozc(isize) - frhozc(isize-1))/2.0_rspec ! linear interpolation
       if(frhozc(isize).ge.0.0_rspec) then
          frhob(isize+1)=max(0.1_rspec*frhozc(isize),frhob(isize+1))
       else
          frhob(isize+1)=min(0.1_rspec*frhozc(isize),frhob(isize+1))
       endif
    endif


    do iii=2,isize
       frhob(iii)=0.5_rspec*(frhozc(iii-1)+frhozc(iii))
    enddo
!
! calculate or apply 'edge_ratio'
!
    if(edge_ratio.gt.0.0_rspec) then
      frhob(isize+1)=edge_ratio* frhob(isize)
    else
       if(frhob(isize).gt.0.0_rspec) then
          edge_ratio=frhob(isize+1)/frhob(isize)
       else
          edge_ratio=0.0_rspec
       endif
    endif

    call eq_ganum(rhonam,id_rho)

    call eqm_rhofun(101,id_rho,znam,frhob,1,zbc,0,zbc,id,ierloc)

    iersum=iersum+abs(ierloc)

  end subroutine mkpro

  subroutine cksym(nphia,wt_nphia,iok)

    integer, dimension(:), intent(in) :: nphia
    real*8, dimension(:), intent(in) :: wt_nphia

    logical, intent(out) :: iok

    !-------------------------------
    !  check if nphi spectrum is symmetric
    !  if so return iok=.TRUE. else return iok=.FALSE.
    !-------------------------------

    integer :: inum,inum2,ii
    real*8 :: wdiff
    real*8, parameter :: wtol=1.0d-5

    inum = size(nphia)
    inum2 = inum/2
    if(inum2*2.ne.inum) then
       write(lunmsg,*) ' spectrum of odd length ',inum,' cannot be symmetric.'
       iok=.FALSE.
       return
    endif

    do ii=1,inum2
       if(nphia(inum2+ii).ne.-nphia(inum2+1-ii)) then
          iok=.FALSE.
          return
       endif
       
       wdiff = abs(wt_nphia(inum2+ii)-wt_nphia(inum2+1-ii))
       if(wdiff.gt.wtol) then
          iok=.FALSE.
          return
       endif
    enddo

    iok=.TRUE.
  end subroutine cksym

end subroutine rfxqlo_state_input
