subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, &
    rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, &
    powin, dpdv, jcd, pec, icd, ierr)

  use itm_types, only : r8
  implicit none

  integer, intent(in) :: ijetto, mr, mz, nrho, nbnd
  real(r8), intent(in) :: r(mr), z(mz), psin(mr,mz)
  real(r8), intent(in) :: psiax, psibnd, rax, zax
  real(r8), intent(in), dimension(nbnd) :: rbnd, zbnd
  real(r8), intent(in), dimension(nrho) :: psijet, f, qsf, te, dne, zeff
  real(r8), intent(in) :: powin
  real(r8), intent(out), dimension(nrho) :: dpdv, jcd
  real(r8), intent(out) :: pec, icd
  integer, intent(out) :: ierr

  integer :: istop, iercom, igrad, iopmin, iowmin, index_rt, ipass
  real(r8) :: sox, p0mw, powrfl, taumn, taumx, pabstot, currtot
  real(r8) :: p0mw1, powtr, pabstott, currtott

  common/istop/istop
  common/ierr/iercom
  common/igrad/igrad
  common/iovmin/iopmin,iowmin
  common/mode/sox
  common/p0/p0mw
  common/powrfl/powrfl
  common/index_rt/index_rt
  common/taumnx/taumn,taumx,pabstot,currtot
  common/ipass/ipass
      
! read data plus initialization
  index_rt=1
  call prfile
  call paraminit
  call read_data(ijetto, mr, mz, r, z, psin, psiax, psibnd, rax, zax, &
    nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, powin)
  call vectinit
  if(igrad.eq.0) call ic_rt
  if(igrad.gt.0) call ic_gb
  if(iercom.gt.0) then
    ierr=iercom
    write(*,*) '    IERR = ', ierr
    return
  end if

! beam/ray propagation
  call gray_integration

! postprocessing
  call after_gray_integration

  pabstott=pabstot
  currtott=currtot
  powtr=p0mw-pabstot

  if (iowmin==2 .and. ipass>1) then
! second pass into plasma
    p0mw1=p0mw
    igrad=0
    
    index_rt=2
    p0mw=p0mw1*powrfl
    call prfile
    call vectinit2
    call paraminit
    call ic_rt2
    call gray_integration
    call after_gray_integration
    pabstott=pabstott+pabstot
    currtott=currtott+currtot
  
    index_rt=3
    sox=-sox
    p0mw=p0mw1*(1._r8-powrfl)
    call prfile
    call vectinit2
    call paraminit
    call ic_rt2
    call gray_integration
    call after_gray_integration
    pabstott=pabstott+pabstot
    currtott=currtott+currtot
  end if
  write(*,*)
  write(*,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.e3_r8
  pec=pabstott*1.e6_r8
  icd=currtott*1.e6_r8
! NOT YET COMPLETE! 
  dpdv=0._r8
  jcd=0._r8
! =================
  ierr=iercom
end subroutine gray_main