!   Fortran 77 interface to JETTO
      subroutine gray(ijetto, mr, mz, mrd, r, z, psin, psiax, psibnd,
     .    rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff,
     .    qsf, nbeam, powin, alphin, betain, dpdv, jcd, pec, icd, ierr)
      implicit none
!   input arguments
      integer ijetto, mr, mz, mrd, nbnd, nrho, nbeam
      real*8 r(mr), z(mz), psin(mrd,mz)
      real*8 psiax, psibnd, rax, zax
      real*8 rbnd(nbnd), zbnd(nbnd)
      real*8 psijet(nrho), f(nrho), qsf(nrho), te(nrho), dne(nrho)
      real*8 zeff(nrho)
      real*8 powin(nbeam), alphin(nbeam), betain(nbeam)
!   output arguments
      real*8 dpdv(nrho), jcd(nrho), pec, icd
!   gray_main output arguments
      real*8 dpdvloop(nrho), jcdloop(nrho), pecloop, icdloop
      integer ierr
!   local variables
      integer i,j
      real*8 tescal(nrho), dnescal(nrho)
      real*8 p0mw

! === input arguments ==================================================
!
!     ijetto  Equilibrium source (1 EFIT, 2 ESCO)
!             If IJETTO=2, then PSIN values are valid only inside
!             plasma boudary (PSIN=0 outside)
!     mr      Size of flux map grid in R direction
!     mz      Size of flux map grid in Z direction
!     mrd     Leading dimension of the psin(:,:) array, mrd>mr
!     r       R coordinates of flux map grid points [m]
!     z       Z coordinates of flux map grid points [m]
!     psin    Normalised poloidal flux psin=(psi-psiax)/(psibnd-psiax)
!             on the (R, Z) grid.
!     psiax   Poloidal flux on axis [Wb rad-1]
!     psibnd  Poloidal flux on boundary [Wb rad-1]
!     rax     R coordinate of magnetic axis [m]
!     zax     Z coordinate of magnetic axis [m]
!     nbnd    Number of points in plasma boundary contour
!     rbnd    R coordinates of plasma boundary contour [m]
!     zbnd    Z coordinates of plasma boundary contour [m]
!
!     nrho    Number of points in JETTO rho grid    -
!     psijet  Normalised poloidal flux on JETTO radial grid
!     f       Poloidal current stream function f=Bphi*R on JETTO
!             radial grid [T m]
!     te      Electron temperature on JETTO radial grid [eV]
!     dne     Electron density on JETTO radial grid [m-3]
!     zeff    Effective nuclear charge Zeff on JETTO radial grid
!     qsf     Safety factor on JETTO radial grid
!
!     nbeam   Total number of injected beams
!     powin   Input ECRH power array [W] (powin(i) =< 0 means i-th beam is unused)
!     alphin  Beams poloidal injection angles array [rad]
!     betain  Beams toroidal injection angles array [rad]
!
! === output arguments =================================================
!
!     dpdv    Absorbed EC power density on JETTO radial grid [W m-3]
!     jcd     EC driven flux averaged current density on JETTO
!             radial grid [A m-2]
!     pec     Total absorbed EC power [W]
!     icd     Total EC driven current [A]
!     ierr    Return code. IERR>0 on error
!             ierr = 90-93: error computing integrals for current drive
!             ierr = 94: absorption coefficient alpha < 0
!             ierr = 97: parallel comp. refract. idx N//>0.99 (warning)
!             ierr = 98: parallel comp. refract. idx N//>1.05
!
! === Note =============================================================
!
!   JETTO coordinate system assumes toroidal angle increasing CW
!   in GRAY toroidal angle increases CCW --> adapt signs on input data
!
!   f   is passed as -f
!   qsf is passed as -qsf
!
!   jcd is returned as -jcd
!   icd is returned as -icd
!
! ======================================================================

!   set output variables to 0
      do i=1,nrho
        dpdv(i) = 0.d0
        jcd(i) = 0.d0
      end do
      pec = 0.d0
      icd = 0.d0
      
      dnescal = dne * 1.d-19
      tescal = te * 1.d-3
!     loop over beams with power>0
      do j=1,nbeam
        if (powin(j).le.0.0d0) cycle
        p0mw=powin(j)*1.d-6

!       read j-th beam properties from file
!       and adjust alpha/beta if out of the allowed range

!       call main subroutine for the j-th beam
        call gray_jetto1beam(ijetto, mr, mz, r, z, psin(1:mr,:),
     .    psibnd-psiax, rax, zax, nbnd, rbnd, zbnd, nrho, psijet, -f,
     .    tescal, dnescal, zeff, -qsf, j, p0mw, alphin(j), 
     .    betain(j), dpdvloop, jcdloop, pecloop, icdloop, ierr)

!       add contribution of j-th beam to the total
!       adapting output data to JETTO convention on toroidal angle
        do i=1,nrho
          dpdv(i) = dpdv(i) + dpdvloop(i)*1.d6
          jcd(i)  = jcd(i)  - jcdloop(i)*1.d6
        end do
        pec = pec + pecloop*1.d6
        icd = icd - icdloop*1.d6

!     end of loop over beams with power>0
      end do

      return
      end