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 units,                only: set_active_units, close_units
  use gray_params,          only: gray_parameters, gray_data, gray_results
  use gray_core,            only: gray_main

  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(gray_results)    :: res
  logical, save         :: firstcall = .true.

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

    ! Activate some debug units
    set_units: block
      use units
      call set_active_units([         &
        ucenr, usumm, uprj0, uprj0+1, &
        uwbm, udisp, ubres, ucnt,     &
        uoutr, umaps, uprfin, uflx,   &
        upec])
    end block set_units

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

      ! Override some parameters
      params%misc%rwall         = r(1)
      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

      ! Set the global variables of `gray_params`
      call set_globals(params)
    end block init_params

    ! Set a simple limiter
    simple_limiter: block
      use reflections, only : range2rect
      use limiter,     only : limiter_set_globals=>set_globals

      real(wp_) :: r0m
      r0m = norm2(params%antenna%pos(1:2)) * 0.01_wp_
      call range2rect(params%misc%rwall, max(r0m, r(mr)), z(1), z(mz), &
                      data%equilibrium%rlim, 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))

  ! 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

  call close_units

end subroutine gray_jetto1beam