subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
    nbnd, rbnd, zbnd, nrho, psrad, fpol, te, dne, zeff, qpsi, ibeam,   &
    p0mw, alphain, betain, dpdv, jcd, pabs, icd, err)
  use const_and_precisions, only: wp_
  use gray_params,          only: gray_parameters, gray_data, gray_results
  use gray_core,            only: gray_main
  use gray_plasma,          only: numeric_plasma

  implicit none

  ! subroutine arguments
  integer, intent(in) :: ijetto, mr, mz, nbnd, nrho, ibeam
  real(wp_), dimension(mr),    intent(in)  :: r
  real(wp_), dimension(mz),    intent(in)  :: z
  real(wp_), dimension(mr,mz), intent(in)  :: psin
  real(wp_),                   intent(in)  :: psia, rax, zax, p0mw, alphain, betain
  real(wp_), dimension(nbnd),  intent(in)  :: rbnd, zbnd
  real(wp_), dimension(nrho),  intent(in)  :: psrad, fpol, te, dne, zeff, qpsi
  real(wp_), dimension(nrho),  intent(out) :: dpdv, jcd
  real(wp_),                   intent(out) :: pabs, icd
  integer,                     intent(out) :: err

  ! local variables
  type(gray_parameters) :: params
  type(gray_data)       :: data
  type(numeric_plasma)  :: plasma
  type(gray_results)    :: res
  logical, save         :: firstcall = .true.

  ! Initialisation tasks for the first call only
  first_call: if (firstcall) then
    firstcall = .false.

    ! Read parameters from external file
    init_params: block
      use gray_params, only : read_gray_params
      call read_gray_params('gray.data', params, err)
      if (err /= 0) return

      ! Override some parameters
      params%misc%rwall         = r(1)
      params%misc%active_tables = [4, 7, 8, 9, 12, 17, 33, 70, 71, 55, 48]
      params%antenna%filenm     = 'graybeam.data'
      params%equilibrium%filenm = 'JETTO'
      params%equilibrium%iequil = ijetto + 1
      params%profiles%filenm    = 'JETTO'
      params%profiles%iprof     = 1
      params%output%ipec        = 1
      params%output%nrho        = nrho
    end block init_params

    ! Set a simple limiter following the boundary of the data grid
    simple_limiter: block
      use const_and_precisions, only : cm
      use types,                only : contour

      ! Avoid clipping out the launcher
      real(wp_) :: R_launch, R_max
      R_launch = norm2(params%antenna%pos(1:2)) * cm
      R_max = max(R_launch, R(mr))

      ! Convert to a closed contour
      data%equilibrium%limiter = contour( &
        Rmin=params%misc%rwall, Rmax=R_max, zmin=z(1), zmax=z(mz))
    end block simple_limiter

  end if first_call

  ! Set MHD equilibrium data
  init_equilibrium: block
    use equilibrium, only : set_equil_spline
    use types,       only : contour

    ! Copy argument arrays
    data%equilibrium%rv    = r
    data%equilibrium%zv    = z
    data%equilibrium%rax   = rax
    data%equilibrium%rvac  = rax
    data%equilibrium%zax   = zax
    data%equilibrium%psinr = psrad
    data%equilibrium%fpol  = fpol
    data%equilibrium%psia  = psia
    data%equilibrium%psin  = psin
    data%equilibrium%qpsi  = qpsi
    data%equilibrium%boundary = contour(rbnd, zbnd)

    ! Compute splines
    call set_equil_spline(params%equilibrium, data%equilibrium, err)
    if (err /= 0) return
  end block init_equilibrium

  ! Compute splines of kinetic profiles
  call plasma%init(params, psrad, te, dne, zeff, err)
  if (err /= 0) return

  ! Set wave launcher parameters
  init_antenna: block
    use beams, only: read_beam2

    ! Copy argument variables
    params%antenna%alpha = alphain
    params%antenna%beta  = betain
    params%antenna%power = p0mw

    ! Read beam description file
    call read_beam2(params%antenna, beamid=ibeam, err=err)
    if (err /= 0) return
  end block init_antenna

  ! Call main subroutine for the ibeam-th beam
  call gray_main(params, data, plasma, res, err, rhout=sqrt(psrad))

  write_debug_files: block
    integer :: i, err
    ! Write results to file
    do i = 1, size(res%tables%iter)
      associate (tbl => res%tables%iter(i)%ptr)
      call tbl%save(err, filepath='gray_'//trim(tbl%title)//'.txt')
      end associate
    end do
  end block write_debug_files

  ! Free memory
  free_memory: block
    use equilibrium,  only : unset_equil_spline

    ! Unset global variables of the `equilibrium` module
    call unset_equil_spline
  end block free_memory

  ! Copy over the results
  pabs = res%pabs
  icd  = res%icd
  dpdv = res%dpdv
  jcd  = res%jcd

end subroutine gray_jetto1beam