!   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)
      use units, only : ucenr,usumm,uprj0,uprj0+1,uwbm,udisp,ubres,
     .    ucnt,uoutr,ueq,uprfin,uflx,upec
      use gray_params, only : read_params
      use beams, only : read_beam2
      use graycore, only : gray_main
      implicit none
!   input arguments
      integer ijetto, mr, mz, 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
      real*8 rlim(5),zlim(5)
      logical firstcall=.true.
      save firstcall

! === 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
!
! ======================================================================

!   if first call read parameters from external file
      if (firstcall) then
        call read_params('gray.data',rtrp,hcdp,antp,eqp,rwallm,
     .    prfp,outp,uprm)
        antp%filenm='graybeam.data'
        eqp%filenm='JETTO'
        eqp%iequil=ijetto+1
        prfp%filenm='JETTO'
        firstcall=.false.
      end if

!   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
      
!     loop over beams with power>0
      do j=1,nbeam
        if (powin(j).gt.0.0d0) cycle

!       read j-th beam properties from file
!       and adjust alpha/beta if out of the allowed range
        alpha0=alphin(j)
        beta0=betain(j)
        p0mw=powin(j)*1.d-6
        call read_beam2(antp%filenm,j,alpha0,beta0,fghz,antp%iox,
     .                  x0,y0,z0,w1,w2,ri1,ri2,phiw,phir,ubeam)
        psipol0=antp%psi
        chipol0=antp%chi
        iox0=antp%iox

!       set simple limiter
        r0m=sqrt(x0**2+y0**2)*0.01d0
        call range2rect(rwallm,max(r0m,rv(mr)),zv(1),zv(mz),rlim,zlim)

!       call main subroutine for the j-th beam
        subroutine gray_main(r,z,psin(1:mr,:),psibnd-psiax,
     .      psijet,-f,-qsf,rax,rax,zax,rbnd,zbnd,eqp,
     .      psijet,te,dne,zeff,prfp,rlim,zlim,
     .      p0mw,fghz,alpha0,beta0,(/x0,y0,z0/),
     .      w1,w2,ri1,ri2,phiw,phir,iox0,psipol0,chipol0,
     .      dpdvloop,jcdloop,pecpool,icdloop,outp,rtrp,hcdp,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)
          jcd(i)  = jcd(i)  - jcdloop(i)
        end do
        pec = pec + pecloop
        icd = icd - icdloop

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

!     close output (debug) files
      close(ucenr,usumm,uprj0,uprj0+1,uwbm,udisp,ubres,ucnt,uoutr,
     .    ueq,uprfin,uflx,upec)

      return
      end