!********************************************************
!	this subroutine sets up fixed beam and plasma parameters

subroutine parms(ierr)
  !  tbt jun 1993 - thermal species upgrade - from 2 to 5.
  !
  !  dmc mar 1985 - geometry generalization
  !   ** afac kluge deleted ** beam-neutral interactions will take
  !   place on the 2d beam dist fcn grid
  !      spatially 2d volume source thermal neutral density profile
  !   evaluated
  !
  !  dmc feb 1987 - neutral angular velocity profiles processed
  !
  !  dmc mar 2008 - revised handling of halo neutral sources; khalo option added.
  !
  !	neutral beam internal common blocks
  !
  !x	include 'nbcom.for/nolist'
  use nbi_com
  use orbrzv, only :orbrzv_bnw_lrf,nlmcrfgcon,mcrf_on,orbrzv_rfdv,&
       orbrzv_out,rfdbg_pdiag
  implicit none
  integer :: ig,j,is,ivf,inumi,jz2,iz2,izg,jz1,iz1,ien,ib
  integer :: isb,iz,iopstat,iopstat1,idi,jmid
  !============
  real(fp) :: zn00,ze0,zom0,zvol,zsrbn0,zn0sum
  !============
  real(fp) :: zn0_inner,zn0_outer
  !
  real(fp) :: zn0(mig),zn0w(mig),zp0(mig),zomg0(mig),zstot(mig)&
       ,ze_rel,zrms_rel,zbnb
  character(len=120) :: zfile
  !
  integer :: ihalo(mis)
  !
  !-------------------------------
  !
  integer :: iersum,ierr,iz2_min,ino
  integer, save :: imcrf_first=0
  !
  !-------------------------------
  !     run parameters
  !
  !     set up
  !
  !%	nbplot=1
  !%	ndep0=1
  !
  !  initialize the field description for the beam code ***
  !   part 1 (part 2 requires binv calculation, below)
  !
  iersum=0

  call bparms(1,ierr)
  iersum=iersum+abs(ierr)
  !
  !
  !  thermal neutral densities *****************
  !
  !    vol sce neutral densities are summed in the old 1d arrays
  !    but will be mapped to the 2d arrays proportional to source fcn
  !
  !  first: determine which sources are vol sce or beam halo sources...

  ihalo = khalo   ! (1:mis)
  do is=1,mis
     if((ihalo(is).lt.-1).or.(ihalo(is).gt.1)) then
        write(nonlin,*) &
             ' ?illegal halo flag value for neutral source #',is
        write(nonlin,*) '  expected -1, 0, or 1.'
        iersum=iersum+1
        ihalo=0
        exit
     end if

     if(ihalo(is).eq.-1) then
        ! auto-detect halo or edge source
        zn0_inner=0.0_fp
        zn0_outer=0.0_fp
        jmid=(lcentr+ledge)/2
        do j=lcentr,ledge
           zn0sum=0.0_fp
           do ig=1,ng
              zn0sum=zn0sum + den0mn(ig,j,is)
           end do
           if(j.le.jmid) then
              zn0_inner=max(zn0sum,zn0_inner)
           else
              zn0_outer=max(zn0sum,zn0_outer)
           end if
        end do

        if(zn0_inner.gt.zn0_outer/3) then
           ihalo(is)=1
        else
           ihalo(is)=0
        end if
     end if
  end do

  !  number of neutral sources-- mis
  !
  !  clear and define volume and wall neut. densities and
  !  avg energies
  do j=1,mzones
     do ig=1,ng             ! tbt
        voln0(ig,j)=0.0_fp
        walln0(ig,j)=0.0_fp
        evol0(ig,j)=0.0_fp
        ewall0(ig,j)=0.0_fp
        ovol0(ig,j)=0.0_fp
        owall0(ig,j)=0.0_fp
     end do
  end do
  !
  do is=1,mis
     ivf=ihalo(is)
     do j=lcentr,ledge
        do ig=1,ng
           zn00=den0mn(ig,j,is)
           ze0=en0mn(ig,j,is)
           zom0=omg0mn(ig,j,is)

           if(ivf.eq.0) then
              !  wall (edge source) neutrals modeled as cylindrically symmetric
              walln0(ig,j)=walln0(ig,j)+zn00
              ewall0(ig,j)=ewall0(ig,j)+zn00*ze0
              owall0(ig,j)=owall0(ig,j)+zn00*zom0
           else
              !  volume sce neutrals modeled as in vicinity of beam
              voln0(ig,j)=voln0(ig,j)+zn00
              evol0(ig,j)=evol0(ig,j)+zn00*ze0
              ovol0(ig,j)=ovol0(ig,j)+zn00*zom0
           end if

        end do
     end do
  end do

  !  divide out density to get density-weighted avg
  !  energy of neutrals
  do j=lcentr,ledge
     do ig=1,ng
        if(walln0(ig,j).gt.0._fp) then
           ewall0(ig,j)=ewall0(ig,j)/walln0(ig,j)
           owall0(ig,j)=owall0(ig,j)/walln0(ig,j)
        end if
        if(voln0(ig,j).gt.0._fp) then
           evol0(ig,j)=evol0(ig,j)/voln0(ig,j)
           ovol0(ig,j)=ovol0(ig,j)/voln0(ig,j)
        end if
     end do
  end do
  !
  !  **> all new dmc mar 1985 <**
  !  volume source thermal neutrals on 2d grid ***
  !
  inumi=nznbmr-nxtbzn
  !  loop over radial groups of beam dist fcn zones inside plasma
  jz2=lcm1
  iz2=nthzsm(lcm1)
  !
  do izg=1,inumi
     !  transp zone index limits
     jz1=jz2+1
     jz2=jz1+nxskpb-1
     !  beam dist fcn zone index limits
     iz1=iz2+1
     iz2=nthzsm(izg+lcm1)
     !
     !  loop over species
     do ig=1,ng
        !  add up vol sce neutral densities in transp zones encompassed by
        !  beam dist fcn zone row; add up 1d sources
        !
        !  also, compute wall neutral density on same grid but assume that
        !  it is evenly distribued poloidally (dmc dec 2001)
        !
        zomg0(ig)=0.0_fp
        zp0(ig)=0.0_fp
        zn0(ig)=0.0_fp
        zn0w(ig)=0.0_fp
        !  add up 1d (recombination) vol sce and total volume source
        zstot(ig)=0.0_fp
        zvol=0.0_fp
        do j=jz1,jz2
           zvol=zvol+dvolzc(j)
           zn0w(ig)=zn0w(ig)+dvolzc(j)*walln0(ig,j)
           zn0(ig)=zn0(ig)+dvolzc(j)*voln0(ig,j)
           zp0(ig)=zp0(ig)+dvolzc(j)*voln0(ig,j)*evol0(ig,j)
           zomg0(ig)=zomg0(ig)+dvolzc(j)*voln0(ig,j)*ovol0(ig,j)
           !  beam sorbn0, srbn02 already include *dv
           do isb=1,nsbeam
              zstot(ig)=zstot(ig)+sorbn0s(j,ig,isb)
           end do
        end do
        !
        !  fat zone average...
        walln0x(ig,jz1:jz2)=zn0w(ig)/zvol
        walln02(ig,iz1:iz2)=zn0w(ig)/zvol
        !
        !  loop thru beam dist fcn zones (in theta direction, this xi row)
        !
        !  mod dmc aug 1991 -- breakup loops (70,77) inside if-then-else clause
        !  to satisfy request from jet; avoid problem with ibm compiler
        !
        if(zn0(ig).eq.0.0_fp) then
           do iz=iz1,iz2
              evol02(ig,iz)=0.0_fp
              voln02(ig,iz)=0.0_fp
              ovol02(ig,iz)=0.0_fp
           end do
        else
           do iz=iz1,iz2
              evol02(ig,iz)=zp0(ig)/zn0(ig)
              ovol02(ig,iz)=zomg0(ig)/zn0(ig)
              !  2d neutral density presumed proportional to local source
              !  (neutral transport was evaluated but only in 1d
              if(zstot(ig).gt.0.0_fp) then
                 zsrbn0=0.0_fp
                 do isb=1,nsbeam
                    zsrbn0=zsrbn0+srbn02(iz,ig,isb)
                 end do
                 voln02(ig,iz)=zn0(ig)*(zsrbn0)/ &
                      (zstot(ig)*bmvol(iz))
              else
                 voln02(ig,iz)=zn0(ig)/zvol
              end if
           end do
        end if
        !  loop end ...
     end do                  ! ig
     !co-injected beam populations
     iz2_min=min(iz2,mimxbzf)
     do isb=1,nsbeam
        do ib=1,nbeam
           do ien=1,3
              ze_rel=0.0_fp
              zrms_rel=0.0_fp
              zvol=0.0_fp
              do iz=iz1,iz2_min
                 zvol=zvol+bmvol(iz)
                 ze_rel=ze_rel+bmvol(iz)* &
                      ebb_rel_co(iz,ien,ib,isb)
                 zrms_rel=zrms_rel+bmvol(iz)* &
                      rms_ebb_co(iz,ien,ib,isb)
              end do
              ebb_relx_co(jz1:jz2,ien,ib,isb)=ze_rel/zvol
              rms_ebbx_co(jz1:jz2,ien,ib,isb)=zrms_rel/zvol
           end do
        end do
     end do
     !counter-injected beam populations
     do isb=1,nsbeam
        do ib=1,nbeam
           do ien=1,3
              ze_rel=0.0_fp
              zrms_rel=0.0_fp
              zvol=0.0_fp
              do iz=iz1,iz2_min
                 zvol=zvol+bmvol(iz)
                 ze_rel=ze_rel+bmvol(iz)* &
                      ebb_rel_cnt(iz,ien,ib,isb)
                 zrms_rel=zrms_rel+bmvol(iz)* &
                      rms_ebb_cnt(iz,ien,ib,isb)
              end do
              ebb_relx_cnt(jz1:jz2,ien,ib,isb)=ze_rel/zvol
              rms_ebbx_cnt(jz1:jz2,ien,ib,isb)=zrms_rel/zvol
           end do
        end do
     end do
     !density
     do isb=1,nsbeam
        do ino=1,2
           zbnb=0.0_fp
           zvol=0.0_fp
           do iz=iz1,iz2_min
              zvol=zvol+bmvol(iz)
              zbnb=zbnb+bmvol(iz)* bnb_ien2(iz,ino,isb)
           end do
           bnb_ienx(jz1:jz2,ino,isb)=zbnb/zvol
        end do
     end do

     !fusioun product populations
     do isb=1,nsbeam
        do ib=1,nbeam
           do ien=1,3
              ze_rel=0.0_fp
              zrms_rel=0.0_fp
              zvol=0.0_fp
              do iz=iz1,iz2_min
                 zvol=zvol+bmvol(iz)
                 ze_rel=ze_rel+bmvol(iz)* &
                      ebb_rel_fn(iz,ien,ib,isb)
                 zrms_rel=zrms_rel+bmvol(iz)* &
                      rms_ebb_fn(iz,ien,ib,isb)
              end do
              ebb_relx_fn(jz1:jz2,ien,ib,isb)=ze_rel/zvol
              rms_ebbx_fn(jz1:jz2,ien,ib,isb)=zrms_rel/zvol
           end do
        end do
        zbnb=0.0_fp
        zvol=0.0_fp
        do iz=iz1,iz2_min
           zvol=zvol+bmvol(iz)
           zbnb=zbnb+bmvol(iz)*bnb_fn(iz,isb) !density by fusion product species
        end do
        bnb_fnx(jz1:jz2,isb)=zbnb/zvol
     end do
  end do                     ! fat zones

  !
  !--------------------------------------------------------------------
  !
  !  update the bilinear inverse map (r,y)-->(xi,theta) for flr, etc.
  !
  call nubeam_xstraln_init(nzp1,id_mcgrid,nonlin,ierr)
  iersum=iersum+abs(ierr)
  !
  call eq_gfnum('R',    id_r)
  call eq_gfnum('Z',    id_z)
  !
  call eq_gfnum('BPHI', id_bphi)
  call eq_gfnum('BR',   id_br)
  call eq_gfnum('BZ',   id_bz)
  if(min(id_bphi,id_br,id_bz).eq.0) then
     write(nonlin,*) '??parms: xplasma id missing for br bz bphi'
     iersum=iersum+1
  end if
  !
  call eq_gfnum('PSI',id_psi)
  if(id_psi.eq.0) then
     write(nonlin,*) '??parms: xplasma id missing for psi'
     iersum=iersum+1
  end if
  !
  !     --- define phi_orb profile on r,z grid ---
  !
  if(id_epot.gt.0) then
     call eqm_rzfun(id_epot, 1.e6_fp*raxis, 102, 0.0_fp, ierr)
     ! no decay from boundary
     if (ierr.ne.0) then
        iersum=iersum+1
        write(nonlin,*) &
             '?xpgload:  eqm_rzfun(id_phi,...) returned ierr=',ierr
     end if
  end if
  !
  !------------------------------------------
  if(nrip.eq.2) then
     call eq_gfnum('LOG_TF_RIPPLE',idi)
     if(idi.gt.0) id_ripple=idi
     if(id_ripple.eq.0) then
        iersum=iersum+1
        write(nonlin,*) &
             ' ?nubeam/parms:  nrip=2 set but no tf ripple data found.'
     end if

     if(ncoils_2.gt.0) then
        call eq_gfnum('LOG_TF_RIPPLE2',idi)
        if(idi.gt.0) id_ripple2=idi
        if(id_ripple2.eq.0) then
           iersum=iersum+1
           write(nonlin,*) &
                ' ?nubeam/parms:  nrip=2 & ncoils_2 > 0 but:'// &
                ' no tf ripple (2) data found.'
        end if
     end if
  end if
  !------------------------------------------
  !  store field at midplane; find max field around each flux surface...
  !
  call bparms(2,ierr)
  iersum=iersum+abs(ierr)
  !
  orbrzv_bnw_lrf=0.0_fp
  nlmcrfgcon=.true.
  mcrf_on=nlmcrf_on
  if(orbrzv_option .eq. 1) then
     ! initialize splines for guiding center rz-orbit integrator
     orbrzv_rfdv=0.0_fp
     call orbrz_init_spline_coeffs

     if(orbrzv_rf_option .eq. 1) then
        if(mcrf_on) then

           ! initialize updated rf-field renormalizer
           if(imcrf_first.eq.0) then
              orbrzv_toric_ufrnm(1:mibs) = orbrzv_toric_frnm
           end if
           imcrf_first=nbsteps

           ! get rfxqlo info: wavefields from toric...
           ! these may need to be mpi broadcast...

           call mcrf_toric_pack2(ierr)
           iersum=iersum+abs(ierr)
           call orbrz_set_rf_info2 !mcrf_on could be forced to be .false.

#ifndef __MPI
           !#ifdef __DEBUG
           if(mcrf_on) then
              rfdbg_pdiag = 0.0_fp

              zfile = 'prfh_'//trim(runid)//'.txt'

              orbrzv_out = 11111
              open(unit=orbrzv_out,file=trim(zfile),status="old", &
                   action="write",position="append",iostat = iopstat)

              if(iopstat > 0) then
                 print *, " %mcrf : prf-his.txt is being created"
                 open(unit=orbrzv_out,file=trim(zfile),status="new", &
                      action="write",position="append",iostat = iopstat1)
                 if(iopstat1 > 0) then
                    print *, " ?mcrf : prf-his.txt open failed"
                    call bad_exit('prf-his.txt new open fail')
                 end if
              end if

              write(unit=orbrzv_out,fmt=*) ' '
              write(unit=orbrzv_out,fmt=*) '######################'
              write(unit=orbrzv_out,fmt=*) '# time :', tbm1, tbm2
              write(unit=orbrzv_out,fmt=*) '# mcrf_on :', mcrf_on
           end if
           !     #endif
#endif
        else
           !perhaps there is not enough rf power, see trfpi/nubeam_component_step
           !  nlmcrf_on=.false.
           orbrzv_prftot = 0.0_fp
           orbrzv_prf = 0.0_fp
           orbrzv_toric_ufrnm(1:mibs) = 1.0_fp
           orbrzv_toric_prftot = 0.0_fp
           orbrzv_toric_prf = 0.0_fp
        end if
     end if
  end if
  iersum=iersum+abs(ierr)
  !
  !-----------------------------------
  !  mod dmc mar 2008 -- broadcast difb & velb to separate arrays for
  !  beam ions and fusion product ions
  !
  if(nmdifb.eq.1) then

     !  control what kind of ions diffuse anomalously:
     !  nkdifb=1 -- beam ions  nkdifb=2 -- fusion product ions
     !  nkdifb=3 -- both kinds
     !  and... all ions for x<x_qlim0

     if(x_qlim0.eq.0.0_fp) x_qlim0=-1.0_fp

     do j=lcentr,lep1
        if((nkdifb.eq.1).or.(nkdifb.eq.3).or. &
             (xi(j,1).le.x_qlim0)) then
           difb_bi(j)=difb(j)
           velb_bi(j)=velb(j)
        end if
        if((nkdifb.eq.2).or.(nkdifb.eq.3).or. &
             (xi(j,1).le.x_qlim0)) then
           difb_fi(j)=difb(j)
           velb_fi(j)=velb(j)
        end if
     end do
  end if
  !-----------------------------------
  !
  ierr=iersum
  !
  return
end subroutine parms
!******************** end file parms.for ; group parms *****************
