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


  ! 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(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 utils,                only : range2rect
      use limiter,              only : limiter_set_globals=>set_globals
      use const_and_precisions, only : cm

      ! 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 list of R,z
      call range2rect(xmin=params%misc%rwall, xmax=R_max, ymin=z(1), ymax=z(mz), &
                      x=data%equilibrium%rlim, y=data%equilibrium%zlim)

      ! Set the global variables of `limiter`
      call limiter_set_globals(data%equilibrium)
    end block simple_limiter

  end if first_call

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

    ! 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%rbnd  = rbnd
    data%equilibrium%zbnd  = zbnd
    data%equilibrium%psia  = psia
    data%equilibrium%psin  = psin
    data%equilibrium%qpsi  = qpsi

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

  ! Set plasma kinetic profiles
  init_profiles: block
    use coreprofiles, only : set_profiles_spline

    ! Copy argument arrays
    data%profiles%derad = dne
    data%profiles%terad = te
    data%profiles%zfc   = zeff

    ! Compute splines
    call set_profiles_spline(params%profiles, data%profiles, err)
    if (err /= 0) return
  end block init_profiles

  ! 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, 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
    use coreprofiles, only : unset_profiles_spline

    ! Unset global variables of the `equilibrium` module
    call unset_equil_spline

    ! Unset global variables of the `coreprofiles` module
    call unset_profiles_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