module graycore
  use const_and_precisions, only : wp_
  implicit none

contains

subroutine gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp,   &
            psrad,terad,derad,zfc,prfp,  rlim,zlim,                       &
            p0,fghz,alpha0,beta0,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir,iox0, &
            psipol0,chipol0,  dpdv,jcd,pec,icd,  outp,rtrp,hcdp,ierr)
  use coreprofiles, only : set_prfan, set_prfspl
  use gray_params, only : eqparam_type, prfparam_type, outparam_type, &
     rtrparam_type, hcdparam_type, set_codepar, igrad, iequil, iprof
  use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff
  use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl
  use beamdata, only : init_rtr, dealloc_beam, nstep, dst
  use reflections, only : set_lim
  implicit none
  type(eqparam_type), intent(in) :: eqp
  type(prfparam_type), intent(in) :: prfp
  type(outparam_type), intent(in) :: outp
  type(rtrparam_type), intent(in) :: rtrp
  type(hcdparam_type), intent(in) :: hcdp
  real(wp_), dimension(:), allocatable, intent(in) :: psrad, terad, derad, zfc
  real(wp_), dimension(:), allocatable, intent(in) :: rv, zv, psinr, fpol, qpsi
  real(wp_), dimension(:), allocatable, intent(in) :: rbnd, zbnd, rlim, zlim
  real(wp_), dimension(:,:), allocatable, intent(in) :: psin
  real(wp_), intent(in) :: psia, rvac, rax, zax
  integer, intent(in) :: iox0
  real(wp_), intent(in) :: p0, fghz, psipol0, chipol0
  real(wp_), intent(in) :: alpha0,beta0, x0,y0,z0, w1,w2, ri1,ri2, phiw,phir
  real(wp_), intent(out) :: pec,icd
  real(wp_), dimension(:), intent(out) :: dpdv,jcd
  integer, intent(out) :: ierr
  real(wp_), dimension(:), allocatable :: rhotn
  
  real(wp_) :: sox,ak0,bres,xgcn,anx0,any0,anz0
  integer :: i,iox,istop,index_rt

      integer :: istep
      real(wp_) :: st
      common/ss/st
      common/istep/istep

! ======= set environment BEGIN ======
  call set_codepar(eqp,prfp,outp,rtrp,hcdp)

  call set_lim(rlim,zlim)

  if(iequil<2) then
    call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3))
    call flux_average_an
  else
    call set_eqspl(rv,zv,psin, psia, psinr,fpol, eqp%ssplps,eqp%ssplf, rvac,  &
       rax,zax, rbnd,zbnd, eqp%ixp)
! compute rho_pol/rho_tor mapping
    allocate(rhotn(size(qpsi)))
    call setqphi_num(psinr,abs(qpsi),abs(psia),rhotn)
    call set_rhospl(sqrt(psinr),rhotn)
    deallocate(rhotn)
! compute flux surface averaged quantities
    call flux_average ! requires frhotor for dadrhot,dvdrhot
! print  psi surface for q=1.5 and q=2
    call surfq(psinr,qpsi,size(qpsi),1.5_wp_)
    call surfq(psinr,qpsi,size(qpsi),2.0_wp_)
  end if

  if(iprof==0) then
    call set_prfan(terad,derad,zfc)
  else
    call set_prfspl(psrad, terad, derad, zfc, prfp%sspld, prfp%psnbnd)
  end if

  call xgygcoeff(fghz,ak0,bres,xgcn)
  call launchangles2n(alpha0,beta0,x0,y0,z0,anx0,any0,anz0)
  call init_rtr(rtrp)

! ======= set environment  END  ======

! ======= pre-proc prints BEGIN ======
! print Btot=Bres
! print ne, Te, q, Jphi versus psi, rhop, rhot
  if (iequil==1) then
    call bres_anal(bres)
    call print_prof_an
  else
    call bfield_res(rv,zv,size(rv),size(zv),bres)
    call print_prof
  end if
! ======= pre-proc prints  END  ======

! ======= main loop BEGIN ======
  iox=iox0
  sox=-1.0_wp_
  if(iox.eq.2) sox=1.0_wp_
  index_rt=1
  call prfile
  call paraminit
  call vectinit
  if(igrad==0) then
    call ic_rt(x0,y0,z0,anx0,any0,anz0,ak0,xgcn,bres, &
       w1,w2,ri1,ri2,phiw,phir,sox,psipol0,chipol0)
  else
    call ic_gb(x0,y0,z0,anx0,any0,anz0,ak0,xgcn,bres, &
       w1,w2,ri1,ri2,phiw,phir,sox,psipol0,chipol0)
  end if
!  if(ierr==0) return
!  beam/ray propagation
! st0=0.0_wp_
! if(index_rt.gt.1) st0=strfl11
  do i=1,nstep
    istep=i
    st=i*dst !+st0
!  advance one step
    call rkint4(sox,xgcn,bres)
!  calculations after one step:
    call after_onestep(p0,sox,ak0,xgcn,bres,i,istop,ierr)
    if(istop.eq.1) exit
  end do
!  postprocessing
  call after_gray_integration(pec,icd,dpdv,jcd)

! ======= main loop  END  ======

! ======= free memory BEGIN ======
  call dealloc_beam
!  call unset...
! ======= free memory  END  ======
end subroutine gray

end module graycore