subroutine nubeam_ctrl_psinit(ss,psinit,ierr)

  !  set up initialization block based on plasma state contents.
  !  this is called ONCE ONLY during nubeam initialization sequence,
  !  when driven from plasma state in the component interface.

  use iso_c_binding, only: fp => c_double
  use physconst_mod, only: twopi
  use nbi_pstypes
  use plasma_state_mod

  implicit NONE

  !---------------------------------------------
  !  arguments:

  type (plasma_state) :: ss
  type (nbi_ps_init) :: psinit

  integer, intent(out) :: ierr

  !---------------------------------------------
  !  local

  integer :: ii,ij,ing,inx,inb,iz,ia,inum,jccw,iap2,iberr
  real(fp) :: zrtcen,zdivcon,zEmax
  character*32 :: zword
  character*120 zmsg
  !---------------------------------------------

  ierr = 0
  call nbi_ps_init_dflt(psinit)

  psinit%runid = ss%runid    ! runid

  psinit%nthets = ss%nth_eq  ! poloidal angle grid size (from equilibrium)

  jccw = ss%kccw_Jphi        ! +1 if J_phi flows CCW viewed from above
  !  ccw = "counter-clockwise"

  !  conversion factor for beam divergences: degrees -> radians
  !  (1/e)divergence / sqrt(2) to match NUBEAM definition

  zdivcon = twopi/(360.0_fp*sqrt(2.0_fp))

  !-----------------------------------------------------
  ! thermal plasma species: main gasses and impurities

  ing = 0
  inx = 0

  do ii=1,ss%nspec_th
     if(ss%s_type(ii).ne.ps_impurity) then
        ing = ing + 1
        psinit%backz(ing) = ss%q_s(ii)/ps_xe   ! Z of target ion species
        psinit%aplasm(ing) = ss%m_s(ii)/ps_mp  ! A of target ion species
     else
        inx = inx + 1
        psinit%xzimpx(inx) = ss%q_s(ii)/ps_xe  ! Z of impurity species
        psinit%aimpx(inx) = ss%m_s(ii)/ps_mp   ! A of impurity species
     endif
  enddo

  psinit%ngmax = ing
  psinit%nrhix = inx

  !-----------------------------------------------------
  ! RF minority species

  inx = ss%nspec_rfmin
  psinit%nmini = inx

  do ii=1,inx
     psinit%xzmini(ii)=ss%q_rfmin(ii)/ps_xe    ! Z of RF minority
     psinit%amini(ii)=ss%m_rfmin(ii)/ps_mp     ! A of RF minority
  enddo

  !-----------------------------------------------------
  ! Fusion product species

  psinit%nlfsp = inbi_logical(.FALSE.)
  psinit%nlfsp14 = inbi_logical(.FALSE.)
  psinit%nlfst = inbi_logical(.FALSE.)
  psinit%nlfhe3= inbi_logical(.FALSE.)
  psinit%nlfhe4= inbi_logical(.FALSE.)

  inx = ss%nspec_fusion

  do ii=1,inx
     iz = ss%qatom_sfus(ii)/ps_xe + 0.1_fp
     ia = ss%m_sfus(ii)/ps_mp + 0.1_fp
     if(iz.eq.1) then
        if(ia.eq.1) then
           if(ss%sfus_name(ii).eq."H_fusn") then
              psinit%nlfsp = inbi_logical(.TRUE.) ! proton
           elseif(ss%sfus_name(ii).eq."H14_fusn") then
              psinit%nlfsp14 = inbi_logical(.TRUE.) ! proton, 14MeV
           endif
        else if(ia.eq.3) then
           psinit%nlfst = inbi_logical(.TRUE.) ! triton
        endif
     else if(iz.eq.2) then
        if(ia.eq.3) then
           psinit%nlfhe3 = inbi_logical(.TRUE.) ! He3
        else if(ia.eq.4) then
           psinit%nlfhe4 = inbi_logical(.TRUE.) ! He4 (alpha)
        endif
     endif
  enddo

  !-----------------------------------------------------
  !  beams

  inb = ss%nbeam
  psinit%nbeam = inb

  !  species injected by each beam...

  iberr=0
  do ii=1,inb
     psinit%ntrace(ii)=0
     if(ss%nbion(ii).eq.'H') then
        psinit%xzbeama(ii)=1.0_fp
        psinit%abeama(ii)=1.0_fp

     else if(ss%nbion(ii).eq.'D') then
        psinit%xzbeama(ii)=1.0_fp
        psinit%abeama(ii)=2.0_fp

     else if(ss%nbion(ii).eq.'T') then
        psinit%xzbeama(ii)=1.0_fp
        psinit%abeama(ii)=3.0_fp

     else if((ss%nbion(ii).eq.'He3').or.(ss%nbion(ii).eq.'HE3')) then
        psinit%xzbeama(ii)=2.0_fp
        psinit%abeama(ii)=3.0_fp

     else if((ss%nbion(ii).eq.'He4').or.(ss%nbion(ii).eq.'HE4')) then
        psinit%xzbeama(ii)=2.0_fp
        psinit%abeama(ii)=4.0_fp

     else if((ss%nbion(ii).eq.'Ne').or.(ss%nbion(ii).eq.'NE')) then
        psinit%xzbeama(ii)=10.0_fp
        psinit%abeama(ii)=20.02_fp

     else if((ss%nbion(ii).eq.'Ar').or.(ss%nbion(ii).eq.'AR')) then
        psinit%xzbeama(ii)=18.0_fp
        psinit%abeama(ii)=39.635_fp

     else if((ss%nbion(ii).eq.'Kr').or.(ss%nbion(ii).eq.'KR')) then
        psinit%xzbeama(ii)=36.0_fp
        psinit%abeama(ii)=83.14_fp

     else if((ss%nbion(ii).eq.'Xe').or.(ss%nbion(ii).eq.'XE')) then
        psinit%xzbeama(ii)=54.0_fp
        psinit%abeama(ii)=130.27_fp

     else
        if(iberr.eq.0) iberr=ii
        ierr = ierr + 1
     endif
  enddo

  if(ierr.gt.0) then
     ierr=99
     call nbi_warning('?nubeam_ctrl_psinit: beam species in ss%nbion(:) not recognized:')
     call nbi_warning(' must be "H", "D", "T", "He3", "He4", or...')
     call nbi_warning(' a noble gas: "Ne", "Ar", "Kr", "Xe".')
     call nbi_warning(' instead the value received is: '//ss%nbion(iberr))
     return
  endif

  !  add in trace beams, if any...
  do ii=1,inb
     if(ss%nbion_trace(ii).eq.' ') cycle

     psinit%nbeam=psinit%nbeam + 1
     ij = psinit%nbeam
     psinit%ntrace(ij)=ii

     if(ss%nbion_trace(ii).eq.'H') then
        psinit%xzbeama(ij)=1.0_fp
        psinit%abeama(ij)=1.0_fp

     else if(ss%nbion_trace(ii).eq.'D') then
        psinit%xzbeama(ij)=1.0_fp
        psinit%abeama(ij)=2.0_fp

     else if(ss%nbion_trace(ii).eq.'T') then
        psinit%xzbeama(ij)=1.0_fp
        psinit%abeama(ij)=3.0_fp

     else if((ss%nbion_trace(ii).eq.'He3').or.(ss%nbion_trace(ii).eq.'HE3')) then
        psinit%xzbeama(ij)=2.0_fp
        psinit%abeama(ij)=3.0_fp

     else if((ss%nbion_trace(ii).eq.'He4').or.(ss%nbion_trace(ii).eq.'HE4')) then
        psinit%xzbeama(ij)=2.0_fp
        psinit%abeama(ij)=4.0_fp

     else if((ss%nbion_trace(ii).eq.'Ne').or.(ss%nbion_trace(ii).eq.'NE')) then
        psinit%xzbeama(ij)=10.0_fp
        psinit%abeama(ij)=20.02_fp

     else if((ss%nbion_trace(ii).eq.'Ar').or.(ss%nbion_trace(ii).eq.'AR')) then
        psinit%xzbeama(ij)=18.0_fp
        psinit%abeama(ij)=39.635_fp

     else if((ss%nbion_trace(ii).eq.'Kr').or.(ss%nbion_trace(ii).eq.'KR')) then
        psinit%xzbeama(ij)=36.0_fp
        psinit%abeama(ij)=83.14_fp

     else if((ss%nbion_trace(ii).eq.'Xe').or.(ss%nbion_trace(ii).eq.'XE')) then
        psinit%xzbeama(ij)=54.0_fp
        psinit%abeama(ij)=130.27_fp

     else
        if(iberr.eq.0) iberr=ii
        ierr = ierr + 1
     endif
  enddo

  if(ierr.gt.0) then
     ierr=99
     call nbi_warning('?nubeam_ctrl_psinit: beam species in ss%nbion_trace(:) not recognized:')
     call nbi_warning(' must be "H", "D", "T", "He3", "He4", or...')
     call nbi_warning(' a noble gas: "Ne", "Ar", "Kr", "Xe".')
     call nbi_warning(' instead the value received is: '//ss%nbion_trace(iberr))
     return
  endif

  !  beam geometry...

  do ii=1,inb
     zrtcen = 100.0_fp*ss%sRtcen(ii)  ! cm (signed: + means ccw momentum inj.)

     psinit%rtcena(ii)=abs(zrtcen)   ! tangency radius (absolute value)

     if(jccw*zrtcen.gt.0.0_fp) then
        psinit%nlco(ii) = inbi_logical(.TRUE.)    ! inj. w/ plasma current
     else
        psinit%nlco(ii) = inbi_logical(.FALSE.)   ! inj. against plasma current
     endif

     psinit%xlbtna(ii) = 100.0_fp*ss%Lbsctan(ii)   ! dist., sce to tangency pt.
     psinit%xybsca(ii) = 100.0_fp*ss%Zbsc(ii)      ! height, sce above midplane
     psinit%xbzeta(ii) = ss%Phibsc(ii)            ! toroidal angle of sce

     psinit%xlbapa(ii) = 100.0_fp*ss%Lbscap(ii)    ! dist., sce to aperture
     psinit%xybapa(ii) = 100.0_fp*ss%Zbap(ii)      ! height, aperture center

     zword = ss%nbshape(ii)                       ! shape of source
     call uupper(zword)  ! uppercase convert

     if(zword.eq.'CIRCLE') then
        ! circular ion source grid
        psinit%nbshapa(ii) = 2
        psinit%bmwidra(ii) = 100.0_fp*ss%b_halfwidth(ii)     ! radius
        psinit%foclra(ii) = 100.0_fp*ss%b_Hfocal_length(ii)  ! focal length
        psinit%divra(ii) = zdivcon*ss%b_Hdivergence(ii)     ! divergence
     else if(zword.eq.'RECTANGLE') then
        ! rectangular ion source grid
        psinit%nbshapa(ii) = 1
        psinit%bmwidra(ii) = 100.0_fp*ss%b_halfwidth(ii)     ! half-width
        psinit%bmwidza(ii) = 100.0_fp*ss%b_halfheight(ii)    ! half-height
        psinit%foclra(ii) = 100.0_fp*ss%b_Hfocal_length(ii)  ! horiz. focal len
        psinit%foclza(ii) = 100.0_fp*ss%b_Vfocal_length(ii)  ! vert. focal len
        psinit%divra(ii) = zdivcon*ss%b_Hdivergence(ii)     ! horiz. diverg.
        psinit%divza(ii) = zdivcon*ss%b_Vdivergence(ii)     ! vert.  diverg.
     else
        call nbi_warning('?nubeam_ctrl_psinit: beam shape in state (ss%nbshape(:) unrecognized: '//trim(zword))
        ierr = ierr + 1
     endif

     zword = ss%nbap_shape(ii)  ! aperture shape
     call uupper(zword)  ! uppercase convert

     if(zword.eq.'CIRCLE') then
        psinit%nbapsha(ii) = 2
        psinit%rapedga(ii) = 100.0_fp*ss%ap_halfwidth(ii)
     else if(zword.eq.'RECTANGLE') then
        psinit%nbapsha(ii) = 1
        psinit%rapedga(ii) = 100.0_fp*ss%ap_halfwidth(ii)
        psinit%xzpedga(ii) = 100.0_fp*ss%ap_halfheight(ii)
     else
        call nbi_warning('?nubeam_ctrl_psinit: aperture shape in state (ss%nbap_shape(:) unrecognized: '//trim(zword))
        ierr = ierr + 1
     endif

     psinit%xrapoffa(ii) = 100.0_fp*ss%ap_horiz_offset(ii)
     psinit%xzapoffa(ii) = 100.0_fp*ss%ap_vert_offset(ii)

     zword = ss%nbap2_shape(ii)  ! 2nd aperture shape
     call uupper(zword)  ! uppercase convert

     iap2=0
     if((zword.eq.' ').or.(zword.eq.'NONE')) then
        continue         ! no 2nd aperture
     else if(zword.eq.'CIRCLE') then
        iap2=1
        psinit%nbapsh2(ii) = 2
        psinit%xlbapa2(ii) = 100.0_fp*ss%Lbscap2(ii)
        psinit%rapedg2(ii) = 100.0_fp*ss%ap2_halfwidth(ii)
     else if(zword.eq.'RECTANGLE') then
        iap2=1
        psinit%nbapsh2(ii) = 1
        psinit%xlbapa2(ii) = 100.0_fp*ss%Lbscap2(ii)
        psinit%rapedg2(ii) = 100.0_fp*ss%ap2_halfwidth(ii)
        psinit%xzpedg2(ii) = 100.0_fp*ss%ap2_halfheight(ii)
     else
        call nbi_warning('?nubeam_ctrl_psinit: aperture shape in state (ss%nbap_shape(:) unrecognized: '//trim(zword))
        ierr = ierr + 1
     endif

     if(iap2.eq.1) then
        psinit%xrapoff2(ii) = 100.0_fp*ss%ap2_horiz_offset(ii)
        psinit%xzapoff2(ii) = 100.0_fp*ss%ap2_vert_offset(ii)
     endif

  enddo
  if(ierr.ne.0) return

  zEmax = 0.0_fp
  do ii=1,inb
     zEmax = max(zEmax, 1000.0_fp*ss%Einj_max(ii))   ! max Einj of beam
     psinit%einjmin(ii) = 1000.0_fp*ss%Einj_min(ii)  ! min Einj of beam

     ! place holders for beam power & voltage
     psinit%pinja(ii) = 0.0_fp
     psinit%einja(ii) = psinit%einjmin(ii)
     psinit%ffulla(ii) = 1.0_fp
     psinit%fhalfa(ii) = 0.0_fp
  enddo

  psinit%ebdmax = 1.5_fp * zEmax   ! set max(E) of beam ion dist.fcns.

  !  use Psi(R,Z) (MHD equilibrium grid sizes) here...
!MG 11/16 psinit%orbrzv_nr = ss%nR
!MG 11/16 psinit%orbrzv_nz = ss%nZ

  !  atomic physics settings

  ing = ss%ngsc0
  psinit%nsmax = ing + 1    ! no. of neutral sources (add 1 for halos)
  psinit%khalo(1:ing) = 0
  psinit%khalo(ing+1) = 1

  psinit%nmsigx = 2         ! Phaneuf set impurity stopping cross sections

  !  TF ripple

  psinit%ncoils = ss%ntf_coils
  psinit%ncoils_2 = ss%ntf_coil2

  !  anomalous fast ion transport

  psinit%x_qlim0 = 0.0_fp

  psinit%nkdifb = 3         ! applicable to all species (separate profiles
  !  of v & D provided for each class of species).

  inum = ss%nefi_anom
  if(inum.gt.size(psinit%edifbe)) then

     zmsg=' '
     write(zmsg,*) ' ss%nefi_anom = ',ss%nefi_anom,'; NUBEAM max ndifbe: ', &
          size(psinit%edifbe)
     call nbi_warning('ss%nefi_anom > NUBEAM internal limit, amomalous diffusivity energy variation.')
     ierr=1

  else

     psinit%ndifbe = inum
     if(inum.gt.0) then
        psinit%edifbe(1:inum) = 1000.0_fp*ss%E_anom(1:inum)
     endif

  endif

  ! new DMC Mar 2010:  orbit topology dependent diffusivities vs. (x,E)

  ! set the type of fast ion diffusivities (xyuan, Mar, 2012)
  psinit%nmdifb = ss%nmdifb

  if(ss%nefi_anom2.gt.0) then
     psinit%ndifbe_2d = ss%nefi_anom2
     psinit%dif2d_conv = 1.0e+4_fp  ! m**2/s -> cm**2/s
  endif

end subroutine nubeam_ctrl_psinit
